clear all;

% nmph = number of minutes per hour (30x2min WHI, 12x5min AGU Dec06)
% dmin = number of minutes between CMIT files 

% Read in the first hdf file to obtain grid and other attributes and set date for plot label
%    sd_id = hdfsd( 'start', [filenames(1,:)], 'rdonly');
%  sd_id = hdfsd( 'start', '/hao/arc2/schmitt/data/LTR-2_1_4/Agu2006/CMITs/AGU_CMITs_mix_2006-12-15T00-00-00Z.hdf', 'rdonly');

  iyear = 2006
% iyear = 2008
% iyear = 2007
% ifrd1 = 1 when read 'outcmit_single_time.txt'
  ifrd1 = 0

if (iyear==2006)
  dmin = 5
  daynum0 = 348
  daynum = 348
%  daynum = 349
  month = 12
  daymon = daynum - 334
  uthr = 22
  utmn = 20
% 06 348 Dec14, 2100-2225, (94kV max 2110UT w By~-9.6, Bz~14.2)
%  sd_id = hdfsd( 'start', '/hao/aim4/schmitt/data/LTR-2_1_4/Agu2006/CMITs/AGU_CMITs_mix_2006-12-14T22-20-00Z.hdf', 'rdonly');
  mix_name = sprintf('/hao/aim4/schmitt/data/LTR-2_1_4/Agu2006/CMITs/AGU_CMITs_mix_%2.2d-%2.2d-%dT%2.2d-%2.2d-00Z.hdf',fix(iyear),fix(month),fix(daymon),fix(uthr),fix(utmn));
   sd_id = hdfsd( 'start', mix_name, 'rdonly');
end

if (iyear==2007)
  dmin = 2
  daynum0 = 226
%  sd_id = hdfsd( 'start', '/hao/aim3/emery/cmit/CR2060_CMIT/CR2060-CMIT_mix_2007-08-26T21-56-00Z.hdf', 'rdonly');
%  daynum = 238
%  uthr = 21
%  utmn = 56
% 07239 8/27 0024,0052,0054,0100(max 98kV)-0110 w By~-10, Bz~+4
%  sd_id = hdfsd( 'start', '/hao/aim3/emery/cmit/CR2060_CMIT/CR2060-CMIT_mix_2007-08-27T01-00-00Z.hdf', 'rdonly');
  daynum = 239
  uthr = 1
  utmn = 0
% 07238 8/26 2018 79S,100NkV(crad=26.43deg) but both min<12! w By~-12, Bz~+7
  sd_id = hdfsd( 'start', '/hao/aim3/emery/cmit/CR2060_CMIT/CR2060-CMIT_mix_2007-08-26T20-18-00Z.hdf', 'rdonly');
  daynum = 238
  uthr = 20
  utmn = 18
end

if (iyear==2008)
% 08095 Apr 4, 4 1216-1310UT(max 72kV) for By+8,Bz+5 
  sd_id = hdfsd( 'start', '/hao/aim2/wiltbemj/mhd_data/WHI-CMIT-final/WHI-CMIT-final_mix_2008-04-04T12-56-00Z.hdf', 'rdonly');
  dmin = 2
  daynum0 = 80
  daynum = 95
  uthr = 12
  utmn = 56
end

nmph = 60/dmin
nskp = (daynum-daynum0)*nmph*24 + uthr*nmph + utmn/dmin

% Open IMF output
if (iyear==2006)
% filimf = sprintf('imf_out_dec06')	% Feb 2012 up to day 349 0455UT
% filloc = sprintf('outcmit_agu2006b')  % Feb 2012 up to day 349 0455UT
 filimf = sprintf('imf_out_dec06b')	% Aug 2012 all of 348-349
%filloc = sprintf('outcmit_agu2006b_2355ut')   % Aug 2012 all of 348-349 SH noon=midnit
 filloc = sprintf('outcmit_dec2006_testc')    % Sep 2012

end
if (iyear==2007)
 filimf = sprintf('imf_out_07226_253')
%filloc = sprintf('outcmit_cr2060_2007b')
 filloc = sprintf('outcmit_2007cr2060_all_testc')    % Sep 2012
end
if (iyear==2008)
 filimf = sprintf('imf_out_08080_106')
%filloc = sprintf('outcmit_whi2008b')
 filloc = sprintf('outcmit_whi08_all_2testc')    % Sep 2012
end
if (ifrd1==1)
 filloc = sprintf('outcmit_single_time.txt')
end

imfin = fopen(filimf,'r');

% Open driver output
locin = fopen(filloc,'r');
% Skip first 2 lines
linloc = fgetl(locin)
linloc = fgetl(locin)
% Skip over nskp times (nskp*2 lines)
 for nl=1:nskp
 linimf = fgetl(imfin);
  if (ifrd1==0)
   linloc1 = fgetl(locin);
   linloc2 = fgetl(locin);
  end
 end
