%%
%% 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, February 2012 from /aim3/emery/cmit/plot_crad_cp.m (based on plot_cloc.m)

clear all; close all

% Set ifcheck
% ifflag=1 if read new files with ifbad
% ifflag=2 set center to 0 always, and set crad=12 if both hems are bad
ifflag = 1

cpx = 450;
radx = 22;

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

% Open output from driver+calccloc.f90 and imf from plot_cloc.m
% 1851 bad 1 hem (4.6%) + 2947 both hem bad (7.4%, 36x>60kV, 16.85kV median) = 4798 times bad (12.0% 3645SH+4100NH) 
% num = 39917
% filloc = sprintf('outcmit_2006-7-8b')
% fileimf = sprintf('imf_out_2006-7-8b')
% Missing is 31 so 39948-31=39917 ntru
% num = 39948
% filloc = sprintf('outcmit_2006-7-8_ifbad')
% fileimf = sprintf('imf_out_2006-7-8_ifbad')
 num = 40176
% filloc = sprintf('outcmit_2006b-7-8_ifbad')	% Aug 2012 version (SH noon=midnight)
 filloc = sprintf('outcmit_2006-7-8c_ifbadc')	% Sep 2012 version
 fileimf = sprintf('imf_out_2006b-7-8_ifbad')

IMF = fopen(fileimf,'r');
locin = fopen(filloc,'r');
% Skip first 2 lines
linloc = fgetl(locin)
linloc = fgetl(locin)

n = 0;
ntru = 0;
nmiss = 0;
nrlo = 0;
nrlo2 = 0;
nzpyb = 0;
nzybad = 0;
nrlobad = 0;
nrlo2bad = 0;
nzybad0 = 0;
nrlobad0 = 0;
nrlo2bad0 = 0;
nlo = 0;
nyn = 0;
nyp = 0;
nog = 0;
nog1 = 0;
nogs = 0;
nogn = 0;
nog60 = 0;
nog50 = 0;
nog40 = 0;
nog30 = 0;
n10 = 0;
n15 = 0;
n20 = 0;
n25 = 0;
n30 = 0;
n40 = 0;
n60 = 0;
n80 = 0;
n90 = 0;

for nt=1:num

% Read IMF data
 linimf = fgetl(IMF);
 inimf = sscanf(linimf,'%f');
 bxint = inimf(5);
 byint = inimf(6);
 bzint = inimf(7);
 vswint = inimf(8);
 dswint = inimf(9);

% 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.);

 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
   for i=7:12
    in2(i) = 0;
   end
% Set auroral center to -2.5 deg (2.5 deg towards dawn)
    in2(13) = -2.5;
    in2(14) = -2.5;
  end	% ifflag=2
% Check if both hems are bad and set crad=12 if ifflag=2
    if (in2(25)>0 & in2(26)>0)
     if (ifflag==2) 
     in2(3) = 12.;
     in2(4) = 12.;
     end
     nogood = 1;
     nogoods = 1;
     nogoodn = 1;
    else
%  Both are good or only 1 is good
      if (in2(25)>0)
        nogoods = 1;
	nogs = nogs + 1;
	nog1 = nog1 + 1;
      end
      if (in2(26)>0)
        nogoodn = 1;
	nogn = nogn + 1;
	nog1 = nog1 + 1;
      end
    end

 end	% ifflag=0,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);
% average CP and crad
 cradb(n) = 0.5*(in2(3)+in2(4));
 aradb(n) = 0.5*(in2(5)+in2(6));
 amcradb(n) = aradb(n)-cradb(n);
 cpb(n) = 0.5*(in2(1)+in2(2));
 offcb(n) = 0.5*(in2(7)+in2(8));
% Save arad-crad (ifarad=1 arad found from eflux instead of arad=crad+2)
 amcrads(n) = arads(n)-crads(n);
 amcradn(n) = aradn(n)-cradn(n);
% 9/12 latest files with HP S,N
 hpb(n) = 0.5*(in1(21)+in1(22));	% only in files Sep 2012 and later
 hps(n) = in1(21);
 hpn(n) = in1(22);

% Sort by 1 bad hem
   if (nogoods>0 | nogoodn>0)
    cp1bad(nog1) = cpb(n);
    crad1bad(nog1) = cradb(n);
   end
