%%
%% This program is written for plotting the mean TEC data at 3,15LT for all models at 1 of 8 glons
%% BA Emery, May 2012 from /aim3/emery/dec06/igs/plot_lt.m
clear all; close all

% imod = (0)1-5 for (NRL IGS) CODE, ESOC, JPL, UPC, and IGS 2h data
% imod = 6,7 for MIT,JPL 15-min data

ndays = 3;

% Read 71 or 36 latitudes, but only save every other of 71 lats for plotting
dlat = 5;		% 36*5deg=180 deg from 90S-90N		
%lat5deg=-90:dlat:85;
lat5deg = -87.5:5:87.5;
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;

dlon = 5;	% dlon = 360/24=15deg for hourly; dlon=15/4=2.75 for 15 min
nlon = 360/dlon;


% Data
        numtec=NaN(length(lat5deg),6,7);
        tecarrmed=NaN(length(lat5deg),6,7);

 vartit = ' TEC '
% Data
 loncutsE=[25,90,140,175,200,250,285,345]; % specify climatology longitudes to plot
for jlon = 1:8

% Calculate UT n for 3 and 15 LT for 15 min and 2h files (LT=UT+glon/15)
ut03 = mod(3-loncutsE(jlon)/15,24)
ut15 = mod(15-loncutsE(jlon)/15,24)
n03_15m = round(ut03*4)+1
n03_2h = round(ut03/2)+1
n15_15m = n03_15m + 12*4;
irev = 0;
if (n15_15m > 96)
 n15_15m = n15_15m - 96;
 irev = 1;
end
n15_2h = n03_2h + 6;
if (n15_2h > 12)
 n15_2h = n15_2h - 12;
 irev = 1;
end


% Read MIT and JPL 15min, 5deg data first
 nt = 96;
 for imod=6:7
  if (imod==6)
   gpsloc=' MIT';
   DTECnam = sprintf('/hao/aim3/emery/dec06/igs/TEC_lon_%3.3d_MIT_2006.347-2006.349.dat',loncutsE(jlon))
  else
   gpsloc=' JPL';
   DTECnam = sprintf('/hao/aim3/emery/dec06/igs/TEC_lon_%3.3d_JPL_2006.347-2006.349.dat',loncutsE(jlon))
  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);
m6 = 0;
for nh=1:nt*ndays
nskp = 1;
if (nh==n03_15m | nh==n15_15m | nh==n03_15m+nt | nh==n15_15m+nt | nh==n03_15m+2*nt | nh==n15_15m+2*nt)
 nskp = 0;
 m6 = m6 + 1;
 if (irev==0)
 n6 = m6;
 else
  if (round((m6-1)/2.)*2==m6-1)
   n6 = m6 + 1;
  else
   n6 = m6 - 1;
  end
 end
end
 for nl=1:36
 linegps = fgetl(DTECin);
 ingps = sscanf(linegps,'%f');
 if (nskp==0)
  numtec(nl,n6,imod) = ingps(12);
  if (numtec(nl,n6,imod) > 0)
% 8=tecavwt and 9=tecmed
   if (ingps(9) > 0.) 
    tecarrmed(nl,n6,imod) = ingps(9);
   end
  else
    numtec(nl,n6,imod) = 0;
  end
 end	% if nskp==0
 end	% for 1:nlat
end	% for nh=1:nt*ndays
end	% for imod=6:7


% Read IGS TEC every 2 hours (deleted every other glat to have 5 deg glat) at 2.5 deg but save 5deg
for imod=0:4
if (imod==0)
gpsloc = ' IGS';
% 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(jlon))
  jt = 4;
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(jlon))
 end
 if (imod==2)
  gpsloc = 'ESOC';
  DTECnam = sprintf('/hao/aim3/emery/dec06/igs/cddis/06347-9_long_%3.3d_tec_ESOC_2p5deg.dat',loncutsE(jlon))
 end
 if (imod==3)
  gpsloc = 'cJPL';
  DTECnam = sprintf('/hao/aim3/emery/dec06/igs/cddis/06347-9_long_%3.3d_tec_cJPL_2p5deg.dat',loncutsE(jlon))
 end
 if (imod==4)
  gpsloc = 'cUPC';
  DTECnam = sprintf('/hao/aim3/emery/dec06/igs/cddis/06347-9_long_%3.3d_tec_cUPC_2p5deg.dat',loncutsE(jlon))
 end
 if (imod==5)
  gpsloc = 'cIGS';
  DTECnam = sprintf('/hao/aim3/emery/dec06/igs/cddis/06347-9_long_%3.3d_tec_cIGS_2p5deg.dat',loncutsE(jlon))
 end