% Read the IMF file
 linimf = fgetl(imfin)
 inimf = sscanf(linimf,'%f');
 bxint = inimf(5);
 byint = inimf(6);
 bzint = inimf(7);
 vswint = inimf(8);
 dswint = inimf(9);
%  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')") 
% 19-24 added for offcen, dskofcen, radcen
%  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.
% Read file from driver.f90 and calccloc.f90
 linloc1 = fgetl(locin);
 in1 = sscanf(linloc1,'%f');
 linloc2 = fgetl(locin);
 in2 = sscanf(linloc2,'%f');
%  NOTE SH MLT wrong in driver.f90 since noon=midnight, but 6 and 18 MLT OK
%    Need to correct MLT and offc for Aug 2012 outputs (NOT for Sep 2012!)
% in2(19)=-in2(19);
% in2(7) = -in2(7);
% All values reset to defaults for various ifbad>0 (in2 values 25 and 26 for SH and NH)
% Set crad (theta0) convection radius and center
 xcs = -in2(11);	% -dskofc_deg_SH
 ycs = -in2(7);		% -offc_deg_SH
 Rcs = in2(3);		% crad_deg_SH
 xcn = -in2(12);	% -dskofc_deg_NH
 ycn = -in2(8);		% -offc_deg_NH
 Rcn = in2(4);		% crad_deg_NH
% Convert back to polar (LAT in deg, MLT in radians)
  [MLTcents,LATcents]=cart2pol(xcs,ycs);
  [MLTcentn,LATcentn]=cart2pol(xcn,ycn);
%LATcents,LATcentn,MLTcents,MLTcentn
% Look at default (offc=0.1, dskofc=0.0) positions
  xc0 = 0.;
  yc0 = -0.4;
% Calculate Rc10, the av(crads+cradn) (but not less than 10deg) for crit1(+5),2(+20) 
 Rc10 = max(10.0,0.5*(in2(3)+in2(4)));
% Revised to 12 deg if ifbad>0 both hems
 if (in2(25)>0 & in2(26)>0)
   Rc10 = 12;
 end
% Convert back to polar (LAT in deg, MLT in radians)
  [MLT0,LAT0]=cart2pol(xc0,yc0);
% Look at values from min-max line (cen values 19-24)
 xcens = -in2(21);	% -dskofc_deg_SH
 ycens = -in2(19);		% -offc_deg_SH
 Rcens = in2(23);		% crad_deg_SH
 xcenn = -in2(22);	% -dskofc_deg_NH
 ycenn = -in2(20);		% -offc_deg_NH
 Rcenn = in2(24);		% crad_deg_NH
% Convert back to polar
  [MLTcens,LATcens]=cart2pol(xcens,ycens);
  [MLTcenn,LATcenn]=cart2pol(xcenn,ycenn);

% Set max/min potential points SH 7-13, NH 14-20
%  NOTE SH MLT wrong in driver.f90 since noon=midnight, but 6 and 18 MLT OK
%    Need to correct MLT and offc
%   if (in1(10) <= 12.)
%     d06 = 6.-in1(10);
%     in1(10) = 6.+d06;
%   else
%     d18 = 18.-in1(10);
%     in1(10) = 18.+d18;
%   end
%   if (in1(13) <= 12.)
%     d06 = 6.-in1(13);
%     in1(13) = 6.+d06;
%   else
%     d18 = 18.-in1(13);
%     in1(13) = 18.+d18;
%   end 
   rhoxs = in1(10)*pi/12. - pi/2.;
   rhons = in1(13)*pi/12. - pi/2.;
   rhoxn = in1(17)*pi/12. - pi/2.;
   rhonn = in1(20)*pi/12. - pi/2.;
   colatxs = 90.-abs(in1(9));
   colatns = 90.-abs(in1(12));
   colatxn = 90.-in1(16);
   colatnn = 90.-in1(19);

    alpha = hdfsd('readattr', sd_id, 5);
    beta  = hdfsd('readattr', sd_id, 6);
    R     = hdfsd('readattr', sd_id, 7);
   
    % Read in the X(rows27,cols181) grid (0 to ~19.9Re? for 181 columns)
    sds_id = hdfsd( 'select', sd_id,hdfsd('nametoindex',sd_id,'Grid X'));
    X = hdfsd('readdata',sds_id,zeros(1,2),[],[181; 27]);
    
    % Read in the Y(27,181) grid
    sds_id = hdfsd( 'select', sd_id,hdfsd('nametoindex',sd_id,'Grid Y'));
    Y = hdfsd('readdata',sds_id,zeros(1,2),[],[181 27]);

    %x_gridded(1,i,:,:,1) = X;
    %x_gridded(2,i,:,:,1) = X;
    %y_gridded(1,i,:,:,1) = Y;
    %y_gridded(2,i,:,:,1) = Y;
   
    % Calculate Radius(27,181) 181 columns all the same, rows 1-27 from 0-.6947(0 to 39.8deg co-lat?)
    radius = sqrt(X.^2 + Y.^2);