% Sort by n/x (minpot/maxpot)
 if (in1(8)/in1(11)<-2. | in1(15)/in1(18)<-2.)
  nrlo2 = nrlo2 + 1;
  cprlo2(nrlo2) = cpb(n);
  cradrlo2(nrlo2) = cradb(n);
 end
 if (in1(8)/in1(11)<-3. | in1(15)/in1(18)<-3.)
  nrlo = nrlo + 1;
  cprlo(nrlo) = cpb(n);
  cradrlo(nrlo) = cradb(n);
 end
% Sort by x/n (maxpot/minpot)
 if (in1(11)/in1(8)<-2. | in1(18)/in1(15)<-2.)
  nrlo2 = nrlo2 + 1;
  cprlo2(nrlo2) = cpb(n);
  cradrlo2(nrlo2) = cradb(n);
 end
 if (in1(11)/in1(8)<-3. | in1(18)/in1(15)<-3.)
  nrlo = nrlo + 1;
  cprlo(nrlo) = cpb(n);
  cradrlo(nrlo) = cradb(n);
 end
% Sort by Bz>0 and |By|>1.2Bz
  if (bz2(n)>0. & abs(by2(n))>1.2*bz2(n))
   nzpyb = nzpyb + 1;
   cpzpyb(nzpyb) = cpb(n);
   cradzpyb(nzpyb) = cradb(n);
  end
% Sort by CP
  if (cpb(n)<10.)
   n10 = n10 + 1;
   hp10(n10) = hpb(n);
   cp10(n10) = cpb(n);
   crad10(n10) = cradb(n);
   arad10(n10) = aradb(n);
   amcrad10(n10) = amcradb(n);
   offc10(n10) = offcb(n);
  end
  if (cpb(n)>=10. & cpb(n)<15.)
   n15 = n15 + 1;
   hp15(n15) = hpb(n);
   cp15(n15) = cpb(n);
   crad15(n15) = cradb(n);
   arad15(n15) = aradb(n);
   amcrad15(n15) = amcradb(n);
   offc15(n15) = offcb(n);
  end
  if (cpb(n)>=15. & cpb(n)<20.)
   n20 = n20 + 1;
   hp20(n20) = hpb(n);
   cp20(n20) = cpb(n);
   crad20(n20) = cradb(n);
   arad20(n20) = aradb(n);
   amcrad20(n20) = amcradb(n);
   offc20(n20) = offcb(n);
  end
  if (cpb(n)>=20. & cpb(n)<25.)
   n25 = n25 + 1;
   hp25(n25) = hpb(n);
   cp25(n25) = cpb(n);
   crad25(n25) = cradb(n);
   arad25(n25) = aradb(n);
   amcrad25(n25) = amcradb(n);
   offc25(n25) = offcb(n);
  end
  if (cpb(n)>=25. & cpb(n)<30.)
   n30 = n30 + 1;
   hp30(n30) = hpb(n);
   cp30(n30) = cpb(n);
   crad30(n30) = cradb(n);
   arad30(n30) = aradb(n);
   amcrad30(n30) = amcradb(n);
   offc30(n30) = offcb(n);
  end
  if (cpb(n)>=30. & cpb(n)<40.)
   n40 = n40 + 1;
   hp40(n40) = hpb(n);
   cp40(n40) = cpb(n);
   crad40(n40) = cradb(n);
   arad40(n40) = aradb(n);
   amcrad40(n40) = amcradb(n);
   offc40(n40) = offcb(n);
  end
  if (cpb(n)>=40. & cpb(n)<60.)
   n60 = n60 + 1;
   hp60(n60) = hpb(n);
   cp60(n60) = cpb(n);
   crad60(n60) = cradb(n);
   arad60(n60) = aradb(n);
   amcrad60(n60) = amcradb(n);
   offc60(n60) = offcb(n);
  end
  if (cpb(n)>=60. & cpb(n)<100.)
   n80 = n80 + 1;
   hp80(n80) = hpb(n);
   cp80(n80) = cpb(n);
   crad80(n80) = cradb(n);
   arad80(n80) = aradb(n);
   amcrad80(n80) = amcradb(n);
   offc80(n80) = offcb(n);
  end
  if (cpb(n)>=100.)
   n90 = n90 + 1;
   hp90(n90) = hpb(n);
   cp90(n90) = cpb(n);
   crad90(n90) = cradb(n);
   arad90(n90) = aradb(n);
   amcrad90(n90) = amcradb(n);
   offc90(n90) = offcb(n);
  end
