clear all;

   % Read in ascii epots,n file made in driver.f90 and plot with calccloc values
% nlat=27, nmlt=181 (181=1), trho=40 to 0 (colat), nmlt=-90 to +270 or 0 to 24MLT in rad

% open 2 files 
% filloc = sprintf('outcmit_06348_0300utb')
%  filepot = sprintf('outcmit_06348_3ut_epots_n.txt')
 filloc = sprintf('outcmit_single_time.txt')
  filepot = sprintf('single_time_epot.txt')
  locin = fopen(filloc,'r');
  potin = fopen(filepot,'r');
% Skip 3 extra lines in new potin
 linpot = fgetl(potin)
 linpot = fgetl(potin)
 linpot = fgetl(potin)

% Skip first 2 lines
linloc = fgetl(locin)
linloc = fgetl(locin)
%  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.

n = 1;
% Read file from driver.f90 and calccloc.f90
 linloc1 = fgetl(locin);
 in1 = sscanf(linloc1,'%f');
 linloc2 = fgetl(locin);
 in2 = sscanf(linloc2,'%f');
 fracdaymo(n) = in1(4)+in1(5)/24.+in1(6)/(24.*60.);
 fracdaymo(n) = in1(2)+in1(5)/24.+in1(6)/(24.*60.);
 year(n) = in1(1);
 daynum(n) = in1(2);
 uthr(n) = in1(5);
 utmn(n) = in1(6);
 cps(n) = in1(7);
 cpsx(n) = in1(8);
 cpsn(n) = in1(11);
 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);
 cpnx(n) = in1(15);
 cpnn(n) = in1(18);
 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);
% These are cartesian coordinates:
  xfitCs(n) = -dskofcs(n);
  yfitCs(n) = -offcs(n);
  RfitCs(n) = crads(n);
  xfitCn(n) = -dskofcn(n);
  yfitCn(n) = -offcn(n);
  RfitCn(n) = cradn(n);
% Convert back to polar
  [MLTcentCs(n),LATcentCs(n)]=cart2pol(xfitCs(n),yfitCs(n));
  [MLTcentCn(n),LATcentCn(n)]=cart2pol(xfitCn(n),yfitCn(n));


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

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

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

% Get number for 5 min OMNI data to skip over
nd0 = (daynum(n)-1)*12*24 + uthr(n)*12 + utmn(n)/5;
byimf = byyr(nd0+1)
bzimf = bzyr(nd0+1)
vsw = vswyr(nd0+1)
ydfrac = dateyr(nd0+1)
iyr= fix(ydfrac/1000.)
nda = fix(ydfrac-iyr*1000)
frac = ydfrac-iyr*1000-nda
uthr = frac*24.

% Read epots,n

%  write (12,"(a,' colatdeg, mlthr epotSH(181,27),NH')") datehrmn 
%  write (12,"(10f7.2)") colatdeg
%  write (12,"(10f7.2)") mlts
% do j=1,njlat
%  write (12,"(10f7.1)") (epotSH(i,j),i=1,nimlt)
% enddo
% do j=1,njlat
%  write (12,"(10f7.1)") (epotNH(i,j),i=1,nimlt)
% enddo

% Skip date and descr line
 linpot = fgetl(potin)
% Read 27 colatdeg 40-0
 linpot1 = fgetl(potin);
 inp1 = sscanf(linpot1,'%f');
 for n=1:10
  colatdeg(n) = inp1(n);
 end
 linpot1 = fgetl(potin);
 inp1 = sscanf(linpot1,'%f');
 for n=11:20
  colatdeg(n) = inp1(n-10);
 end
 linpot1 = fgetl(potin);
 inp1 = sscanf(linpot1,'%f');
 for n=21:27
  colatdeg(n) = inp1(n-20);
 end
% Read 181 mlt 0-24 and convert to radians -pi/2 to +3pi/2
 for n10=1:18
 linpot1 = fgetl(potin);
 inp1 = sscanf(linpot1,'%f');
 n2 = n10*10;
 n1 = n2 - 9;
  for n=n1:n2
   mlt(n) = inp1(n1-n2+10);
   theta(n) = mlt(n)*pi/12. - pi/2.;
  end
 end
 linpot1 = fgetl(potin);
 inp1 = sscanf(linpot1,'%f');
  mlt(181) = inp1(1);
  theta(181) = 3*pi/2.;
% Read epots,n
 for j=1:27
 for n10=1:18
 linpot2 = fgetl(potin);
 inp2 = sscanf(linpot2,'%f');
 n2 = n10*10;
 n1 = n2 - 9;
  for n=n1:n2
   epotSH(n,j) = inp2(n1-n2+10);
  end
 end
 linpot2 = fgetl(potin);
 inp2 = sscanf(linpot2,'%f');
   epotSH(181,j) = inp2(1);
 end
 for j=1:27
 for n10=1:18
 linpot2 = fgetl(potin);
 inp2 = sscanf(linpot2,'%f');
 n2 = n10*10;
 n1 = n2 - 9;
  for n=n1:n2
   epotNH(n,j) = inp2(n1-n2+10);
  end
 end
 linpot2 = fgetl(potin);
 inp2 = sscanf(linpot2,'%f');
   epotNH(181,j) = inp2(1);
 end

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

k = 1;
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,theta,(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 2006 NH MIX %dkV=%d-%d Bz=%5.1f By=%5.1f V=%d',uthr(k),utmn(k),daynum(k),fix(cpn(k)),fix(cpnx(k)),fix(cpnn(k)),bzimf,byimf,fix(vsw));
     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 convection circle red 
      rectangle('position',[xfitCn(k)-RfitCn(k),yfitCn(k)-RfitCn(k),RfitCn(k)*2,RfitCn(k)*2],...
           'curvature',[1,1],'linestyle','-','edgecolor','r');
      polar(MLTcentCn(k),LATcentCn(k),'r.');
      axis equal;    

subplot('position',[.1,.05,.45,.45]);
       h = polarcontf(colatdeg,theta,(epotSH'),-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 2006 SH MIX %dkV=%d-%d Bz=%5.1f By=%5.1f V=%d',uthr(k),utmn(k),daynum(k),fix(cps(k)),fix(cpsx(k)),fix(cpsn(k)),bzimf,byimf,fix(vsw));
     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 convection circle red 
      rectangle('position',[xfitCs(k)-RfitCs(k),yfitCs(k)-RfitCs(k),RfitCs(k)*2,RfitCs(k)*2],...
           'curvature',[1,1],'linestyle','-','edgecolor','r');
      polar(MLTcentCs(k),LATcentCs(k),'r.');
      axis equal;    
hold off;
      axis off;