% make colatdeg 39.8 to 0 deg (reverse order) and colat in radians
for j=1:27
 colat(j) = radius(j,2);
 colatdeg(j) = radius(28-j,2)*180/pi;
 mlat2pol(j) = 90. - colatdeg(j);
 mlat2polx(j+1) = 90.-colatdeg(j);
end
% The dmlat is irregular from 50.1990,51.6616,..89.75,89.90,90.0 mlat
 mlat2polx(1) = mlat2polx(2)-(mlat2polx(3)-mlat2polx(2));
mlat2polx
    
    % Calculate co-latitude
% 181 theta is not co-latitude, but is 0 to +pi, then -pi to 0 so is mlon or ~MLT
% Want to plot with midnight at bottom, or -pi/2 to 3pi/2 (-90 to 270 deg)
%  theta=0 is plotted at 6MLT rhs usually so add +pi/2 for NH and -pi/2 for SH in plot.py
%  epotn first theta is at noon, while epots first theta is at midnight
%  09/12:  This is wrong.  theta=0 is noon NH, theta=90 is 18 MLT NH
%   while theta=0 is noon SH, theta=90 is 6 MLT SH
    theta = atan2(Y,X);
    % Note: the following syntax might not work... May need to explicitly loop over theta for periodic point:
%    theta[theta < 0] = theta[theta < 0] + 2*pi
%length(theta),length(radius)
for i=1:length(theta)
  if theta(i)<0
    theta(i) = theta(i) + 2*pi;
  end
end
%theta*180/pi
htr = pi/12.;
dtr = pi/180.;
nimlt = 181;
for i=1:nimlt
  mltn(i) = 12.+theta(2,i)/htr;
   if (mltn(i)>=24.) 
     mltn(i) = mltn(i)-24.;
   end
% 09/12 corrected mlts time going backwards from noon
  mlts(i) = 12.-theta(2,i)/htr;
   if (mlts(i)<0.) 
     mlts(i) = mlts(i)+24.;
   end
   if (mlts(i)>=24.) 
     mlts(i) = mlts(i)-24.;
   end
% Want mlt=0-24
   mlt(i) = theta(2,i)/htr;
% Change mlt so is always positive (goes 0-12, then -12 to 0)
   if mlt(i)<0
     mlt(i) = mlt(i) + 24.;
   end
   rho(i) = mlt(i)*pi/12. - pi/2.;
end


% Get xp(27,181),yp(27,181) coordinates for contourf (from polarcontf)
for j=1:27
 for i=1:181
  xp(j,i) = colat(j)*cos(theta(i));
  yp(j,i) = colat(j)*sin(theta(i));
 end
end

%    hdfsd('end',sd_id);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Now allocate & read all the data stored in a list of MIX output files:     
%    f_gridded = zeros(2,length(filenames),27,181);
    f_gridded = zeros(2,1,27,181);
    x_gridded = f_gridded; 
    y_gridded = f_gridded;
   
%    for i=1:length(filenames)
i = 1;
    
%    sd_id = hdfsd( 'start', [filenames(i,:)], 'rdonly');
%   sd_id = hdfsd( 'start', '/hao/arc2/schmitt/data/LTR-2_1_4/Agu2006/CMITs/AGU_CMITs_mix_2006-12-15T00-00-00Z.hdf', 'rdonly');
% stime in seconds?
    stime(i) = datenum(double(hdfsd('readattr', sd_id, 2)));
%stime1 = stime(i)

    % Read in electric potential: V (volt)
    sds_id = hdfsd( 'select', sd_id,hdfsd('nametoindex',sd_id,'Potential North [V]'));
% This makes epotn(181mlt,27nlat)???? but size epotn = (27,181)!!!! or (nlat,nmlt)
% first time is noon for epotn
    epotn = hdfsd('readdata',sds_id,zeros(1,2),[],[181 27]);
    sds_id = hdfsd( 'select', sd_id,hdfsd('nametoindex',sd_id,'Potential South [V]'));
% firt time is midnight for epots
    epots = hdfsd('readdata',sds_id,zeros(1,2),[],[181 27]);
%length(epots)		% This gives 181, the first dimension
%  Change epotn,s from 0-40 colat to 40-0 colat;  add 180 deg to epotn so starts 0 MLT;
%    change epotn so goes from 0-24 not noon-0-noon backwards in time;convert from V to kV
% Switch SH to later times is i(noon_later)=182-i(noon_earlier)
%  Add 12 hours is +90 or -90

for i=1:181
for j=1:27
 if (i<91)
  epotNH(i,j) = double(epotn(28-j,i+90)*0.001);
  epotSH(i,j) = double(epots(28-j,92-i)*0.001);
 else
  epotNH(i,j) = double(epotn(28-j,i-90)*0.001);
  epotSH(i,j) = double(epots(28-j,272-i)*0.001);
 end