else
% nogood = 1 so save CP to see what default situations are!
% ifflag=1 and both are 4, this is for missing data so skip and revise ntru etc
 ifskip = 0;
 if (ifflag==1)
  if (in2(25)==4 & in2(26)==4)
   ifskip = 1;
   nmiss = nmiss + 1;
   ntru = ntru - 1;
  end
 end
if (ifskip==0)
 nog = nog + 1;
 cpbad(nog) = 0.5*(in2(1)+in2(2));
 daynbad(nog) = fracdaynum;
 off24bad(nog) = 0.1;	% 9/12 Change from 0 since cannot use 0 offc+dskofc in tiegcm
 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);
% Sort by Bz>0 and |By|>1.2Bz
  if (bzint>0. & abs(byint)>1.2*bzint)
   nzybad0 = nzybad0 + 1;
  end
% Sort by n/x or x/n
 if (in1(8)/in1(11)<-3. | in1(15)/in1(18)<-3. | in1(11)/in1(8)<-3. | in1(18)/in1(15)<-3.)
  nrlobad0 = nrlobad0 + 1;
 end
 if (in1(8)/in1(11)<-2. | in1(15)/in1(18)<-2. | in1(11)/in1(8)<-2. | in1(18)/in1(15)<-2.)
  nrlo2bad0 = nrlo2bad0 + 1;
 end

 if (cpbad(nog)>30)
   nog30 = nog30 + 1;
% Sort by Bz>0 and |By|>1.2Bz
  if (bzint>0. & abs(byint)>1.2*bzint)
   nzybad = nzybad + 1;
  end
% Sort by n/x or x/n
 if (in1(8)/in1(11)<-3. | in1(15)/in1(18)<-3. | in1(11)/in1(8)<-3. | in1(18)/in1(15)<-3.)
  nrlobad = nrlobad + 1;
 end
 if (in1(8)/in1(11)<-2. | in1(15)/in1(18)<-2.| in1(11)/in1(8)<-2. | in1(18)/in1(15)<-2.)
  nrlo2bad = nrlo2bad + 1;
 end
 end
 if (cpbad(nog)>40)
   nog40 = nog40 + 1;
 end
 if (cpbad(nog)>50)
   nog50 = nog50 + 1;
 end
 if (cpbad(nog)>60)
   nog60 = nog60 + 1;
 end
end 	% if ifskip==0
end	% if nogood==0 for all

end	% do nt=1:num
nmiss,ntru,n,nzpyb,nrlo,nrlo2,nog1,nogs,nogn,nog,nzybad0,nrlobad0,nrlo2bad0,nog30,nzybad,nrlobad,nrlo2bad,nog40,nog50,nog60

% Recalc medians for whole period
 medby2 = median(by2)
 medcpn = median(cpn)
 medcps = median(cps)
 medoff24n = median(offcn)
 medoff24s = median(offcs)
 aveoff24n = mean(offcn)
 aveoff24s = mean(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)
 medhpn = median(hpn)
 medhps = median(hps)
 if (nog>0)
  medcpbad = median(cpbad)
 end
 medcp10 = median(cp10)
 medcrad10 = median(crad10)
 medcp15 = median(cp15)
 medcrad15 = median(crad15)
 medcp20 = median(cp20)
 medcrad20 = median(crad20)
 medcp25 = median(cp25)
 medcrad25 = median(crad25)
 medcp30 = median(cp30)
 medcrad30 = median(crad30)
 medcp40 = median(cp40)
 medcrad40 = median(crad40)
 medcp60 = median(cp60)
 medcrad60 = median(crad60)
 medcp80 = median(cp80)
 medcrad80 = median(crad80)
 medcp90 = median(cp90)
 medcrad90 = median(crad90)
 medhp10 = median(hp10)
 medarad10 = median(arad10)
 medamcrad10 = median(amcrad10)
 medhp15 = median(hp15)
 medarad15 = median(arad15)
 medamcrad15 = median(amcrad15)
 medhp20 = median(hp20)
 medarad20 = median(arad20)
 medamcrad20 = median(amcrad20)
 medhp25 = median(hp25)
 medarad25 = median(arad25)
 medamcrad25 = median(amcrad25)
 medhp30 = median(hp30)
 medarad30 = median(arad30)
 medamcrad30 = median(amcrad30)
 medhp40 = median(hp40)
 medarad40 = median(arad40)
 medamcrad40 = median(amcrad40)
 medhp60 = median(hp60)
 medarad60 = median(arad60)
 medamcrad60 = median(amcrad60)
 medhp80 = median(hp80)
 medarad80 = median(arad80)
 medamcrad80 = median(amcrad80)
 medhp90 = median(hp90)
 medarad90 = median(arad90)
 medamcrad90 = median(amcrad90)
 medoffc10 = median(offc10)
 medoffc15 = median(offc15)
 medoffc20 = median(offc20)
 medoffc25 = median(offc25)
 medoffc30 = median(offc30)
 medoffc40 = median(offc40)
 medoffc60 = median(offc60)
 medoffc80 = median(offc80)
 medoffc90 = median(offc90)
 medhpnps = median(0.5*(hpn+hps))
 medamcradnps = median(0.5*(amcradn+amcrads))
 medoff24nps = median(0.5*(offcn+offcs))
 medoff18nps = median(0.5*(dskofcn+dskofcs))
 avoff24nps = mean(0.5*(offcn+offcs))
 avoff18nps = mean(0.5*(dskofcn+dskofcs))
 



