%%
%% This program is written for plotting the calculated convection location from
%%  CMIT MIX NH and SH electric potentials found in driver.f90 and calccloc.f90
%% BA Emery, December 2011 from /aim3/emery/cmit/plot_cloc.m

clear all; close all

fileout = sprintf('imf_out');
IMFOUT = fopen(fileout,'w');

filebad = sprintf('bad_crad');
BADOUT = fopen(filebad,'w');
fprintf(BADOUT, 'yr daynm mo da hr mn Bx,y,z Vsw Dsw \n');
fprintf(BADOUT,'SthenN:cp x/n x,lat,mlt n,lat,mlt \n');
fprintf(BADOUT,'second line: SthenN:cp,theta0,arad,offc,offa,dskofc,dskofa,phid,phin \n'); 

% nmph = number of minutes per hour (30x2min WHI, 12x5min AGU Dec06)
% dmin = number of minutes between CMIT files 
% nhpd = number of hours per day (usually 24, but is 29 for Dec 14-15 2006)

% Set year and ifcheck
year1 = 2008

% ifflag=1 if read new files with ifbad
% ifflag=2 set center to -0.4 deg always, and set crad=12 if both hems are bad or to arad-4.0 deg
% Found 14 nogoods (both hem) with ifflag=0, while 16 were read with ifflag=1 for Dec06
% Found 1477 nogoods (both hem) with ifflag=0, while 1526 were read with ifflag=1 for CR07
% Found 1372 nogoods (both hem) with ifflag=0, while 1432 were read with ifflag=1 for WHI08
% TOTAL 3014-61=2953 bad both hems of 39918-61=39857 times (7.6% or 7.4% instead of 12% bad with 1 or both hem - got 4.5% back!)
ifflag = 2

% Set daynumbers for study
if (year1==2006)
% Dec 14-15 2006 5-min MIX files
dayn0 = 348;
dayn1 = 348;
%dayn2 = 349;	% for 1.94.1 version
dayn2 = 350;	% for 1.94.2 version
x1 = 14.;
x2 = 15.25;
x1 = dayn1;
x2 = dayn1+1.25;	% for 1.94.1 version
x2 = dayn2;		% for 1.94.2 version
cpx = 550;
radx = 21;
bmag = 20;
dmin = 5;
nhpd = 29;		% for 1.94.1 version
nhpd = 24;
end

if (year1==2007)
% CR2060 Aug 14 (1UT start rev to 0UT) - Sep 10 2007 (226-253) 2-min MIX files
dayn0 = 226;
dayn1 = 226
dayn2 = 254
x1 = dayn1;
x2 = dayn2;
cpx = 225;
radx = 19;
bmag = 10;
dmin = 2;
nhpd = 24;
end

if (year1==2008)
% WHI Mar 20 - Apr 15 2008 (80-106) 2-min MIX files
dayn0 = 80;
dayn1 = 89;
dayn1 = 98;
dayn2 = 107;
dayn1 = 98
dayn2 = dayn1 + 9
dayn1 = 80
dayn2 = 107
x1 = dayn1;
x2 = dayn2;
cpx = 225;
radx = 17;
bmag = 10;
dmin = 2;
nhpd = 24;
end


nmph = 60/dmin

% Set if plot S+N (for dynamo, ifboth=0) or separate (ifboth=1)
ifboth=1

of18x = 14;
of18x = 9;
of24x = 7;
if ifboth==0
 of18x = 3;
 of24x = 4.5;
end


% Get number for 5 min OMNI data to skip over
nd0 = (dayn1-1)*12*24;


% Open 15-min lagged every 5 min (t=0 uses t=-20 to -5min) .nc IMF OMNI file
if (year1==2006)
 Fname = sprintf('/hao/aim3/tgcm/data/imf_OMNI_2006001-2006365.nc')
end
if (year1==2007)
 Fname = sprintf('/hao/aim3/tgcm/data/imf_OMNI_2007001-2007365.nc')
end
if (year1==2008)
 Fname = sprintf('/hao/aim3/tgcm/data/imf_OMNI_2008001-2008366.nc')
end

% open .nc file Fname
    nc=netcdf(Fname,'nowrite');

    date=nc{'date'}(:);        % (year+daynum+fracday)
    bzyr=nc{'bz'}(:);
    byyr=nc{'by'}(:);
    bxyr=nc{'bx'}(:);
    vswyr=nc{'swvel'}(:);
    dswyr=nc{'swden'}(:);
    close(nc) 
    clear nc

n = 0;
nyn = 0;
nyp = 0;
% Get 1 extra day for IMF
%for nd=dayn1:dayn2-1
for nd=dayn1:dayn2
for nh=1:24
%for nh=1:24+5
% Read every 5 min for OMNI
 for nl=1:12
% Find medians for 0-14 UT 14 Dec 2006:CPn,s=74.1,60.6kV; off24n,s=0.50,0.97deg;
%  off18n,s=+1.96,-1.96deg; cradn,s=10.94,11.96deg; By=-1.45; Bz=-1.02
%  ie, off18s~+By and off18n~-By
% Find medians for 08080-085 Mar 20-25 2008:CPn,s=33.2,32.6kV; off24n,s=0.50,0.97deg;
%  off18n,s=-3.88,+1.95deg; cradn,s=13.87,13.89deg; By=1.08; Bz=0.94
% Find medians for 08080-085 Mar 20-25 2008 #2:CPn,s=33.2,32.6kV; off24n,s=-0.97,0.97deg;
%  off18n,s=-1.94,+0.97deg; cradn,s=13.47,13.03deg (w Nan=0 so low?); By=1.08; Bz=0.94
%  Neg By: -0.88 vs -0.86; off24 -0.97N,+1.91S; off18 +0.97N,-0.97S   default 08080-85
%  Pos By: +1.79 vs +1.80; off24 -0.97N,+0.02S; off18 -3.91N,+2.94S   default
%  ie, off18s~+By and off18n~-By
%  if nh==14 & nl==1
  if nd== 86 & nh==1 & nl==1
   medby = median(by)
   medbz = median(bz)
   medbyneg = median(byneg)
   medbypos = median(bypos)
  end
 n = n + 1;