end	% if imod==0 else

DTECin = fopen(DTECnam,'r');
nt = 12;
m6 = 0;
for nh=1:nt*ndays
nskp = 1;
if (nh==n03_2h | nh==n15_2h | nh==n03_2h+nt | nh==n15_2h+nt | nh==n03_2h+2*nt | nh==n15_2h+2*nt)
 nskp = 0;
 m6 = m6 + 1;
 if irev==0
  n6 = m6;
 else
  if (round((m6-1)/2.)*2==m6-1)
   n6 = m6 + 1;
  else
   n6 = m6 - 1;
  end
 end
end
 nli = 0;
 for nl=1:71
 linegps = fgetl(DTECin);
 ingps = sscanf(linegps,'%f');
   if nskp==0 & round((nl-1)/2.)*2==nl-1
    nli = nli + 1;
    tecarrmed(nli,n6,imod+1) = ingps(jt);
   end
 end	% for 1:nlat
end	% for nh=1:nt

end	% for imod=0:4


%% Make plots of 3LT on left column and 15LT on right column for 3 days

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

% tecx=50-70 TECU 
if (jlon>3 & jlon<7)
 tecx = 70;
end
if (jlon<4)
 tecx = 50;
end
if (jlon>6)
 tecx = 60;
end
titx = 0.93*tecx;
titx2 = 0.85*tecx;

for j=1:6
            subplot(3,2,j)   
        hold on;
 if j==1
  tectit = sprintf('GPS TEC av 03 LT 13 Dec 2006, LON=%s',num2str(loncutsE(jlon))); 
 end
 if j==2
  tectit = sprintf('GPS TEC av 15 LT 13 Dec 2006, LON=%s',num2str(loncutsE(jlon))); 
 end
 if j==3
  tectit = sprintf('GPS TEC av 03 LT 14 Dec 2006, LON=%s',num2str(loncutsE(jlon))); 
 end
 if j==4
  tectit = sprintf('GPS TEC av 15 LT 14 Dec 2006, LON=%s',num2str(loncutsE(jlon))); 
 end
 if j==5
  tectit = sprintf('GPS TEC av 03 LT 15 Dec 2006, LON=%s',num2str(loncutsE(jlon))); 
 end
 if j==6
  tectit = sprintf('GPS TEC av 15 LT 15 Dec 2006, LON=%s',num2str(loncutsE(jlon))); 
 end
            title([tectit]);
            
            xlabel('latitude');
            ylabel('TECU');            
	    set (gca,'xlim',[-90 90],'ylim',[0 tecx]);
            plot (lat5deg,tecarrmed(:,j,2),'-r');
            plot (lat5deg,tecarrmed(:,j,3),'-m');
            plot (lat5deg,tecarrmed(:,j,5),'-y');      
            plot (lat5deg,tecarrmed(:,j,1),'-g');
            plot (lat5deg,tecarrmed(:,j,6),'-xk');
            plot (lat5deg,tecarrmed(:,j,7),'-ob');
            plot (lat5deg,tecarrmed(:,j,4),'--b');


	text (-85,titx,'MIT','Color','k');
	text (-70,titx,'JPL,--CDDIS','Color','b');
	text (35,titx,'CODE','Color','r');
	text (65,titx,'ESOC','Color','m');
  if (round((j-1)/2)*2==j-1)
	text (-15,titx,'IGS','Color','g');
	text (10,titx,'UPC','Color','y');
  else
	text (-85,titx2,'IGS','Color','g');
	text (65,titx2,'UPC','Color','y');
  end
    end				% end of plotting 6 LTs

hold off;

end		% for jlon=1:8