fsize = 12;

% 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);
ravcpr=polyfit(cprlo,cradrlo,1)
xcprlo = polyval(ravcpr,cpt);
medradcprlo = median(cradrlo)
ravcpr2=polyfit(cprlo2,cradrlo2,1)
xcprlo2 = polyval(ravcpr2,cpt);
medradcprlo2 = median(cradrlo2)
ravcpzy=polyfit(cpzpyb,cradzpyb,1)
xcpzpyb = polyval(ravcpzy,cpt);
medradcpzy = median(cradzpyb)
rbcpt=polyfit(cpb,cradb,1)
xbcpt = polyval(rbcpt,cpt);


% Plot histogram of cpbad and good (includes 1bad)
% n_elem =   CP<10kV    10-15    15-20    20-25    25-30    30-40     40-60    60-100   100-400kV
%n_elem cpb= 688+559=1247 1760+676=2436 2515+573=3088 3159+442=3601 2647+313=2960 5112+253=5365 7704+93=7797 8973+36=9009 4642+0=4642(>100kV) = 40145 total (correct!)
% median CP and crad for good
%     8.6/12.635  12.9/11.765  17.4/12.060  22.7/11.820  27.4/11.9  34.65/11.845  49.35/12.20
% 2247 cpbad<25kV (560*12.64+676*11.76+573*12.06+438*11.82)/2247 = 27115.7/2247 = 12.06 deg
% 2911 cpbad<60kV (27115.7+314*11.9+255*11.84+95*12.20)/2911 = (27115.7+7914.8)/2911 = 35030.5/2911 = 12.03 deg - so choose 12.0 deg for default!
edges = [0.,10.,15.,20.,25.,30.,40.,60.,100.,400.];

n_bad = histc(cpbad,edges)
n_good = histc(cpb,edges)
n_1bad = histc(cp1bad,edges)
% Revise n_good to delete the n_1bad points
for n=1:10
 n_good(n) = n_good(n) - n_1bad(n);
end
%n_good

% Set size to be 8.5x11 inches
set(gcf,'position',[142    99   647   822])
H = gcf;
 subplot('Position',[0.12,0.66,0.92,0.28]) 
 hold on;
%title('Good calccloc CMIT-MIX for 06348-349, 07226-253, 08080-106','FontSize',15,'FontWeight','bold');
title('Calccloc in CMIT-MIX 06348-349, 07226-253, 08080-106','FontSize',15,'FontWeight','bold');
bar(edges,n_good);
 ylabel ('# both good','FontSize',fsize,'Color','k');
 grid;

 subplot('Position',[0.12,0.38,0.92,0.28]) 
 hold on;
%title('One Bad calccloc CMIT-MIX for 06348-349, 07226-253, 08080-106','FontSize',15,'FontWeight','bold');
bar(edges,n_1bad);
 ylabel ('# one bad (made good)','FontSize',fsize,'Color','k');
 grid;

 subplot('Position',[0.12,0.10,0.92,0.28]) 
 hold on;
%title('Both Bad calccloc CMIT-MIX for 06348-349, 07226-253, 08080-106','FontSize',15,'FontWeight','bold');
bar(edges,n_bad);
 ylabel ('# both bad','FontSize',fsize,'Color','k');
 xlabel ('Lower bound for CP (kV) bin','FontSize',fsize,'Color','k');
 grid;