%  Assign 5-min IMF data
 daynfrac(n) = date(nd0+n)-year1*1000.;
 bx(n) = bxyr(nd0+n);
 by(n) = byyr(nd0+n);
 bz(n) = bzyr(nd0+n);
 vsw(n) = vswyr(nd0+n);
 dsw(n) = dswyr(nd0+n);
% Sort by By sign
  if (by(n)<0) 
    nyn = nyn + 1;
    byneg(nyn) = by(n);
  else
    nyp = nyp + 1;
    bypos(nyp) = by(n);
  end
 end
end
end	% do nd
% Had to add 1 at end for interpolation with 2min data at end until went 1 day more
% Recalc medians for whole periodoutcmit_cr2060_2007b
   medby = median(by)
   medbz = median(bz)
   medbyneg = median(byneg)
   medbypos = median(bypos)

% Open file of outputs 
if (year1==2006)
% filloc = sprintf('outcmit_agu2006_dec14-15_2006_0UT-5UT')
% filloc = sprintf('outcmit_agu2006')
% filloc = sprintf('outcmit_agu2006b')
% filloc = sprintf('outcmit_agu2006_test')
 filloc = sprintf('outcmit_decb2006_test')	% Aug 2012 version (SH noon=midnight)
%filloc = sprintf('outcmit_dec2006_testc')	% Sep 2012 version w arad from eflux and HP
 filloc = sprintf('outcmit_dec2006_testd')	% Sep 2012 w Rc=Ra+4, Ra=max diff>4deg 
end
if (year1==2007)
% filloc = sprintf('outcmit_cr2060_2007b')
%filloc = sprintf('outcmit_aug-sep_2007_test')
%filloc = sprintf('outcmit_2007cr2060_all_testc')	% Sep 2012 version w arad from eflux
 filloc = sprintf('outcmit_2007cr2060_all_testd') % Sep 2012 w Rc=Ra+4, Ra=max diff>4deg
end
if (year1==2008)
% filloc = sprintf('outcmit_whi2008_mar20-apr15')
% filloc = sprintf('outcmit_whi2008_mar-apr2')
% filloc = sprintf('outcmit_whimar08')
% filloc = sprintf('outcmit_whi2008')
% filloc = sprintf('outcmit_whi2008b')
%filloc = sprintf('outcmit_mar-apr_2008_test')
%filloc = sprintf('outcmit_whi2008_all_testc')	% Sep 2012 version w arad from eflux
 filloc = sprintf('outcmit_whi2008_all_testd') % Sep 2012 w Rc=Ra+4, Ra=max diff>4deg
end
locin = fopen(filloc,'r');
% Skip first 2 lines
linloc = fgetl(locin)
linloc = fgetl(locin)
 for nd=dayn0:dayn1-1
 for nh=1:24
 for nl=1:30
 linloc1 = fgetl(locin);
 linloc2 = fgetl(locin);
 end
 end
 end	% for nd

%  write (11,"('yr daynm mo da hr mn SthenN:cp x,lat,mlt n,lat,mlt')")
%  write (11,"('second line: SthenN:cp,theta0,arad,offc,offa,dskofc,dskofa,phid,phin')") 
%  write (11,"(a20,2x,2f7.1,2f7.2,f7.1,2f7.2,2x,2f7.1,2f7.2,f7.1,2f7.2)") datehrmn, cps,potsx,latsx,mltsx,potsn,latsn,mltsn, cpn,potnx,latnx,mltnx,potnn,latnn,mltnn
%  write (11,"(22x,2f7.1,16f7.2)")  weictpoten,theta0*rtd,rrad*rtd,offc*rtd,offa*rtd,dskofc*rtd,dskofa*rtd,phid*rtd/15.+12.,phin*rtd/15.+12.

nbadab = 0;
nbadas = 0;
nbadan = 0;
nbadcb = 0;
n = 0;
ntru = 0;
nrlo = 0;
nlo = 0;
nyn = 0;
nyp = 0;
nog = 0;
nogs = 0;
nogn = 0;
nog1 = 0;
nog60 = 0;
nog50 = 0;
nog40 = 0;
nog30 = 0;
for nd=dayn1:dayn2-1
nd,n
%for nh=1:24
%for nh=1:24+5
 for nh=1:nhpd
% Read every 5 min for Dec 14-15, 2006
% Read every 2 min for WHI 2008 Mar 20-Apr 15
 for nl=1:nmph
% Find medians for 0-14 UT 14 Dec 2006:CPn,s=74.1,60.6kV; off24n,s=0.50,0.97deg;
%  off18n,s=+1.96,-1.96deg; cradn,s=10.94,11.96deg; By=-1.45; Bz=-1.02
%  ie, off18s~+By and off18n~-By
% if nh==14 & nl==1
  if nd==86 & nh==1 & nl==1
 medby2 = median(by2)
 medcpn = median(cpn)
 medcps = median(cps)
 medoff24n = median(offcn)
 medoff24s = median(offcs)
 medoff18n = median(dskofcn)
 medoff18s = median(dskofcs)
 fp = isnan(cradn);
 cradn(fp) = 0;
 fp = isnan(crads);
 crads(fp) = 0;
 medcradn = median(cradn)
 medcrads = median(crads)
 medby2neg = median(byneg2)
 medby2pos = median(bypos2)
 medoff24nneg = median(offcnneg)
 medoff24sneg = median(offcsneg)
 medoff24npos = median(offcnpos)
 medoff24spos = median(offcspos)
 medoff18nneg = median(dskofcnneg)
 medoff18sneg = median(dskofcsneg)
 medoff18npos = median(dskofcnpos)
 medoff18spos = median(dskofcspos)
n,nm,np,fracdaymo(n),daynfrac(nm),daynfrac(np),frac1,frac2,by(nm),by(np),by2(n)
 end

