%%
%% This program is written for plotting the mean IGS TEC data at 8 lons
%% BA Emery, May 2012 from /aim3/emery/dec06/igs/plot_8lon.m
clear all; close all


% if5deg = 1,0 for 5 deg or 2.5 deg dlat (only possible for IGS and TIE-WT)
if5deg = 0;
if (if5deg==1)
 nl36=36;
else
 nl36=71;
end

ndays = 3;

% ifmit=1,0 for MIT/JPL or IGS TEC
ifmit = 1;
if (ifmit==0)
 nt = 12;
 ifjpl = 0;
else
 nt = 96;
% ifjpl = 1,0 for JPL or other TEC (MIT or IGS)
 ifjpl = 1;
% ifjpl = 0;
end

 vartit = ' TEC '


% Data
 loncutsE=[25,90,140,175,200,250,285,345]; % specify climatology longitudes to plot
dlat = 5;		% 36*5deg=180 deg from 90S-90N		
lat5deg=-90:dlat:85;
lat2p5deg = -87.5:2.5:87.5;	

%min120 = 1:2:23;
% Make 13 values, and then duplicate last as first in arrays or otherwise plot only 11 values - not 12!
x2 = 24*ndays;
min120 = 0:2:24*ndays;
min60=0:0.25:24*ndays;
if (ifmit==0)
 mingps=min120;
 if (if5deg==0)
  dlat = 2.5;
  latgps=lat2p5deg;
 else
  latgps = lat5deg;
 end
else
% MIT and JPL data
 mingps=min60;	
 latgps = lat5deg;
end

dlon = 5;	% dlon = 360/24=15deg for hourly; dlon=15/4=2.75 for 15 min
nlon = 360/dlon;
nlat = length(latgps)

% calculate global density (and tec) and compare with Emmert at 3 altitudes
re=6371.e5;   %  earth radius in cm
dtr = pi/180.;
areatot = 0;
for ii=1:nlat
% Use midpt latitude (+dlat/2.)
 da(ii)=re*cos((latgps(ii)+dlat/2.)*dtr)*(dlon*dtr)*(re*dlat*dtr);
 areatot = areatot+da(ii)*nlon;
end

% imod = (0)1-5 for (NRL IGS) CODE, ESOC, JPL, UPC, and IGS (where ..rg files are older than ...g files)
%imod = 0
%for imod=0:4

% Data
        numtec=NaN(length(mingps),length(latgps),8);
        tecarrmed=NaN(length(mingps),length(latgps),8);

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

if ifmit==1
 if ifjpl==0
   gpsloc=' MIT';
   DTECnam = sprintf('/hao/aim3/emery/dec06/igs/TEC_lon_%3.3d_MIT_2006.347-2006.349.dat',loncutsE(j))
 else
   gpsloc=' JPL';
   DTECnam = sprintf('/hao/aim3/emery/dec06/igs/TEC_lon_%3.3d_JPL_2006.347-2006.349.dat',loncutsE(j))
 end

% year daynum beginning_uthr min sec beg_glat beg_glon tecavwt tecmed averr tecstd num
% -999.0 for missing with num=0
DTECin = fopen(DTECnam,'r');
%  Skip first line
 linegps = fgetl(DTECin);
datarea(j) = 0;
gdat(j) = 0;
for nh=1:nt*ndays
 for nl=1:nlat
 linegps = fgetl(DTECin);
 ingps = sscanf(linegps,'%f');
  numtec(nh,nl,j) = ingps(12);
  if (numtec(nh,nl,j) > 0)
% 8=tecavwt and 9=tecmed
   if (ingps(9) > 0.) 
    tecarrmed(nh,nl,j) = ingps(9);
   end
  else
    numtec(nh,nl,j) = 0;
  end
% Calculate 'global' mean TEC from 24h of 1 glon
% Weight with area of globe and find global gdat values
	if numtec(nh,nl,j)>0 & tecarrmed(nh,nl,j)>0.
         datarea(j) = datarea(j) + da(nl);
         gdat(j) = gdat(j) + da(nl)*tecarrmed(nh,nl,j);
        end
 end	% for 1:nlat
end	% for nh=1:nt
  gdat(j) = gdat(j)/datarea(j);