hold off;

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('CMIT-MIX for 06348-349, 07226-253, 08080-106','FontSize',15,'FontWeight','bold');
%plot (cps,crads,'xb');
%plot (cpn,cradn,'+r');
 plot (cpb,cradb,'.k');
 plot (cpt,xscpt,'-b');
 plot (cpt,xncpt,'-r');
 plot (cpt,xbcpt,'-k');
 plot (cpt,xcprlo,'-m');
 plot (cpt,xcprlo2,'-g');
 plot (cpt,xcpzpyb,'-c');
 plot (medcp10,medcrad10,'dr');
 plot (medcp15,medcrad15,'dr');
 plot (medcp20,medcrad20,'dr');
 plot (medcp25,medcrad25,'dr');
 plot (medcp30,medcrad30,'dr');
 plot (medcp40,medcrad40,'dr');
 plot (medcp60,medcrad60,'dr');
 plot (medcp80,medcrad80,'dr');
 plot (medcp90,medcrad90,'dr');
%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');
 plot (0.5*cpx,11.0,'dr');
 text (0.51*cpx,11.0,'medians in CP bins for histograms','Color','r','FontSize',fsize,'FontWeight','bold');
text (0.60*cpx,11.9,'median N,S=47.9,48.0 kV','Color','k','FontSize',fsize,'FontWeight','bold');
 text (0.60*cpx,10.1,'median N,S=12.75,12.87 deg','Color','k','FontSize',fsize,'FontWeight','bold');
 text (0.3*cpx,9.0,'Use 12.0deg default radius from histogram of bad pts','Color','k','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('CMIT-MIX for 06348-349, 07226-253, 08080-106','FontSize',15,'FontWeight','bold');
%plot (cplo,cradlo,'xk');
%plot (cpt,xcplo,'-k');
 plot (cpzpyb,cradzpyb,'xc');
 plot (cpt,xcpzpyb,'-c');
 plot (cprlo2,cradrlo2,'xg');
 plot (cpt,xcprlo2,'-g');
 plot (cprlo,cradrlo,'xm');
 plot (cpt,xcprlo,'-m');
 ylabel ('Convection Radius (deg)','FontSize',fsize,'Color','k');
 xlabel ('Big By CP (kV)','FontSize',fsize,'Color','k');
% Use numbers of nog/ntru=2945/40145,nzybad0=280/nog,nrlo2bad0=2028/nog,nrlobad0=1234/nog
% Bad>30 is nog30=381 (ratio nog30/nog=13.0%, use nzy(nrlo(2))bad/nog
text (0.09*cpx,20.5,'7.3% bad (12.9%>30kV, 9.5% Bz+ |By|>1.2Bz, 69%,42% cpn/x,x/n<-2,-3)','Color','k','FontSize',fsize,'FontWeight','bold');
% Do not know what t means (1bad?)
%text (0.46*cpx,10.5,'26%g,25%t (29% bad>30kV) Bz+ |By|>1.2Bz','Color','c','FontSize',fsize,'FontWeight','bold');
%text (0.46*cpx,9.5,'11%g,12%t (51% bad>30kV) cpn/x<-2','Color','g','FontSize',fsize,'FontWeight','bold');
%text (0.46*cpx,8.5,'6%g,6%t (37% bad>30kV) cpn/x<-3','Color','m','FontSize',fsize,'FontWeight','bold');
% Use ratios of nzpyb/n(good)=9647/37200,nrlo2=9249/n, nrlo=3523/n
text (0.46*cpx,10.5,'25.5% good Bz+ |By|>1.2Bz','Color','c','FontSize',fsize,'FontWeight','bold');
text (0.46*cpx,9.5,'24.9% good cpn/x,x/n<-2','Color','g','FontSize',fsize,'FontWeight','bold');
text (0.46*cpx,8.5,'9.5% good cpn/x,x/n<-3','Color','m','FontSize',fsize,'FontWeight','bold');
grid;
set (gca,'xlim',[0 cpx],'ylim',[8 radx],'FontSize',fsize,'box','on');

hold off;

% Plot arad vs CP and arad-crad vs CP and HP vs CP
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('CMIT-MIX for 06348-349, 07226-253, 08080-106','FontSize',15,'FontWeight','bold');
 plot (cpb,aradb,'.k');
 plot (medcp10,medarad10,'dr');
 plot (medcp15,medarad15,'dr');
 plot (medcp20,medarad20,'dr');
 plot (medcp25,medarad25,'dr');
 plot (medcp30,medarad30,'dr');
 plot (medcp40,medarad40,'dr');
 plot (medcp60,medarad60,'dr');
 plot (medcp80,medarad80,'dr');
 plot (medcp90,medarad90,'dr');
 plot (0.5*cpx,11.2,'dr');
 text (0.51*cpx,11.2,'medians in CP bins for histograms','Color','r','FontSize',fsize,'FontWeight','bold');
 text (0.60*cpx,10.0,'median N,S=16.74,16.98 deg','Color','k','FontSize',fsize,'FontWeight','bold');
 ylabel ('Auroral 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('CMIT-MIX for 06348-349, 07226-253, 08080-106','FontSize',15,'FontWeight','bold');
 plot (cpb,amcradb,'.k');
 plot (medcp10,medamcrad10,'dr');
 plot (medcp15,medamcrad15,'dr');
 plot (medcp20,medamcrad20,'dr');
 plot (medcp25,medamcrad25,'dr');
 plot (medcp30,medamcrad30,'dr');
 plot (medcp40,medamcrad40,'dr');
 plot (medcp60,medamcrad60,'dr');
 plot (medcp80,medamcrad80,'dr');
 plot (medcp90,medamcrad90,'dr');
 plot (0.5*cpx,0.5,'dr');
 text (0.51*cpx,0.5,'medians in CP bins for histograms','Color','r','FontSize',fsize,'FontWeight','bold');
 text (0.60*cpx,-0.5,'median N,S=3.75,4.20 deg','Color','k','FontSize',fsize,'FontWeight','bold');
 ylabel ('Auroral-Convection Radius (deg)','FontSize',fsize,'Color','k');
 xlabel ('CP (kV)','FontSize',fsize,'Color','k');
grid;
set (gca,'xlim',[0 cpx],'ylim',[-1 10],'FontSize',fsize,'box','on');

hold off;

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('CMIT-MIX for 06348-349, 07226-253, 08080-106','FontSize',15,'FontWeight','bold');
 plot (cpb,hpb,'.k');
 plot (medcp10,medhp10,'dr');
 plot (medcp15,medhp15,'dr');
 plot (medcp20,medhp20,'dr');
 plot (medcp25,medhp25,'dr');
 plot (medcp30,medhp30,'dr');
 plot (medcp40,medhp40,'dr');
 plot (medcp60,medhp60,'dr');
 plot (medcp80,medhp80,'dr');
 plot (medcp90,medhp90,'dr');
 plot (0.5*cpx,57.,'dr');
 text (0.51*cpx,57.,'medians in CP bins for histograms','Color','r','FontSize',fsize,'FontWeight','bold');
 text (0.50*cpx,52.,'median N,S 1.12,1.36 GW','Color','k','FontSize',fsize,'FontWeight','bold');
 ylabel ('Hemispheric Power (GW)','FontSize',fsize,'Color','k');
 xlabel ('Cross-tail Potential Drop (kV)','FontSize',fsize,'Color','k');
grid;
set (gca,'xlim',[0 cpx],'ylim',[0 60],'FontSize',fsize,'box','on');
 subplot('Position',[0.07,0.54,0.85,0.40]) 
 hold on;
title('CMIT-MIX for 06348-349, 07226-253, 08080-106','FontSize',15,'FontWeight','bold');
 plot (cpb,offcb,'.k');
 plot (medcp10,medoffc10,'dr');
 plot (medcp15,medoffc15,'dr');
 plot (medcp20,medoffc20,'dr');
 plot (medcp25,medoffc25,'dr');
 plot (medcp30,medoffc30,'dr');
 plot (medcp40,medoffc40,'dr');
 plot (medcp60,medoffc60,'dr');
 plot (medcp80,medoffc80,'dr');
 plot (medcp90,medoffc90,'dr');
 plot (0.5*cpx,5.0,'dr');
 text (0.51*cpx,5.0,'medians in CP bins for histograms','Color','r','FontSize',fsize,'FontWeight','bold');
 text (0.50*cpx,3.5,'median,av -0.345,-0.430 deg','Color','k','FontSize',fsize,'FontWeight','bold');
 ylabel ('Offset towards 24MLT (deg)','FontSize',fsize,'Color','k');
 xlabel ('Cross-tail Potential Drop (kV)','FontSize',fsize,'Color','k');
grid;
set (gca,'xlim',[0 cpx],'ylim',[-7 7],'FontSize',fsize,'box','on');
hold off;