% Read file from driver.f90 and calccloc.f90
 linloc1 = fgetl(locin);
 in1 = sscanf(linloc1,'%f');
 linloc2 = fgetl(locin);
 in2 = sscanf(linloc2,'%f');
 ntru = ntru + 1;

 fracdaynum = in1(2)+in1(5)/24.+in1(6)/(24.*60.);
if (year1==2006)
 bxint = bx(ntru);
 byint = by(ntru);
 bzint = bz(ntru);
 vswint = vsw(ntru);
 dswint = dsw(ntru);
else
% Interpolate 5min IMF data to 2min CMIT times (1=1, 6=3, 11=5, 16=7 ...)
nm = fix((ntru-1)*0.40+1);
np = nm + 1;
% daynfrac-dayn1+x1 = fracdaymo
 frac1 = (fracdaynum-daynfrac(nm)+dayn1-x1)/(daynfrac(np)-daynfrac(nm));
 frac2 = (daynfrac(np)-dayn1+x1-fracdaynum)/(daynfrac(np)-daynfrac(nm));
 bxint = frac2*bx(nm) + frac1*bx(np);
 byint = frac2*by(nm) + frac1*by(np);
 bzint = frac2*bz(nm) + frac1*bz(np);
 vswint = frac2*vsw(nm) + frac1*vsw(np);
 dswint = frac2*dsw(nm) + frac1*dsw(np);
end
% Print out IMF for use in driver.f90
fprintf(IMFOUT, '%d %3.3d %2.2d %2.2d %6.1f %6.1f %6.1f %6.1f %6.1f \n',in1(1),in1(2),in1(5),in1(6),bxint,byint,bzint,vswint,dswint);


%nogood = 0;
% Check problems w mltx,n (reverse convection? often radii~20deg)
% if (in1(10)>12. | in1(13)<12. | in1(17)>12. | in1(20)<12.)
%  nogood = 1;
% end
% Check for problems with large crad (crad+5=crit(1) in range 15-30, or crad 10-25 (limit 5-25))
% if (in2(3)>45. | in2(3)<5. | in2(4)>45. | in2(4)<5.)
% Found no neg crad except for bad default values in corrupt (8-18T13-24) or missing files
%  Found crad<5 only for 07230 1324ut (which was missing anyway and had fake values!)
%  nogood = 1;
% end

 nogood = 0;
 nogoods = 0;
 nogoodn = 0;

% For new files with ifbad flags, don't check for problems
 if (ifflag==0) 

% Check problems w mltx,n (reverse convection? often radii~20deg)
% Check for problems with large crad (crad+5=crit(1) in range 15-30, or crad 10-25 (limit 5-25))
%   Can get large radii when dskofc very large (no +/-3h min/max at 18 or 6 MLT)
%  Check to see if only 1 hem is off
  if (in1(10)>12. | in1(13)<12. | in2(3)>25.)
   nogoods = 1;
   nogs = nogs + 1; 
  end
  if (in1(17)>12. | in1(20)<12. | in2(4)>25.)
   nogoodn = 1;
   nogn = nogn + 1;    
  end
  if (nogoods==1 & nogoodn==1)
   nogood = 1;
  end
   if (nogoods+nogoodn==1)
    nog1 = nog1 + 1;
% Revise bad hemisphere only (should also revise phin,d to be defaults w byimf)
    if (nogoods==1) 
   in2(3) = in2(4);
   in2(5) = in2(6);
   in2(7) = -in2(8);
   in2(9) = -in2(10);
   in2(11) = -in2(12);
   in2(13) = -in2(14); 
%   crads = in2(4)
%   offcs = -in2(8)
%   dskofcs = -in2(12)  
    else
   in2(4) = in2(3);
   in2(6) = in2(5);
   in2(8) = -in2(7);
   in2(10) = -in2(9);
   in2(12) = -in2(11);
   in2(14) = -in2(13);
%   cradn = in2(3)
%   offcn = -in2(7)
%   dskofcn = -in2(11) 
    end
   end

 else
% For ifflag=1, do nothing - but plot out as is
  if (ifflag==2)
% Set center to zero always for CMIT - WRONG! - pole position sets dynamo zero!
   for i=7:10
    in2(i) = -0.4;
   end
% Set convection center at pole (not towards dawn or dusk)
    in2(11) = 0.;
    in2(12) = 0.;
% Set auroral center to -2.5 deg (2.5 deg towards dawn)
    in2(13) = -2.5;
    in2(14) = -2.5;
% Check if both hems are bad and set crad=12 or crad=arad-4
    if (in2(25)>0 & in2(26)>0)
    nbadcb = nbadcb + 1;
%  If good arad, set crad from arad
     if (in2(5)>=14 & in2(6)>=14)
      in2(3) = in2(5)-4.0;
      in2(4) = in2(6)-4.0;
     else
      if (in2(5)<14 & in2(6)>=14)
       nbadas = nbadas + 1;
       in2(3) = in2(4);
       in2(5) = in2(6);
      end
      if (in2(5)>=14 & in2(6)<14)
       nbadan = nbadan + 1;
       in2(4) = in2(3);
       in2(6) = in2(5);
      end
      if (in2(5)<14 & in2(6)<14)
       nogood = 1;
       nbadab = nbadab + 1;
       in2(3) = 12.;
       in2(5) = in2(3) + 4.;
       in2(4) = 12.;
       in2(6) = in2(4) + 4.;
      end
    end		% else
   end		% ifbad>0 both
  end		% ifflag=2
 end		% else

% TEMP See all!
% nogood = 0;


if (nogood==0)
 n = n + 1;
% Assign interpolated IMF data
 bx2(n) = bxint;
 by2(n) = byint;
 bz2(n) = bzint;
 vsw2(n) = vswint;
 dsw2(n) = dswint;
