clear all;

% May 2012  Program to calculate the median/ave of the TEC values from CDDIS GPS satellites
%        from 06347-9 for 8 longitudes for the Emery et al. June 2012 MSU talk and CEDAR poster
%  IGS interpolates data so there are no missing values every 2.5 deg from -87.5 to +87.5 (71 glats)
%       Grid points every 5 deg glon so average 25+30E for 5 deg from 24-30E glon every 2 UT hours
%   year month day UTh glat glon TEC
%  Calculate also the 'global TEC' for each longitude band

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

if imod==1
 gpsmod='CODE'
 modn='cod';
 iskp(1)=348;
 iskp(2)=347;
 iskp(3)=342;
end
if imod==2
 gpsmod='ESOC'
 modn='esa';
 iskp(1)=293;
 iskp(2)=293;
 iskp(3)=294;
end
if imod==3
 gpsmod='cJPL'
 modn='jpl';
 iskp(1)=259;
 iskp(2)=259;
 iskp(3)=260;
end
if imod==4
 gpsmod='cUPC'
 modn='upc';
 iskp(1)=145;
 iskp(2)=144;
 iskp(3)=142;
end
if imod==5
 gpsmod='cIGS'
 modn='igs';
 iskp(1)=392;
 iskp(2)=396;
 iskp(3)=390;
end
 

% 8 begin lons (for 5 deg glon): 25,90,140,175,200,250,285,345
glon=[25,90,140,175,200,250,285,345];

% Number of begin glon from 73 glons from -180 to +180
ngbeg=[42,55,65,72 ,5  ,15 ,22 ,34 ];

ndays = 3;
nda0 = 346;

nt = 12;
nlat = 71;
% calculate daily 'global' tec for each 5 deg longitude slice

re=6371.e5;   %  earth radius in cm
dtr = pi/180.;
% Good for dres 2.5x5.0 deg only!
area5deg = 0;
for ii=1:nlat
 lat(ii) = -90.0 + ii*2.5;
 da(ii)=re*cos(lat(ii)*dtr)*(2.5*dtr)*(re*5.0*dtr);
 area5deg = area5deg+da(ii)*nt;
end
area5deg

% Name global TEC output file for daily global means (from 0-22UT - 12 hours every 2 hrs) 
outglob=sprintf('global_tec_%s_06347-9.dat',gpsmod)
OUTG = fopen(outglob,'w');

for nd=1:ndays
nda = nda0 + nd
 for lon=1:8
  gtec1(lon) = 0.;
  gtec2(lon) = 0.;
 end
% Get name of CDDIS GPS TEC file
 file=sprintf('%sg%3.3d0.06i',modn,nda)
% Open file
IGS = fopen(file,'r');
%Skip header lines until 'END OF HEADER'
 for isk=1:iskp(nd)
  line = fgetl(IGS);
 end

for nh=1:nt
% Skip over END OF (TEC MAP mostly), START of TEC MAP, then read: year mon day UTh,m,s
 line=fgetl(IGS);
 line=fgetl(IGS);
 line=fgetl(IGS);
 indate=sscanf(line,'%f');
 year = indate(1);
 mon = indate(2);
 dmon = indate(3);
 uthbeg = indate(4);
% CDDIS format from +87.5 to -87.5 so have to invert lats!
 for nl1=1:nlat
  nl = nlat-nl1+1;
 line = fgetl(IGS);
 inlat = sscanf(line,'%f');
% Read TEC from -180 to +180 (where #1=#73) 16 lons at a time
 line = fgetl(IGS);
 intec1 = sscanf(line,'%f');
 line = fgetl(IGS);
 intec2 = sscanf(line,'%f');
 line = fgetl(IGS);
 intec3 = sscanf(line,'%f');
 line = fgetl(IGS);
 intec4 = sscanf(line,'%f');
 line = fgetl(IGS);
 intec5 = sscanf(line,'%f');
 for n=1:16
  intec(n) = intec1(n);
  intec(n+16) = intec2(n);
  intec(n+32) = intec3(n);
  intec(n+48) = intec4(n);
 end
 for n=1:9
  intec(n+64) = intec5(n);
 end
%stop

for lon=1:8
% year month day_of_month uthr glat glon TEC
% Accumulate in arrays the values of TEC (read in 0.1 TECU)
tec1(lon,nd,nh,nl) = intec(ngbeg(lon))*0.1;
tec2(lon,nd,nh,nl) = intec(ngbeg(lon)+1)*0.1;

% Calculate global tec (assume always exists and is good)
  gtec1(lon) = gtec1(lon) + da(nl)*tec1(lon,nd,nh,nl);
  gtec2(lon) = gtec2(lon) + da(nl)*tec2(lon,nd,nh,nl);
end	% lon=1:8
end	% nl1=1:nlat
end	% nh=1:nt

for lon=1:8
gtec1(lon) = gtec1(lon)/area5deg;
gtec2(lon) = gtec2(lon)/area5deg;

% Get name for output file of median TEC values for each glon for each day
outfil=sprintf('long_%3.3d_06%3.3d_tec_%s_2p5deg.dat',glon(lon),nda,gpsmod)
OUT = fopen(outfil,'w');
for nh=1:nt
 for nl=1:nlat
    avtec(lon,nh,nl) = 0.5*(tec1(lon,nd,nh,nl)+tec2(lon,nd,nh,nl));
% Print every other latitude so every 5 deg glat (ie, nl=1,3,..71)
  nl2 = fix(nl/2.);
  nl0 = nl2*2;
  iprint = 1;
  if nl0==nl
   iprint=0;
  end
% TEMP print all lats
  iprint = 1;
if iprint==1
% Print out median tec over ndays
  fprintf(OUT,'%4.4d %3.3d %2.2d %2.2d %2.2d %5.2f %5.2f %4.1f %5.2f \n',year,nda,mon,dmon,(nh-1)*2,glon(lon),glon(lon)+5,lat(nl),avtec(lon,nh,nl));
end	% if iprint==1
 end	% nl=1:nlat
end	% nh=1:nt   
end 	% lon=1:8
  fprintf(OUTG,'%4.4d %3.3d %2.2d %2.2d %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f\n',year,nda,mon,dmon,glon(1)+2.5,0.5*(gtec1(1)+gtec2(1)),glon(2)+2.5,0.5*(gtec1(2)+gtec2(2)),glon(3)+2.5,0.5*(gtec1(3)+gtec2(3)), ...
glon(4)+2.5,0.5*(gtec1(4)+gtec2(4)),glon(5)+2.5,0.5*(gtec1(5)+gtec2(5)),glon(6)+2.5,0.5*(gtec1(6)+gtec2(6)),glon(7)+2.5,0.5*(gtec1(7)+gtec2(7)),glon(8)+2.5,0.5*(gtec1(8)+gtec2(8)));
end	% nd=1:ndays





