%%
%% This program is written for plotting the TIEGCM dres TEC data at 8 lonsra
%% BA Emery, October 2011
clear all; close all

iftec = 1;
ifhmf2 = 0;

% Data
 loncutsE=[25,90,140,175,200,250,285,345]; % specify climatology longitudes to plot
if iftec==0
dirname='/hao/aim3/emery/tiegcm_dres/data/cosmic_dec06/'; % directory with COSMIC NmF2,hmF2 data
 loncuts=[25,90,140,175,200,250,285,345]; % specify climatology longitudes to plot
 n15 = 4;
 n1 = 3;
 n2 = 95;
else
%dirname='/hao/aim3/emery/tiegcm_dres/data/mit_climatology/'; % directory with TEC data
dirname='/hao/aim3/emery/tiegcm_dres/data/mit_dec06/wtedonly/'; % directory with TEC data
 loncuts=[25,90,140,175,160,110,75,15]; % specify longitudes to plot 4 E and 4 W
 n15 = 1;
 n1 = 1;
 n2 = 96;
end
latgps=-90:5:90;
min15=0:0.25:23.75;
min60=0.5:1:23.5;
length(min15),length(latgps)
% Dec 13-15, 2006 is daynumbers 347-349
startday=datenum(2006,12,14,0,0,0); % specify first day to plot
    [year1,month1,day1,hrs1,min1,sec1]=datevec(startday);

% Model
% longitude slices at begin lons:  25,90,140,175,-160,-110,-75,-15 for 2.5 deg dres
jj8 = [84 110 130 144 10 30 44 68];
latdres=-90:2.5:87.5;
% year1 = 2006;
 nday1 = 348;
 nday2 = 349;

% Loop over number of days
for ida=nday1:nday2-1
 idastr = sprintf('06%d',ida)
if iftec==1
 ratmd = NaN(length(min15),length(latgps),8);
 pcmd = NaN(length(min15),length(latgps),8);
else
 ratmd = NaN(length(min60),length(latgps),8);
 pcmd = NaN(length(min60),length(latgps),8);
end
% Data
        tecarrmean=NaN(length(min15),length(latgps),8);
        tecarrmed=NaN(length(min15),length(latgps),8);
	hmf2arrmed=NaN(length(min60),length(latgps),8);
	nmf2arrmed=NaN(length(min60),length(latgps),8);
	day1 = day1 + ida-nday1;
figure;
% position [left bottom width height]
% Set size to be 8.5x11 inches
set(gcf,'position',[142    99   647   822])
H = gcf;
    clf
    for j=1:8
        clear year month day hrs mins lat lonmin lonmax tecmean tecmed dtec stderr npts hmf2med nmf2med stdhmf2 stdnmf2
        loncut=loncuts(j);
      if iftec==1
% TEC
        if j<5
        filename=[dirname,'long_',num2str(loncut),...
            'E_',num2str(year1),'_',num2str(month1),'_',num2str(day1),'.dat'];
        else
        filename=[dirname,'long_',num2str(loncut),...
            'W_',num2str(year1),'_',num2str(month1),'_',num2str(day1),'.dat'];
	end
            [year,month,day,hrs,mins,lat,lonmin,lonmax,tecmean,dtec,stderr,npts]=textread(filename,'%d %d %d %f %f %f %f %f %f %f %f %d');
%            [year,month,day,hrs,mins,lat,lonmin,lonmax,tecmean,tecmed,dtec,stderr,npts]=textread(filename,'%d %d %d %f %f %f %f %f %f %f %f %f %d');
            %
            dates=datenum(year,month,day,hrs,mins,0);
            datesun=unique(dates);
            datestart=datenum(year(1),month(1),day(1),hrs(1),mins(1),0);
            datesun1=datesun-datestart;
            %
            for i=1:length(tecmean)
% Emery:  Toss if dtec>4 or stderr>5 TECU - not much difference for 90E 22-Nov-2007 or 16-Dec-2007
	      if dtec(i)<4 & stderr(i)<5
                idates=find(datesun == dates(i));
                ilat=find(latgps == lat(i));
                tecarrmean(idates,ilat,j)=tecmean(i);
%                tecarrmed(idates,ilat,j)=tecmed(i);
              end
            end		% of i=1:length(tecmean)
	else
% COSMIC
        filenamh=[dirname,'hmF2_long_',num2str(loncut),'_',idastr,'.dat'];
        filenamn=[dirname,'NmF2_long_',num2str(loncut),'_',idastr,'.dat'];
            [year,daynum,hrs,mins,secs, lat,lonmin,hmf2med,stdhmf2, npts]=textread(filenamh,'%d %d %d %d %d %f %f %f %f %d');
            [year,daynum,hrs,mins,secs, lat,lonmin,nmf2med,stdnmf2, npts]=textread(filenamn,'%d %d %d %d %d %f %f %f %f %d');
	    i = 0;
	    for n=1:length(min60)
	     for k=1:length(latgps)-1
	      i = i + 1;
	      if (hmf2med(i)>50.)
	       hmf2arrmed(n,k,j) = hmf2med(i);
	      end
	      if (nmf2med(i)>1.)