% Assign read in data
 fracdaymo(n) = in1(4)+in1(5)/24.+in1(6)/(24.*60.);  % day of month
 fracdaymo(n) = in1(2)+in1(5)/24.+in1(6)/(24.*60.);  % daynumber of year
 cps(n) = in1(7);
 latsx(n) = in1(9);
 latsn(n) = in1(12);
 mltsx(n) = in1(10);
 mltsn(n) = in1(13);
 cps2(n) = in2(1);
 crads(n) = in2(3);
 arads(n) = in2(5);
 offcs(n) = in2(7);
 offas(n) = in2(9);
 dskofcs(n) = in2(11);
 dskofas(n) = in2(13);
 phids(n) = in2(15);
 phins(n) = in2(17);
 cpn(n) = in1(14);
 latnx(n) = in1(16);
 latnn(n) = in1(19);
 mltnx(n) = in1(17);
 mltnn(n) = in1(20);
 cpn2(n) = in2(2);
 cradn(n) = in2(4);
 aradn(n) = in2(6);
 offcn(n) = in2(8);
 offan(n) = in2(10);
 dskofcn(n) = in2(12);
 dskofan(n) = in2(14);
 phidn(n) = in2(16);
 phinn(n) = in2(18);
 amcradn(n) = aradn(n)-cradn(n);
 amcrads(n) = arads(n)-crads(n);
% 9/12 add HP from eflux
 hps(n) = in1(21);
 hpn(n) = in1(22);
 hpb(n) = 0.5*(in1(21)+in1(22));
 aradb(n) = 0.5*(in2(5)+in2(6));
% average CP and crad
 cradb(n) = 0.5*(in2(3)+in2(4));
 amcradb(n) = aradb(n)-cradb(n);
 cpb(n) = 0.5*(in2(1)+in2(2));
% Sort by n/x
 if (in1(8)/in1(9)<-3. | in1(15)/in1(18)<-3.)
  nrlo = nrlo + 1;
  cprlo(nrlo) = 0.5*(cps(n)+cpn(n));
  cradrlo(nrlo) = 0.5*(crads(n)+cradn(n));
 end
% Sort by CP
  if (cps(n)<50. | cpn(n)<50.)
   nlo = nlo + 1;
   cplo(nlo) = 0.5*(cps(n)+cpn(n));
   cradlo(nlo) = 0.5*(crads(n)+cradn(n));
  end
% Sort by By sign
  if (by2(n)<0 & nogood==0) 
    nyn = nyn + 1;
    nneg(nyn) = n;
    byneg2(nyn) = by2(n);
    dskofcsneg(nyn) = dskofcs(n);
    dskofcnneg(nyn) = dskofcn(n);
    offcnneg(nyn) = offcn(n);
    offcsneg(nyn) = offcs(n);
  end
  if (by2(n)>0 & nogood==0) 
    nyp = nyp + 1;
    npos(nyp) = n;
    bypos2(nyp) = by2(n);
    dskofcspos(nyp) = dskofcs(n);
    dskofcnpos(nyp) = dskofcn(n);
    offcnpos(nyp) = offcn(n);
    offcspos(nyp) = offcs(n);
  end
else
% nogood = 1 so save CP to see what default situations are!
 nog = nog + 1;
 cpbad(nog) = 0.5*(in2(1)+in2(2));
 daynbad(nog) = fracdaynum;
 off24bad(nog) = 0.;
 off18bad(nog) = 0.;
%  Use formula for crad(CP) from Dec06 both
 cradbad(nog) = 9.5925 + 0.0226*cpbad(nog);
% Use formula for crad(CP) from WHI08 (08080-08106)
 cradbad(nog) = 12.221 + 0.0083*cpbad(nog);
% Use formula for crad(CP) from CR2060 (07226-07353)
 cradbad(nog) = 11.3446 + 0.0136*cpbad(nog);
% Use formula for crad(CP) from 2006-7-8b for n/x<-2 (11.319 at 100kV)
 cradbad(nog) = 9.459 + 0.0186*cpbad(nog);
% Use median 11.6016*0.0152*18.25 = 11.879 = 11.9 or 12.0 from ave of medians in histogram bins
 cradbad(nog) = 12.0;
 if (cpbad(nog)>30)
   nog30 = nog30 + 1;
 end
 if (cpbad(nog)>40)
   nog40 = nog40 + 1;
 end
 if (cpbad(nog)>50)
   nog50 = nog50 + 1;
 end
 if (cpbad(nog)>60)
   nog60 = nog60 + 1;
fprintf(BADOUT, '%d BAD %dkV ave S+N \n',nog60,fix(cpbad(nog)));
 end
fprintf(BADOUT, '%d %d %d %d %d %d ',in1(1),in1(2),in1(3),in1(4),in1(5),in1(6));
fprintf(BADOUT, '%6.1f',bxint,byint,bzint,vswint,dswint);
fprintf(BADOUT,'\n');
fprintf(BADOUT, '%8.2f',in1(7),in1(8)/in1(11),in1(8),in1(9),in1(10),in1(11),in1(12),in1(13),in1(14),in1(15)/in1(18),in1(15),in1(16),in1(17),in1(18),in1(19),in1(20));
fprintf(BADOUT,'\n');
fprintf(BADOUT, '%8.2f',in2);
fprintf(BADOUT,'\n');
end	% if nogood==0 for all
 end
end 
end	% do nd
n,fracdaymo(n)
 if (nog > 0)
  medbadcp = median(cpbad)
 end
% Recalc medians for whole period
 medby2 = median(by2)
 medcpn = median(cpn)
 medcps = median(cps)
 medhpn = median(hpn)
 medhps = median(hps)
 medoff24n = median(offcn)
 medoff24s = median(offcs)
 medoff18n = median(dskofcn)
 medoff18s = median(dskofcs)
 fp = isnan(cradn);
 cradn(fp) = 0;
 fp = isnan(crads);
 crads(fp) = 0;
 medcradn = median(cradn)
 medcrads = median(crads)
 medaradn = median(aradn)
 medarads = median(arads)
 medamcradn = median(amcradn)
 medamcrads = median(amcrads)
 medby2neg = median(byneg2)
 medby2pos = median(bypos2)
 medoff24nneg = median(offcnneg)
 medoff24sneg = median(offcsneg)
 medoff24npos = median(offcnpos)
 medoff24spos = median(offcspos)
 medoff18nneg = median(dskofcnneg)
 medoff18sneg = median(dskofcsneg)
 medoff18npos = median(dskofcnpos)
 medoff18spos = median(dskofcspos) 
 medoff24nps = median(0.5*(offcn+offcs))
 medoff18nps = median(0.5*(dskofcn+dskofcs))
 medamcradnps = median(0.5*(amcradn+amcrads))
 avoff24nps = mean(0.5*(offcn+offcs))
 avoff18nps = mean(0.5*(dskofcn+dskofcs))