end
end	%  for i
% Get vector max, then max and same for min
 maxvsh = max(epotSH);
 maxs = max(maxvsh);
 minvsh = min(epotSH);
 mins = min(minvsh);
 maxvnh = max(epotNH);
 maxn = max(maxvnh);
 minvnh = min(epotNH);
 minn = min(minvnh);
 cpn = maxn-minn
 cps = maxs-mins
   
    % Read in energy: keV (kilo - electron volt)
    sds_id = hdfsd( 'select', sd_id,hdfsd('nametoindex',sd_id,'Average energy North [keV]'));
    ENERGY = hdfsd('readdata',sds_id,zeros(1,2),[],[181 27]);
 
    % Read in flux
    sds_id = hdfsd( 'select', sd_id,hdfsd('nametoindex',sd_id,'Number flux North [1/cm^2 s]'));
    FLUX = hdfsd('readdata',sds_id,zeros(1,2),[],[181 27]);
    
%    f_gridded(1,i,:,:) = ENERGY;
%    f_gridded(2,i,:,:) = FLUX;

% Read SH ENERGY and #FLUX
    sds_id = hdfsd( 'select', sd_id,hdfsd('nametoindex',sd_id,'Average energy South [keV]'));
    ENERGYS = hdfsd('readdata',sds_id,zeros(1,2),[],[181 27]);
    sds_id = hdfsd( 'select', sd_id,hdfsd('nametoindex',sd_id,'Number flux South [1/cm^2 s]'));
    FLUXS = hdfsd('readdata',sds_id,zeros(1,2),[],[181 27]);
%  Change from 0-40 colat to 40-0 colat, then add 180 deg so starts 0
%    and calculate energy flux NH and SH in mW/m2
%  eflux=flux*fac_p2e*alfa   (fac_p2e=2*1.602e-9, so 2*alfa in this factor)
fac_p2e = 1.602e-9 % convert from particle to energy flux with MEAN (not alfa) ENERGY in keV
for i=1:181
for j=1:27
 if (i<91)
  eflxNH(i,j) = double(FLUX(28-j,i+90)*fac_p2e*ENERGY(28-j,i+90));
  eflxSH(i,j) = double(FLUXS(28-j,92-i)*fac_p2e*ENERGYS(28-j,92-i));
 else
  eflxNH(i,j) = double(FLUX(28-j,i-90)*fac_p2e*ENERGY(28-j,i-90));
  eflxSH(i,j) = double(FLUXS(28-j,272-i)*fac_p2e*ENERGYS(28-j,272-i));
 end
end
end	%  for i

% Calculate MIX HP in GW
%  Re = radius of the Earth in m at the level of the ionosphere (use 6370e+3 in MIX)
%  Could use 6816 to increase HP by 7% if use 
      Re = 6370.e+3;
      areatot = 0.;
%  The MIX grid is 2 deg in maglon (irreg in maglat) 
      dmlon = 360./180.;

% 06348 2220 UT:  hemisphericPower.py 22.77S, 17.15N
%  mlat2polx is j=1-28, where j=28 is pole and j=2 is bottom lat in MIX grid, so j is below, and j+1 is at the lat for other j values.
% 06348 2220 UT:  Gives 20.21SH, 15.09NH; areatot = 6.1588e+13, da1=2.3145e+10, da26=4.3145e+6 using area above j
%      for j=1:26
%       da(j) = Re*cos(mlat2pol(j)*dtr)*dmlon*dtr*(Re*dtr*(mlat2pol(j+1)-mlat2pol(j)));
% 06348 2220 UT:  Gives 20.00SH, 14.97NH; areatot = 6.0728e+13, da1=2.3145e+10, da26=6.4719e+6
%  areatot should have been more??
       for j=1:27
%      da(j) = Re*cos(mlat2polx(j+1)*dtr)*dmlon*dtr*(Re*dtr*(mlat2polx(j+1)-mlat2polx(j)));
% 06348 2220 UT:  Gives 21.16SH, 16.02NH j=26 or 27; areatot = 6.3311e+13, da1=2.3499e+10, da26=11.326e+6 using area half below
       da(j) = Re*cos(0.5*(mlat2polx(j+1)+mlat2polx(j))*dtr)*dmlon*dtr*(Re*dtr*(mlat2polx(j+1)-mlat2polx(j)));
% 06348 2220 UT:  Gives 21.16SH, 16.02NH j=26 or 27; areatot = 6.3311e+13, da1=2.3499e+10, da26=11.326e+6 using area below
%       da(j) = Re*cos(mlat2polx(j)*dtr)*dmlon*dtr*(Re*dtr*(mlat2polx(j+1)-mlat2polx(j)));
       areatot = areatot + da(j)*180.;
j,da(j),areatot
      end

areatot

      ghpi_sh_mix = 0.;
      ghpi_nh_mix = 0.;
      for i=1:180
      for j=1:27
%       for j=1:26
        ghpi_sh_mix = ghpi_sh_mix + da(j)*eflxSH(i,j);
        ghpi_nh_mix = ghpi_nh_mix + da(j)*eflxNH(i,j);
       end
      end

% Convert from mW/m2 * m2(da) * 1.e-12GW/mW to GW
      ghpi_sh_mix = ghpi_sh_mix*1.e-12  
      ghpi_nh_mix = ghpi_nh_mix*1.e-12  