else
% Read IGS TEC every 2 hours (deleted every other glat to have 5 deg glat)
if (imod==0)
gpsloc = ' IGS';
  if (if5deg==1)
   DTECnam = sprintf('/hao/aim3/emery/cetic_climatology/igs_tec/long_%3.3d_07325_08020_tec_igs.dat',loncutsE(j))
  else
% ut_beg glat glon1 glon2 tecmed tecmean tecstd wttec averrtec numdays (max 61)
  DTECnam = sprintf('/hao/aim3/emery/dec06/igs/long_%3.3d_06347-9_tec_igs_2p5deg.dat',loncutsE(j))
  jt = 4;
  end
else
% year danum mon mday ut_beg glon1 glon2 glat tecmean 
  jt = 9;
 if (imod==1)
  gpsloc = 'CODE';
  DTECnam = sprintf('/hao/aim3/emery/dec06/igs/cddis/06347-9_long_%3.3d_tec_CODE_2p5deg.dat',loncutsE(j))
 end
 if (imod==2)
  gpsloc = 'ESOC';
  DTECnam = sprintf('/hao/aim3/emery/dec06/igs/cddis/06347-9_long_%3.3d_tec_ESOC_2p5deg.dat',loncutsE(j))
 end
 if (imod==3)
  gpsloc = 'cJPL';
  DTECnam = sprintf('/hao/aim3/emery/dec06/igs/cddis/06347-9_long_%3.3d_tec_cJPL_2p5deg.dat',loncutsE(j))
 end
 if (imod==4)
  gpsloc = 'cUPC';
  DTECnam = sprintf('/hao/aim3/emery/dec06/igs/cddis/06347-9_long_%3.3d_tec_cUPC_2p5deg.dat',loncutsE(j))
 end
 if (imod==5)
  gpsloc = 'cIGS';
  DTECnam = sprintf('/hao/aim3/emery/dec06/igs/cddis/06347-9_long_%3.3d_tec_cIGS_2p5deg.dat',loncutsE(j))
 end
end	% if imod==0 else

DTECin = fopen(DTECnam,'r');
datarea(j) = 0;
gdat(j) = 0;
for nh=1:nt*ndays
 for nl=1:nlat
 linegps = fgetl(DTECin);
 ingps = sscanf(linegps,'%f');
    tecarrmed(nh,nl,j) = ingps(jt);
% Calculate 'global' mean TEC from 24h of 1 glon
% Weight with area of globe and find global gdat values
         datarea(j) = datarea(j) + da(nl);
         gdat(j) = gdat(j) + da(nl)*tecarrmed(nh,nl,j);
 end	% for 1:nlat
end	% for nh=1:nt
  gdat(j) = gdat(j)/datarea(j);
end	% if ifmit==1 else


%% Make plots
gdattit = sprintf(' %5.2f',gdat(j))
            subplot(4,2,j)
	    
    tecarrmed(nt*ndays+1,:,:) = tecarrmed(nt*(ndays-1)+1,:,:);
	    hold on;
            pcolor(mingps,latgps,tecarrmed(:,:,j)')
            caxis([0 50]);		% for Dec 2006
% 2 MArch CANNOT GET gpsloc in if use just title, or gdattit!
  tectit = sprintf('%s GPS TEC av 13-15 Dec 2006, LON=%s %s',gpsloc,num2str(loncutsE(j)),gdattit); 
            title([tectit]);
   
            shading flat
            %xlabel('Date');
            ylabel('Latitude');
%            set(gca,'TickDir','out');
%            datetick('x','HH:MM');
%	    set (gca,'xlim',[0,x2],'Xtick',([0:6:24]),'XTickLabel',{'00:00','06:00','12:00','18:00','24:00'});
	    set (gca,'xlim',[0,x2],'Xtick',([0:24:72]),'XTickLabel',{'347/0UT','348/0UT','349/0UT','350/0UT'},'ylim',[-90 90]);
% Put dots across equator, +/-50 deg
m50 = [-50 -50];
p50 = [50 50];
zero = [0 0];
xtot = [0 x2];
            plot (xtot,zero,':k');
	    plot (xtot,m50,':k');
	    plot (xtot,p50,':k');
% Put dots at 0 UT
	if (ndays>1)
ytot = [-90 90];
	 for n=1:ndays-1
          d0ut = [n*24 n*24];
	  plot (d0ut,ytot,':k');
	 end
	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;
    hc_pos(2)=0.045;
    set(hc,'Position', hc_pos);
hold off;

%end	% for imod=0:4