nbadcb,nbadab,nbadas,nbadan

figure
zero = [0.,0.];
xtot = [x1,x2];

fsize = 12;
% position [left bottom width height]
% Set size to be 8.5x11 inches
set(gcf,'position',[142    99   647   822])
H = gcf;
clf

subplot('position',[0.10,0.82,0.87,0.12])	% 1/4 of 1/2
hold on;
if (year1==2006)
title('Convection Locations Etc for Dec 14-15, 2006 CMIT MIX','FontSize',15,'FontWeight','bold');
end
if (year1==2007)
title('Convection Locations Etc for CR2060 2007 CMIT MIX','FontSize',15,'FontWeight','bold');
end
if (year1==2008)
title('Convection Locations Etc for WHI 2008 CMIT MIX','FontSize',15,'FontWeight','bold');
end
if ifboth==1
if (nog > 0)
 plot (daynbad,cpbad,'+k');
end
plot (fracdaymo,cps,'-b');
plot (fracdaymo,cpn,'-r');
plot (fracdaymo,hps*10.,':b');
plot (fracdaymo,hpn*10.,':r');
text (x1+0.05*(x2-x1),0.90*cpx,'NH','Color','r','FontSize',fsize,'FontWeight','bold');
text (x1+0.10*(x2-x1),0.90*cpx,'SH','Color','b','FontSize',fsize,'FontWeight','bold');
text (x1+0.15*(x2-x1),0.90*cpx,'Bad','Color','k','FontSize',fsize,'FontWeight','bold');
text (x1+0.05*(x2-x1),0.70*cpx,'...HP(GW)*10','Color','k','FontSize',fsize,'FontWeight','bold');
else
if (nog > 0)
 plot (daynbad,cpbad,'+r');
end
plot (fracdaymo,0.5*(cps+cpn),'-k');
plot (fracdaymo,0.5*(hps+hpn)*10.,'-m');
text (x1+0.05*(x2-x1),0.70*cpx,'...HP(GW)*10','Color','m','FontSize',fsize,'FontWeight','bold');
text (x1+0.02*(x2-x1),0.90*cpx,'av(N+S)','Color','k','FontSize',fsize,'FontWeight','bold');
text (x1+0.15*(x2-x1),0.90*cpx,'Bad','Color','r','FontSize',fsize,'FontWeight','bold');
end
ylabel('CP kV','FontSize',fsize,'FontWeight','bold');
xlabel(' ','FontSize',fsize,'FontWeight','bold');
grid;
%set (gca,'xlim',[x1,x2],'ylim',[0,450],'XTick',[14 14.125 14.25 14.375 14.5 14.625 14.75 14.875 15. 15.125 15.25],'XTickLabel',{' ';'3';'6';'9';'12';'15';'18';'21';'0';'3';' ';},'FontSize',fsize,'box','on');
set (gca,'xlim',[x1,x2],'ylim',[0,cpx],'XTickLabel',{' ';},'FontSize',fsize,'box','on');
subplot('position',[0.10,0.70,0.87,0.12])	% 2/4 of 1/2
hold on;
if ifboth==1
if (nog > 0)
 plot (daynbad,off24bad,'+k');
end
plot (fracdaymo,offcs,'-b');
%plot (fracdaymo,offas,'--b');
plot (fracdaymo,offcn,'-r');
%plot (fracdaymo,offan,'--r');
else
plot (fracdaymo,0.5*(offcs+offcn),'-k');
end
text (x1+0.05*(x2-x1),6.0,'midnight','Color','k','FontSize',fsize,'FontWeight','bold');
text (x1+0.05*(x2-x1),-6.0,'noon','Color','k','FontSize',fsize,'FontWeight','bold');
if (ifflag==1)
if (year1==2006)
text (x1+0.80*(x2-x1),6.0,'median -0.69','Color','k','FontSize',fsize,'FontWeight','bold');
end
if (year1==2007)
text (x1+0.80*(x2-x1),6.0,'median -0.01','Color','k','FontSize',fsize,'FontWeight','bold');
end
if (year1==2008)
text (x1+0.80*(x2-x1),6.0,'median -0.47','Color','k','FontSize',fsize,'FontWeight','bold');
end
end	% ifflag=1
ylabel('off24 deg','FontSize',fsize,'FontWeight','bold');
xlabel(' ','FontSize',fsize,'FontWeight','bold');
grid;
%set (gca,'xlim',[x1,x2],'ylim',[0,7],'XTick',[14 14.125 14.25 14.375 14.5 14.625 14.75 14.875 15. 15.125 15.25],'XTickLabel',{' ';'3';'6';'9';'12';'15';'18';'21';'0';'3';' ';},'FontSize',fsize,'box','on');
set (gca,'xlim',[x1,x2],'ylim',[-of24x,of24x],'XTickLabel',{' ';},'FontSize',fsize,'box','on');
subplot('position',[0.10,0.58,0.87,0.12])	% 3/4 of 1/2
hold on;
if ifboth==1
if (nog > 0)
 plot (daynbad,off18bad,'+k');