% Now calculate arad from eflx 
%  Calculate the locations of the maximum eflx and get arad from mlatdeg 12h apart 
%  For some large By cases, get a large eflux at the pole on the dayside, so must
%   find the secondary max on the nightside (2-6 MLT)
       for i=1:180
         maxefxSH = -1.;
         maxefxNH = -1.;
         for j=1:27
           if (eflxSH(i,j) > maxefxSH) 
	     maxefxSH = eflxSH(i,j);
             mlatefxS(i) = mlat2pol(j);
           end
           if (eflxNH(i,j) > maxefxNH) 
	     maxefxNH = eflxNH(i,j);
             mlatefxN(i) = mlat2pol(j);
           end
         end
       end
%  For some large By cases, get a large eflux at the pole on the dayside, so must
%   find the secondary max on the nightside (for 2-6 MLT calc) equatorwards of 75mlat
%    mlat2pol(15)=74.2071
    
%  Calculate for 2-6 MLT only
	 aradSH = 0.;
         narad = 0;
	 aradNH = 0.;
         for i=1:180
	   if (mlt(i) >= 2. & mlt(i) <= 6.) 
             j = i+90;
	     if (j > 180) 
               j = j - 180;
             end
	     narad = narad + 1;
% This next bit of code is unneccessary, since if abs(aradSH-NH)>4deg, set both to max
% Find secondary max equatorwards of 75mlat if 14-18 MLT and 2-6 MLT mlatefx >80mlat
%        if (mlatefxS(i) > 80. | mlatefxS(j) > 80.) 
%	  maxefxSH = -1.
%          for j=1:15
%           if (eflxSH(i,j) > maxefxSH) 
%            maxefxSH = eflxSH(i,j);
%             mlatefxS(i) = mlat2pol(j);
%           end
%          end
%        end
%        if (mlatefxN(i) > 80. | mlatefxN(j) > 80.) 
%	  maxefxNH = -1.;
%          for j=1:15
%           if (eflxNH(i,j) > maxefxNH) 
%	     maxefxNH = eflxNH(i,j);
%             mlatefxN(i) = mlat2pol(j);
%           end
%          end
%        end
	     aradSH = aradSH + 0.5*(180.-mlatefxS(i)-mlatefxS(j));
	     aradNH = aradNH + 0.5*(180.-mlatefxN(i)-mlatefxN(j));
           end
         end
narad
       if (narad > 0) 
        aradSH = max(0.1,aradSH / narad)
        aradNH = max(0.1,aradNH / narad)
% Revise if difference >4deg so smaller arad same as bigger arad
        raSH = aradSH;
	raNH = aradNH;
        if (aradNH-aradSH > 4.) 
           raSH = aradNH
        end
        if (aradSH-aradNH > 4.)
	   raNH = aradSH
        end	
% Calculate Ra10
      Ra10 = max(10.,0.5*(raSH+raNH)-4.);  
       end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 12/2011 Barbara Emery HAO/NCAR: plot MIX convection NH and SH
%  Need epotn,s (MLT 0-24, colat 0-44deg) for plotting

hf=figure
 set(hf,'position',[162    41   700   906]); 

subplot('position',[.1,.5,.45,.45]);
      XM=1.35;