% Convert from cm-3 to 1.e-12 m-3
	       nmf2arrmed(n,k,j) = nmf2med(i)*1.e-6;
	      end
	     end
	    end
	end
%% Make plots
            subplot(4,2,j)
       if iftec==1
            pcolor(datesun1,latgps,tecarrmean(:,:,j)')
%            caxis([0 30]);		% for climatology
            caxis([0 50]);		% for Dec06 storm
            shading flat
            %xlabel('Date');
            ylabel('Latitude');
        if j<5
            title(['GPS TEC wted mean, ',datestr(datestart,'dd-mmm-yyyy'),', LON=',num2str(loncut),'E']);
	else
            title(['GPS TEC wted mean, ',datestr(datestart,'dd-mmm-yyyy'),', LON=',num2str(loncut),'W']);
	end
            set(gca,'TickDir','out');
            datetick('x','HH:MM');
       else
	 if ifhmf2==1
            pcolor(min60,latgps,hmf2arrmed(:,:,j)')
            caxis([0 500]);		% for Dec06 storm
            title(['COSMIC hmF2 (km) ',num2str(year1),' ',num2str(ida),' LON=',num2str(loncuts(j))]);
	 else
            pcolor(min60,latgps,nmf2arrmed(:,:,j)')
            caxis([0 2.5]);		% for Dec06 storm
            title(['COSMIC NmF2*1.e-12 (m^3) ',num2str(year1),' ',num2str(ida),' LON=',num2str(loncuts(j))]);
         end
            shading flat
            %xlabel('Date');
            ylabel('Latitude');
	    set (gca,'xlim',[0,24],'Xtick',([0:6:24]),'XTickLabel',{'00:00','06:00','12:00','18:00','24:00'});
	end
    end				% end of plotting 8 longitudes
    hc=colorbar('horiz');
    hc_pos=get(hc,'Position');
    hc_pos(1)=0.35;
    hc_pos(2)=0.065;
    set(hc,'Position', hc_pos);
if iftec==1
    plotname=['GPSTECwmean_',idastr];
else
  if ifhmf2==1
    plotname=['COSMIC_hmf2med_',idastr];
  else
    plotname=['COSMIC_nmf2med_',idastr];
  end
end
% This worked first time only, next try bombed matlab and linux
%   print('-dpng',plotname);
hold off;

ida
%stop

%%  Read TIEGCM .nc files starting at 0000 UT every 15 min and make lon slice plots
Fname= sprintf('/hao/ftp/pub/emery/dres/dres_s_06%d.nc',ida)

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

     jj = jj8(j);
%  nmax in 1.e-12 m-3
    [nt,dayfrac,nmax,hmax,tec,nz,nlat,ht,htne,lon] = get_teclon_nc(Fname,jj);
%nt,nlat/2+1
%  06347 0 UT is bad tec, nmax, hmax so delete
if (ida==347 | ida==350)
%hmax(1,1)
  for i=1:nlat
    tec(1,i) = NaN;
    nmax(1,i) = NaN;
    hmax(1,i) = NaN;
  end
end
% have length 72 for some reason??? so do not print!
%maxnmf2 = max(nmax(:,:));
%minnmf2 = min(nmax(:,:));
%maxhmf2 = max(hmax(:,:));
%minhmf2 = min(hmax(:,:));

% Calculate ratio of model/data and percent: lat=1:2:nlat; time=0.5*(t+t+1)
     rmodat = 0;
     percentmd = 0;
     ndat = 0;
     for i=1:2:nlat
      i2 = (i+1)/2;
      m = 0;
      for n=n1:n15:n2
       m = m + 1;
       if iftec==1
	dat = tecarrmean(m,i2,j);
        modn = tec(n,i);
	if (n<nt)
	 modnp = tec(n+1,i);
	end
       else
	if ifhmf2==1
         dat = hmf2arrmed(m,i2,j);
	 modn = hmax(n,i);
	 modnp = hmax(n+1,i);
	else 
         dat = nmf2arrmed(m,i2,j);
	 modn = nmax(n,i);
	 modnp = nmax(n+1,i);
        end
       end       
        if (dat>=0.001)
	 ndat = ndat + 1;
         np1 = 0;
         if (ida==347 & n==1)
	   np1 = 1;
	 end
         if (ida==350 & n==1)
	   np1 = 1;
	 end
         if (np1==1)
	   tie = modnp;
	 else
	   if (n==nt)
	    tie = modn;
           else
	    tie = 0.5*(modn+modnp);
	   end
	 end
	ratmd(m,i2,j) = tie/dat;
        pcmd(m,i2,j) =100.*(tie/dat-1.);
        if (tie/dat < 1.)
          pcmd(m,i2,j) = -100.*(dat/tie-1.);
        end
	 rmodat = rmodat + tie/dat;
         percentmd = percentmd + pcmd(m,i2,j);
%if (tie/dat > 2. | tie/dat < 0.5)
%tie/dat,pcmd(m,i2,j)
%end
	end
      end
     end
	rmodat = rmodat/ndat
	percentmd = percentmd/ndat
rmd8(j) = rmodat;
pcmd8(j) = percentmd;
% 53-67 COSMIC npts vs ~1150-2650 TEC GPS npts
ndat
%% Make plots
            subplot(4,2,j)
	if iftec==1
            pcolor(dayfrac,latdres,tec')
%            caxis([0 30]);		% for climatology
            caxis([0 50]);		% for Dec06 storm
            title(['TIEGCM dres TEC ',num2str(year1),' ',num2str(ida),' LON=',num2str(lon(jj)),' M/D=',num2str(rmodat)]);
	else
	  if ifhmf2==1
            pcolor(dayfrac,latdres,hmax')
            caxis([0 500]);		% for Dec06 storm
            title(['TIEGCM dres hmF2 ',num2str(year1),' ',num2str(ida),' LON=',num2str(lon(jj)),' M/D=',num2str(rmodat)]);
	  else
            pcolor(dayfrac,latdres,nmax')
            caxis([0 2.5]);		% for Dec06 storm
            title(['TIEGCM dres NmF2 ',num2str(year1),' ',num2str(ida),' LON=',num2str(lon(jj)),' M/D=',num2str(rmodat)]);
	  end
	end
            shading flat
            %xlabel('Date');
            ylabel('Latitude');
	    set (gca,'xlim',[ida,ida+1],'Xtick',([ida:0.25:ida+1]),'XTickLabel',{'00:00','06:00','12:00','18:00','24:00'});
    end				% end of plotting 8 longitudes
    hc=colorbar('horiz');
    hc_pos=get(hc,'Position');
    hc_pos(1)=0.35;
    hc_pos(2)=0.065;
    set(hc,'Position', hc_pos);
if iftec==1
    plotname=['TIEdres_',idastr,'_TEC'];
else
  if ifhmf2==1
    plotname=['TIEdres_',idastr,'_hmf2'];
  else
    plotname=['TIEdres_',idastr,'_nmf2'];
  end
end
% This works, but matlab dies immediately so skip
%    print('-dpng',plotname);
hold off;

% Plot ratio of model/data
figure;
% position [left bottom width height]
% Set size to be 8.5x11 inches
set(gcf,'position',[142    99   647   822])
H = gcf;
clf
for j=1:8
%% Make plots
            subplot(4,2,j)
	if iftec==1
            pcolor(min15,latgps,ratmd(:,:,j)')		
            title(['model/data TEC ',num2str(year1),' ',num2str(ida),' LON=',num2str(loncutsE(j)),' M/D=',num2str(rmd8(j))]);
	else
            pcolor(min60,latgps,ratmd(:,:,j)')		
	  if ifhmf2==1
            title(['model/data hmF2 ',num2str(year1),' ',num2str(ida),' LON=',num2str(loncuts(j)),' M/D=',num2str(rmd8(j))]);
	  else	
            title(['model/data NmF2 ',num2str(year1),' ',num2str(ida),' LON=',num2str(loncuts(j)),' M/D=',num2str(rmd8(j))]);
	  end
	end
            caxis([0.5 2]);
            shading flat
            %xlabel('Date');
            ylabel('Latitude');
	    set (gca,'xlim',[0,24],'Xtick',([0:6:24]),'XTickLabel',{'00:00','06:00','12:00','18:00','24:00'});
    end				% end of plotting 8 longitudes
    hc=colorbar('horiz');
    hc_pos=get(hc,'Position');
    hc_pos(1)=0.35;
    hc_pos(2)=0.065;
    set(hc,'Position', hc_pos);
hold off;

% Plot percent of model/data
figure;
% position [left bottom width height]
% Set size to be 8.5x11 inches
set(gcf,'position',[142    99   647   822])
H = gcf;
clf
for j=1:8
%% Make plots
            subplot(4,2,j)
	if iftec==1
            pcolor(min15,latgps,pcmd(:,:,j)')		
            title(['%model/data TEC ',num2str(year1),' ',num2str(ida),' LON=',num2str(loncutsE(j)),' %M/D=',num2str(pcmd8(j))]);
	else
            pcolor(min60,latgps,pcmd(:,:,j)')		
	  if ifhmf2==1
            title(['%model/data hmF2 ',num2str(year1),' ',num2str(ida),' LON=',num2str(loncuts(j)),' %M/D=',num2str(pcmd8(j))]);
	  else	
            title(['%model/data NmF2 ',num2str(year1),' ',num2str(ida),' LON=',num2str(loncuts(j)),' %M/D=',num2str(pcmd8(j))]);
	  end
	end
            caxis([-100 100]);
            shading flat
            %xlabel('Date');
            ylabel('Latitude');
	    set (gca,'xlim',[0,24],'Xtick',([0:6:24]),'XTickLabel',{'00:00','06:00','12:00','18:00','24:00'});
    end				% end of plotting 8 longitudes
    hc=colorbar('horiz');
    hc_pos=get(hc,'Position');
    hc_pos(1)=0.35;
    hc_pos(2)=0.065;
    set(hc,'Position', hc_pos);

end	% for ida=nday1:nday2-1