end
plot (fracdaymo,dskofcs,'-b');
%plot (fracdaymo,dskofas,'--b');
plot (fracdaymo,dskofcn,'-r');
%plot (fracdaymo,dskofan,'--r');
else
plot (fracdaymo,0.5*(dskofcs+dskofcn),'-k');
end
plot (xtot,zero,'--k');
text (x1+0.05*(x2-x1),0.8*of18x,'18 MLT','Color','k','FontSize',fsize,'FontWeight','bold');
text (x1+0.05*(x2-x1),-0.8*of18x,'06 MLT','Color','k','FontSize',fsize,'FontWeight','bold');
ylabel('off18 deg','FontSize',fsize,'FontWeight','bold');
xlabel(' ','FontSize',fsize,'FontWeight','bold');
grid;
%set (gca,'xlim',[x1,x2],'ylim',[-14,14],'XTick',[14 14.125 14.25 14.375 14.5 14.625 14.75 14.875 15. 15.125 15.25],'XTickLabel',{' ';'3';'6';'9';'12';'15';'18';'21';'0';'3';' ';},'FontSize',fsize,'box','on');
set (gca,'xlim',[x1,x2],'ylim',[-of18x,of18x],'XTickLabel',{' ';},'FontSize',fsize,'box','on');
subplot('position',[0.10,0.46,0.87,0.12])	% 4/4 of 1/2
hold on;
if ifboth==1
if (nog > 0)
 plot (daynbad,cradbad,'+k');
end
plot (fracdaymo,crads,'-b');
%plot (fracdaymo,arads,'--b');
plot (fracdaymo,cradn,'-r');
%plot (fracdaymo,aradn,'--r');
plot (fracdaymo,arads,':b');
plot (fracdaymo,aradn,':r');
text (x1+0.05*(x2-x1),8.5,'...Ra','Color','k','FontSize',fsize,'FontWeight','bold');
else
 plot (fracdaymo,0.5*(crads+cradn),'-k');
 plot (fracdaymo,0.5*(arads+aradn),'-m');
 text (x1+0.05*(x2-x1),8.5,'...Ra','Color','m','FontSize',fsize,'FontWeight','bold');
end
ylabel('radius deg','FontSize',fsize,'FontWeight','bold');
xlabel(' ','FontSize',fsize,'FontWeight','bold');
grid;
%set (gca,'xlim',[x1,x2],'ylim',[5,30],'XTick',[14 14.125 14.25 14.375 14.5 14.625 14.75 14.875 15. 15.125 15.25],'XTickLabel',{' ';'3';'6';'9';'12';'15';'18';'21';'0';'3';' ';},'FontSize',fsize,'box','on');
set (gca,'xlim',[x1,x2],'ylim',[7,radx],'XTickLabel',{' ';},'FontSize',fsize,'box','on');
subplot('position',[0.10,0.34,0.87,0.12])	% 5
hold on;
% 9/12 Replace x n lat loc with arad-crad
if ifboth==1
plot (fracdaymo,amcrads,'-b');
plot (fracdaymo,amcradn,'-r');
if (year1==2006)
% text (x1+0.80*(x2-x1),9.,'median 4.76','Color','k','FontSize',fsize,'FontWeight','bold'); % for testc
text (x1+0.50*(x2-x1),9.,'Ra=max(RaS,N) if diff>4deg; median 4.80','Color','k','FontSize',fsize,'FontWeight','bold'); % for testd
end
if (year1==2007)
%text (x1+0.80*(x2-x1),9.,'median 4.32','Color','k','FontSize',fsize,'FontWeight','bold');
text (x1+0.50*(x2-x1),9.,'Ra=max(RaS,N) if diff>4deg; median 4.38','Color','k','FontSize',fsize,'FontWeight','bold'); % for testd
end
if (year1==2008)
%text (x1+0.80*(x2-x1),9.,'median 3.65','Color','k','FontSize',fsize,'FontWeight','bold');	% for testc
text (x1+0.50*(x2-x1),9.,'Ra=max(RaS,N) if diff>4deg; median 4.00','Color','k','FontSize',fsize,'FontWeight','bold'); % for testd
end
else
 plot (fracdaymo,0.5*(amcrads+amcradn),'-k');
end
ylabel('c-arad deg','FontSize',fsize,'FontWeight','bold');
xlabel(' ','FontSize',fsize,'FontWeight','bold');
grid;
set (gca,'xlim',[x1,x2],'ylim',[-1 10],'XTickLabel',{' ';},'FontSize',fsize,'box','on');
%if ifboth==1
%plot (fracdaymo,latsx,'-b');
%plot (fracdaymo,latsn,'--b');
%plot (fracdaymo,latnx,'-r');
%plot (fracdaymo,latnn,'--r');
%else
%plot (fracdaymo,0.5*(latsx+latnx),'-r');
%plot (fracdaymo,0.5*(latsn+latnn),'-b');
%end
%ylabel('x n-- lat deg','FontSize',fsize,'FontWeight','bold');
%xlabel(' ','FontSize',fsize,'FontWeight','bold');
%text (x1+0.02*(x2-x1),65.,'---min epot','Color','k','FontSize',fsize,'FontWeight','bold');
%grid;
%set (gca,'xlim',[x1,x2],'ylim',[60,90],'XTick',[14 14.125 14.25 14.375 14.5 14.625 14.75 14.875 15. 15.125 15.25],'XTickLabel',{' ';'3';'6';'9';'12';'15';'18';'21';'0';'3';' ';},'FontSize',fsize,'box','on');
%set (gca,'xlim',[x1,x2],'ylim',[60,90],'XTickLabel',{' ';},'FontSize',fsize,'box','on');
subplot('position',[0.10,0.22,0.87,0.12])	% 6
hold on;
if ifboth==1
plot (fracdaymo,mltsx,'-b');
plot (fracdaymo,mltsn,'--b');
plot (fracdaymo,mltnx,'-r');
plot (fracdaymo,mltnn,'--r');
else
plot (fracdaymo,0.5*(mltsx+mltnx),'-k');
plot (fracdaymo,0.5*(mltsn+mltnn),'-k');
end
ylabel('x n-- mlt h','FontSize',fsize,'FontWeight','bold');
xlabel(' ','FontSize',fsize,'FontWeight','bold');
grid;
%set (gca,'xlim',[x1,x2],'ylim',[0,24],'XTick',[14 14.125 14.25 14.375 14.5 14.625 14.75 14.875 15. 15.125 15.25],'XTickLabel',{' ';'3';'6';'9';'12';'15';'18';'21';'0';'3';' ';},'FontSize',fsize,'box','on');
set (gca,'xlim',[x1,x2],'ylim',[0,24],'XTickLabel',{' ';},'FontSize',fsize,'box','on');
subplot('position',[0.10,0.10,0.87,0.12])	% 7
hold on;
% fracdaymo = dayfracn - 348 + 14
plot (daynfrac-dayn1+x1,by,'-b');
plot (daynfrac-dayn1+x1,bz,'-r');
plot (xtot,zero,'--k');
ylabel('15m IMF nT','FontSize',fsize,'FontWeight','bold');
xlabel('Daynumber','FontSize',fsize,'FontWeight','bold');
text (x1+0.02*(x2-x1),0.8*bmag,'Bz','Color','r','FontSize',fsize,'FontWeight','bold');
text (x1+0.10*(x2-x1),0.8*bmag,'By','Color','b','FontSize',fsize,'FontWeight','bold');
grid;
%set (gca,'xlim',[x1,x2],'ylim',[-19,19],'XTick',[14 14.125 14.25 14.375 14.5 14.625 14.75 14.875 15. 15.125 15.25],'XTickLabel',{'0';'3';'6';'9';'12';'15';'18';'21';'0';'3';'6';},'FontSize',fsize,'box','on');
set (gca,'xlim',[x1,x2],'ylim',[-bmag,bmag],'FontSize',fsize,'box','on');
hold off;