% ~emery/tieconvec/tie2_conv_aur.m h=polarcontf(trho(nlat),ttheta(nmlt),(epot(nmlt,nlat)'),min:del:max) and this WORKS FINE!
% nlat=81, nmlt=48, trho=40 to 0 (colat), nmlt=-90 to +270 or 0 to 24MLT in rad
% THIS FAILS because: ??? Undefined function or method 'contourc' for input arguments of type 'single'. WHY???
       h = polarcontf(colatdeg,rho,(epotNH'),-90:10:70);
%      h = polarcontf(colat,theta,epotn,-90:10:70);
%      h = polarcontf(colat,theta,epotn,17,-90:10:70);
%     [C,h] = contourf(epotn);
%      [C,h,cf] = contourf(xp,yp,epotn);
      caxis([-90 70]);
      axis equal; 
      axis off;
% Print out title etc
    hold on;
 title_label = sprintf('%2.2d%2.2d UT %d %d NH MIX %dkV=%d-%d %4.1fBy %4.1fBz %dV %5.2fGW',uthr,utmn,daynum,iyear,fix(cpn),fix(maxn),fix(minn),byint,bzint,fix(vswint),ghpi_nh_mix);
     tit=3*pi/4;
     x=65*cos(tit);
     y=65*sin(tit);
     text(x,y,title_label,'FontName','Times','FontSize',16,'color','k');
         lab=pi/4;
         x=65*cos(lab);
         y=-65*sin(lab);
         text(x,y,' kV','FontName','Times','FontSize',14,'color','k');
      set(gca,'layer','bottom');
      v={'-90';'-80';'-70';'-60';'-50';'-40';'-30';'-20';'-10';'0  ';'10 ';'20 ';'30 ';'40 ';'50 ';'60 ';'70 '};
      hc=colorbar;      
      pos=get(hc,'position'); 
      pos(1)=1.19*pos(1); pos(3)=0.9*pos(3); pos(4)=0.95*pos(4);
      set(hc,'position',pos,'ylim',[-90 70],'ytick',([-90:10:70]),'yticklabel',v);
      %%%%%%%%%%%%%%%%%%%  to set background (MLAT circle and MLT text)
      PolarBackgroundx;
% Plot max/min potential points SH 7-13, NH 14-20
      polar (rhoxn,colatxn,'mx');
      polar (rhonn,colatnn,'ms');
% Plot convection radius (crad or theta0), crit(1) and crit(2) for offcenter and pole
      axis equal;
      rectangle('position',[xcn-Rcn,ycn-Rcn,Rcn*2,Rcn*2],...
          'curvature',[1,1],'linestyle','--','linewidth',2.5,'edgecolor','y');
      polar(MLTcentn,LATcentn,'y.');
      rectangle('position',[xc0-Rcn,yc0-Rcn,Rcn*2,Rcn*2],...
          'curvature',[1,1],'linestyle','--','linewidth',2.5,'edgecolor','k');
      polar(MLT0,LAT0,'k.');
% Plot crit1,2 from convection (default 12 deg)
      rectangle('position',[xc0-Rc10-5.,yc0-Rc10-5.,(Rc10+5.)*2,(Rc10+5.)*2],...
          'curvature',[1,1],'linestyle','--','edgecolor','k');
      rectangle('position',[xc0-Rc10-20.,yc0-Rc10-20.,(Rc10+20.)*2,(Rc10+20.)*2],...
          'curvature',[1,1],'linestyle','--','edgecolor','k');
% Plot crit1,2 from aurora (default Rc=Ra-4, where Ra revised to max if diff by > 4 deg)
      rectangle('position',[xc0-Ra10-5.,yc0-Ra10-5.,(Ra10+5.)*2,(Ra10+5.)*2],...
          'curvature',[1,1],'linestyle','--','edgecolor','m');
      rectangle('position',[xc0-Ra10-20.,yc0-Ra10-20.,(Ra10+20.)*2,(Ra10+20.)*2],...
          'curvature',[1,1],'linestyle','--','edgecolor','m');
      axis equal;


subplot('position',[.1,.05,.45,.45]);
       h = polarcontf(colatdeg,rho,(epotSH'),-90:10:70);
%      h = polarcontf(colat,theta,epots,17,-90:10:70); 
      caxis([-90 70]);
      axis equal; 
      axis off;
% Print out title etc
    hold on;
 title_label = sprintf('%2.2d%2.2d UT %d %d SH MIX %dkV=%d-%d %4.1fBy %4.1fBz %dV %5.2fGW',uthr,utmn,daynum,iyear,fix(cps),fix(maxs),fix(mins),byint,bzint,fix(vswint),ghpi_sh_mix);
     tit=3*pi/4;
     x=65*cos(tit);
     y=65*sin(tit);
     text(x,y,title_label,'FontName','Times','FontSize',16,'color','k');
         lab=pi/4;
         x=65*cos(lab);
         y=-65*sin(lab);
         text(x,y,' kV','FontName','Times','FontSize',14,'color','k');
      set(gca,'layer','bottom');
      v={'-90';'-80';'-70';'-60';'-50';'-40';'-30';'-20';'-10';'0  ';'10 ';'20 ';'30 ';'40 ';'50 ';'60 ';'70 '};
      hc=colorbar;      
      pos=get(hc,'position'); 
      pos(1)=1.19*pos(1); pos(3)=0.9*pos(3); pos(4)=0.95*pos(4);
      set(hc,'position',pos,'ylim',[-90 70],'ytick',([-90:10:70]),'yticklabel',v);
      %%%%%%%%%%%%%%%%%%%  to set background (MLAT circle and MLT text)
      PolarBackgroundx;
% Plot max/min potential points SH 7-13, NH 14-20
      polar (rhoxs,colatxs,'mx');
      polar (rhons,colatns,'ms');
% Plot convection radius (crad or theta0), crit(1) and crit(2) for offcenter and pole
      axis equal;
      rectangle('position',[xcs-Rcs,ycs-Rcs,Rcs*2,Rcs*2],...
          'curvature',[1,1],'linestyle','--','linewidth',2.5,'edgecolor','y');
      polar(MLTcents,LATcents,'y.');
      rectangle('position',[xc0-Rcs,yc0-Rcs,Rcs*2,Rcs*2],...
          'curvature',[1,1],'linestyle','--','linewidth',2.5,'edgecolor','k');
      polar(MLT0,LAT0,'k.');
% Plot crit1,2 from convection (default 12 deg)
      rectangle('position',[xc0-Rc10-5.,yc0-Rc10-5.,(Rc10+5.)*2,(Rc10+5.)*2],...
          'curvature',[1,1],'linestyle','--','edgecolor','k');
      rectangle('position',[xc0-Rc10-20.,yc0-Rc10-20.,(Rc10+20.)*2,(Rc10+20.)*2],...
          'curvature',[1,1],'linestyle','--','edgecolor','k');
% Plot crit1,2 from aurora (default Rc=Ra-4, where Ra revised to max if diff by > 4 deg)
      rectangle('position',[xc0-Ra10-5.,yc0-Ra10-5.,(Ra10+5.)*2,(Ra10+5.)*2],...
          'curvature',[1,1],'linestyle','--','edgecolor','m');
      rectangle('position',[xc0-Ra10-20.,yc0-Ra10-20.,(Ra10+20.)*2,(Ra10+20.)*2],...
          'curvature',[1,1],'linestyle','--','edgecolor','m');
      axis equal;

hold off;

% Plot energy flux
power = 0.5*(ghpi_nh_mix+ghpi_sh_mix)
% Can only do 11 colors in v to go from blue to red
pwrmax = 15.;
pwrmin = 0.2;
pwrmid = 5.;
if (power<pwrmin)
 maxgw = 0.1;
 dgw = 0.01;
end
if (power>=pwrmin & power<pwrmid)
 maxgw = 0.5;
 dgw = 0.05;
end
if (power>=pwrmid & power<pwrmax)
 maxgw = 2.5;
 dgw = 0.25;
end
if (power>=pwrmax)
 maxgw = 5.;
 dgw = 0.5;
end

hf=figure
 set(hf,'position',[162    41   700   906]); 

subplot('position',[.1,.5,.45,.45]);
      XM=1.35;

       h = polarcontf(colatdeg,rho,(eflxNH'),0:dgw:maxgw);
      caxis([0 maxgw]);
      axis equal; 
      axis off;
% Print out title etc
    hold on;
 title_label = sprintf('%2.2d%2.2d UT %d %d NH MIX %dkV=%d-%d %4.1fBy %4.1fBz %dV %5.2fGW',uthr,utmn,daynum,iyear,fix(cpn),fix(maxn),fix(minn),byint,bzint,fix(vswint),ghpi_nh_mix);
     tit=3*pi/4;
     x=65*cos(tit);
     y=65*sin(tit);
     text(x,y,title_label,'FontName','Times','FontSize',16,'color','k');
         lab=pi/4;
         x=65*cos(lab);
         y=-65*sin(lab);
         text(x,y,'mW/m^2','FontName','Times','FontSize',14,'color','k');
      set(gca,'layer','bottom');
if (power<pwrmin)
     v={'0.00';'0.01';'0.02';'0.03';'0.04';'0.05';'0.06';'0.07';'0.08';'0.09';'0.10'};
end
if (power>=pwrmin & power<pwrmid)
     v={'0';'0.05';'0.10';'0.15';'0.20';'0.25';'0.30';'0.35';'0.40';'0.45';'0.50'};
end
if (power>=pwrmid & power<pwrmax)
      v={'0';'0.25';'0.5';'0.75';'1';'1.25';'1.5';'1.75';'2';'2.25';'2.5'};
end
if (power>=pwrmax)
    v={'0';'0.5';'1';'1.5';'2';'2.5';'3';'3.5';'4';'4.5';'5'};
end
      hc=colorbar;      
      pos=get(hc,'position'); 
      pos(1)=1.19*pos(1); pos(3)=0.9*pos(3); pos(4)=0.95*pos(4);
      set(hc,'position',pos,'ylim',[0 maxgw],'ytick',([0:dgw:maxgw]),'yticklabel',v);
      %%%%%%%%%%%%%%%%%%%  to set background (MLAT circle and MLT text)
      PolarBackgroundw;
% Plot convection radius (crad or theta0), crit(1) and crit(2) for offcenter and pole
      axis equal;
      rectangle('position',[xcn-Rcn,ycn-Rcn,Rcn*2,Rcn*2],...
          'curvature',[1,1],'linestyle','--','linewidth',2.5,'edgecolor','y');
      polar(MLTcentn,LATcentn,'y.');
%      rectangle('position',[xcn-aradNH,ycn-aradNH,aradNH*2,aradNH*2],...
%          'curvature',[1,1],'linestyle','-','linewidth',2.5,'edgecolor','y');
      rectangle('position',[xc0-Rcn,yc0-Rcn,Rcn*2,Rcn*2],...
          'curvature',[1,1],'linestyle','--','linewidth',2.5,'edgecolor','w');
      polar(MLT0,LAT0,'w.');
      rectangle('position',[xc0-aradNH,yc0-aradNH,aradNH*2,aradNH*2],...
          'curvature',[1,1],'linestyle','-','linewidth',2.5,'edgecolor','m');
% Plot crit1,2 from convection (default 12 deg)
      rectangle('position',[xc0-Rc10-5.,yc0-Rc10-5.,(Rc10+5.)*2,(Rc10+5.)*2],...
          'curvature',[1,1],'linestyle','--','edgecolor','w');
      rectangle('position',[xc0-Rc10-20.,yc0-Rc10-20.,(Rc10+20.)*2,(Rc10+20.)*2],...
          'curvature',[1,1],'linestyle','--','edgecolor','w');
% Plot crit1,2 from aurora (default Rc=Ra-4, where Ra revised to max if diff by > 4 deg)
      rectangle('position',[xc0-Ra10-5.,yc0-Ra10-5.,(Ra10+5.)*2,(Ra10+5.)*2],...
          'curvature',[1,1],'linestyle','--','edgecolor','m');
      rectangle('position',[xc0-Ra10-20.,yc0-Ra10-20.,(Ra10+20.)*2,(Ra10+20.)*2],...
          'curvature',[1,1],'linestyle','--','edgecolor','m');
      axis equal;

subplot('position',[.1,.05,.45,.45]);

       h = polarcontf(colatdeg,rho,(eflxSH'),0:dgw:maxgw);
      caxis([0 maxgw]);
      axis equal; 
      axis off;
% Print out title etc
    hold on;
 title_label = sprintf('%2.2d%2.2d UT %d %d SH MIX %dkV=%d-%d %4.1fBy %4.1fBz %dV %5.2fGW',uthr,utmn,daynum,iyear,fix(cps),fix(maxs),fix(mins),byint,bzint,fix(vswint),ghpi_sh_mix);
     tit=3*pi/4;
     x=65*cos(tit);
     y=65*sin(tit);
     text(x,y,title_label,'FontName','Times','FontSize',16,'color','k');
         lab=pi/4;
         x=65*cos(lab);
         y=-65*sin(lab);
         text(x,y,'mW/m^2','FontName','Times','FontSize',14,'color','k');
      set(gca,'layer','bottom');
if (power<pwrmin)
     v={'0.00';'0.01';'0.02';'0.03';'0.04';'0.05';'0.06';'0.07';'0.08';'0.09';'0.10'};
end
if (power>=pwrmin & power<pwrmid)
     v={'0';'0.05';'0.10';'0.15';'0.20';'0.25';'0.30';'0.35';'0.40';'0.45';'0.50'};
end
if (power>=pwrmid & power<pwrmax)
      v={'0';'0.25';'0.5';'0.75';'1';'1.25';'1.5';'1.75';'2';'2.25';'2.5'};
end
if (power>=pwrmax)
    v={'0';'0.5';'1';'1.5';'2';'2.5';'3';'3.5';'4';'4.5';'5'};
end
     hc=colorbar;         
      pos=get(hc,'position'); 
      pos(1)=1.19*pos(1); pos(3)=0.9*pos(3); pos(4)=0.95*pos(4);
      set(hc,'position',pos,'ylim',[0 maxgw],'ytick',([0:dgw:maxgw]),'yticklabel',v);
      %%%%%%%%%%%%%%%%%%%  to set background (MLAT circle and MLT text)
      PolarBackgroundw;
% Plot convection radius (crad or theta0), crit(1) and crit(2) for offcenter and pole
      axis equal;
      rectangle('position',[xcs-Rcs,ycs-Rcs,Rcs*2,Rcs*2],...
          'curvature',[1,1],'linestyle','--','linewidth',2.5,'edgecolor','y');
      polar(MLTcents,LATcents,'y.');
%      rectangle('position',[xcs-aradSH,ycs-aradSH,aradSH*2,aradSH*2],...
%          'curvature',[1,1],'linestyle','-','linewidth',2.5,'edgecolor','y');
      rectangle('position',[xc0-Rcs,yc0-Rcs,Rcs*2,Rcs*2],...
          'curvature',[1,1],'linestyle','--','linewidth',2.5,'edgecolor','w');
      polar(MLT0,LAT0,'w.');
      rectangle('position',[xc0-aradSH,yc0-aradSH,aradSH*2,aradSH*2],...
          'curvature',[1,1],'linestyle','-','linewidth',2.5,'edgecolor','m');
% Plot crit1,2 from convection (default 12 deg)
      rectangle('position',[xc0-Rc10-5.,yc0-Rc10-5.,(Rc10+5.)*2,(Rc10+5.)*2],...
          'curvature',[1,1],'linestyle','--','edgecolor','w');
      rectangle('position',[xc0-Rc10-20.,yc0-Rc10-20.,(Rc10+20.)*2,(Rc10+20.)*2],...
          'curvature',[1,1],'linestyle','--','edgecolor','w');
% Plot crit1,2 from aurora (default Rc=Ra-4, where Ra revised to max if diff by > 4 deg)
      rectangle('position',[xc0-Ra10-5.,yc0-Ra10-5.,(Ra10+5.)*2,(Ra10+5.)*2],...
          'curvature',[1,1],'linestyle','--','edgecolor','m');
      rectangle('position',[xc0-Ra10-20.,yc0-Ra10-20.,(Ra10+20.)*2,(Ra10+20.)*2],...
          'curvature',[1,1],'linestyle','--','edgecolor','m');
      axis equal;

hold off;

      axis off;
    %Close the file buffer when you are finished.
%    hdfsd('end',sd_id);