byt = [-10,10];
%Calculate linear fit for off18 and By: off18 = b*By+c, b=p2(1), const=c=p2(2)
%  Limits good from -6 to +4 off18
% 80-88:  off18s=+0.19+0.645By; off18n=-0.63-0.725*By;  off24s=0.63-0.091*By; off24n=-0.64+0.028*By
% 89-97:  off18s=-0.20+0.590By; off18n=-0.21-0.779*By
% 98-106: off18s=+0.25+0.630By; off18n=-0.93-0.830*By
% 80-106: off18s=+0.08+0.624By; off18n=-0.61-0.782*By; medoff18s=0, medoff18n=-0.97
p2s18t=polyfit(by2,dskofcs,1)
p2s18neg=polyfit(byneg2,dskofcsneg,1)
p2s18pos=polyfit(bypos2,dskofcspos,1)
p2n18t=polyfit(by2,dskofcn,1)
p2n18neg=polyfit(byneg2,dskofcnneg,1)
p2n18pos=polyfit(bypos2,dskofcnpos,1)
ys18t = polyval(p2s18t,byt);
yn18t = polyval(p2n18t,byt);
ys18neg = polyval(p2s18neg,byt);
yn18neg = polyval(p2n18neg,byt);
ys18pos = polyval(p2s18pos,byt);
yn18pos = polyval(p2n18pos,byt);

p2s24t=polyfit(by2,offcs,1)
p2n24t=polyfit(by2,offcn,1)
ys24t = polyval(p2s24t,byt);
yn24t = polyval(p2n24t,byt);


if (year1==2006)
tit_daynum = sprintf('Dec 14-15, 2006 CMIT MIX');
end
if (year1==2007)
tit_daynum = sprintf('CMIT MIX for CR2060 2007 %d-%d',dayn1,dayn2-1);
end
if (year1==2008)
tit_daynum = sprintf('CMIT MIX for WHI 2008 %d-%d',dayn1,dayn2-1);
end

figure
y18=11;
y24=6;


% position [left bottom width height]
% Set size to be 8.5x11 inches
set(gcf,'position',[142    99   647   822])
H = gcf;

 subplot('Position',[0.07,0.54,0.35,0.42])     % upper left in 4
 hold on;
title(tit_daynum,'FontSize',15,'FontWeight','bold');
 plot(by2,dskofcs,'kx');
 med = medoff18s;
 medt = [med,med];
 plot (byt,medt,'--k');
 plot (byt,ys18t,'-k');
 plot (byt,ys18neg,'--b');
 plot (byt,ys18pos,'--r');
 ylabel('SH Off18MLT','FontSize',fsize,'Color','k');
 xlabel('By (nT)','FontSize',fsize,'Color','k');
grid;
set (gca,'xlim',[-10 10],'ylim',[-y18 y18],'FontSize',fsize,'box','on');

 subplot('Position',[0.555,0.54,0.35,0.42])     % upper right in 4
 hold on;
 plot(by2,offcs,'kx');
 med = medoff24s;
 medt = [med,med];
 plot (byt,medt,'--k');
 plot (byt,ys24t,'-k');
 ylabel('SH Off24MLT','FontSize',fsize,'Color','k');
 xlabel('By (nT)','FontSize',fsize,'Color','k');
grid;
set (gca,'xlim',[-10 10],'ylim',[-y24 y24],'FontSize',fsize,'box','on');

 subplot('Position',[0.07,0.04,0.35,0.42])     % bottom left in 4
 hold on;
 plot(by2,dskofcn,'kx');
 med = medoff18n;
 medt = [med,med];
 plot (byt,medt,'--k');
 plot (byt,yn18t,'-k');
 plot (byt,yn18neg,'--b');
 plot (byt,yn18pos,'--r');
 ylabel('NH Off18MLT','FontSize',fsize,'Color','k');
 xlabel('By (nT)','FontSize',fsize,'Color','k');
grid;
set (gca,'xlim',[-10 10],'ylim',[-y18 y18],'FontSize',fsize,'box','on');

 subplot('Position',[0.555,0.04,0.35,0.42])     % bottom right in 4
 hold on;
 plot(by2,offcn,'kx');
 med = medoff24n;
 medt = [med,med];
 plot (byt,medt,'--k');
 plot (byt,yn24t,'-k');
 ylabel('NH Off24MLT','FontSize',fsize,'Color','k');
 xlabel('By (nT)','FontSize',fsize,'Color','k');
grid;
set (gca,'xlim',[-10 10],'ylim',[-y24 y24],'FontSize',fsize,'box','on');
 
hold off;

% 80-88:  off24s=0.59-0.104*Bz; off24n=-0.57+0.082*Bz
% 89-97:  off24s=0.64-0.050*Bz; off24n=-0.75+0.061*Bz
% 98-106: off24s=0.96-0.038*Bz; off24n=-0.99+0.096*Bz
% 80-106: off24s=0.70-0.078*Bz; off24n=-0.75+0.073*Bz; medoff24s=+0.97 medoff24n=-0.97
bzt = [-10 10];
q2s18t=polyfit(bz2,dskofcs,1)
q2n18t=polyfit(bz2,dskofcn,1)
zs18t = polyval(q2s18t,bzt);
zn18t = polyval(q2n18t,bzt);

q2s24t=polyfit(bz2,offcs,1)
q2n24t=polyfit(bz2,offcn,1)
zs24t = polyval(q2s24t,bzt);
zn24t = polyval(q2n24t,bzt);

figure
y18=11;
y24=6;

% position [left bottom width height]
% Set size to be 8.5x11 inches
set(gcf,'position',[142    99   647   822])
H = gcf;

 subplot('Position',[0.07,0.54,0.35,0.42])     % upper left in 4
 hold on;
title(tit_daynum,'FontSize',15,'FontWeight','bold');
 plot(bz2,dskofcs,'kx');
 plot (bzt,zs18t,'-k');
 ylabel('SH Off18MLT','FontSize',fsize,'Color','k');
 xlabel('Bz (nT)','FontSize',fsize,'Color','k');
grid;
set (gca,'xlim',[-10 10],'ylim',[-y18 y18],'FontSize',fsize,'box','on');

 subplot('Position',[0.555,0.54,0.35,0.42])     % upper right in 4
 hold on;
 plot(bz2,offcs,'kx');
 plot (bzt,zs24t,'-k');
 ylabel('SH Off24MLT','FontSize',fsize,'Color','k');
 xlabel('Bz (nT)','FontSize',fsize,'Color','k');
grid;
set (gca,'xlim',[-10 10],'ylim',[-y24 y24],'FontSize',fsize,'box','on');

 subplot('Position',[0.07,0.04,0.35,0.42])     % bottom left in 4
 hold on;
 plot(bz2,dskofcn,'kx');
 plot (bzt,zn18t,'-k');
 ylabel('NH Off18MLT','FontSize',fsize,'Color','k');
 xlabel('Bz (nT)','FontSize',fsize,'Color','k');
grid;
set (gca,'xlim',[-10 10],'ylim',[-y18 y18],'FontSize',fsize,'box','on');

 subplot('Position',[0.555,0.04,0.35,0.42])     % bottom right in 4
 hold on;
 plot(bz2,offcn,'kx');
 plot (bzt,zn24t,'-k');
 ylabel('NH Off24MLT','FontSize',fsize,'Color','k');
 xlabel('Bz (nT)','FontSize',fsize,'Color','k');
grid;
set (gca,'xlim',[-10 10],'ylim',[-y24 y24],'FontSize',fsize,'box','on');
 
hold off;

% 80-88:  crads=12.1+0.0070*CPS, cradn=12.1+0.0080*CPN; medCP<50 12.47deg
% 89-97:  crads=12.2+0.0105*CPS, cradn=11.9+0.0136*CPN; medCP<50 12.29deg
% 98-106: crads=12.4+0.0096*CPS, cradn=12.4+0.0092*CPN; medCP<50 12.78deg
% 80-106: crads=12.3+0.0081*CPS, cradn=12.2+0.0095*CPN; medCP<50 12.48deg

cpt = [0 cpx];
rscpt=polyfit(cps,crads,1)
rncpt=polyfit(cpn,cradn,1)
xscpt = polyval(rscpt,cpt);
xncpt = polyval(rncpt,cpt);
ravcp=polyfit(cplo,cradlo,1)
xcplo = polyval(ravcp,cpt);
medradcplo = median(cradlo)
ravcpr=polyfit(cprlo,cradrlo,1)
xcprlo = polyval(ravcpr,cpt);
medradcprlo = median(cradrlo)
rbcpt=polyfit(cpb,cradb,1)
xbcpt = polyval(rbcpt,cpt);

figure

% position [left bottom width height]
% Set size to be 8.5x11 inches
set(gcf,'position',[142    99   647   822])
H = gcf;
 subplot('Position',[0.07,0.06,0.85,0.40]) 
 hold on;
title(tit_daynum,'FontSize',15,'FontWeight','bold');
 plot (cps,crads,'xb');
 plot (cpn,cradn,'+r');
 plot (cpt,xscpt,'-b');
 plot (cpt,xncpt,'-r');
 plot (cpt,xbcpt,'-k');
 plot (cpt,xcprlo,'-m');
text (0.86*cpx,9.5,'NH','Color','r','FontSize',fsize,'FontWeight','bold');
text (0.93*cpx,9.5,'SH','Color','b','FontSize',fsize,'FontWeight','bold');
 ylabel ('Convection Radius (deg)','FontSize',fsize,'Color','k');
 xlabel ('CP (kV)','FontSize',fsize,'Color','k');
grid;
set (gca,'xlim',[0 cpx],'ylim',[8 radx],'FontSize',fsize,'box','on');
 subplot('Position',[0.07,0.54,0.85,0.40]) 
 hold on;
title(tit_daynum,'FontSize',15,'FontWeight','bold');
 plot (cplo,cradlo,'xk');
 plot (cpt,xcplo,'-k');
 plot (cprlo,cradrlo,'xm');
 plot (cpt,xcprlo,'-m');
 ylabel ('Convection Radius (deg)','FontSize',fsize,'Color','k');
 xlabel ('Low CP (kV)','FontSize',fsize,'Color','k');
text (0.86*cpx,9.5,'cpn/x<-3','Color','m','FontSize',fsize,'FontWeight','bold');
grid;
set (gca,'xlim',[0 cpx],'ylim',[8 radx],'FontSize',fsize,'box','on');

hold off;

nog,nog1,nogs,nogn,nog30,nog40,nog50,nog60



