clear all;

% May 2012  Program to calculate the median/ave of the TEC values from IGS 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 (72 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

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

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

for lon=1:8

% Get name of IGS GPS TEC file
 file1=sprintf('igstec_2006_lon%3.3d.txt',glon(lon))
 file2=sprintf('igstec_2006_lon%3.3d.txt',glon(lon)+5)
% Open file
IGS1 = fopen(file1,'r');
IGS2 = fopen(file2,'r');
%Skip first header line
 line1 = fgetl(IGS1);
 line2 = fgetl(IGS2);

% Name global TEC output file for daily global means (from 0-22UT - 12 hours every 2 hrs) as func glon
outglob=sprintf('global_tec_igs_06347-9_lon%3.3d-%3.3d.dat',glon(lon),glon(lon)+5)
OUTG = fopen(outglob,'w');

for n=1:ndays
nda = nda0 + n
gtec1(n) = 0.;
gtec2(n) = 0.;

for nh=1:nt
 for nl=1:nlat
 line1 = fgetl(IGS1);
 in1 = sscanf(line1,'%f');
 line2 = fgetl(IGS2);
 in2 = sscanf(line2,'%f');

% year month day_of_month uthr glat glon TEC
% Accumulate in arrays the values of TEC
tec1(lon,n,nh,nl) = in1(7);
tec2(lon,n,nh,nl) = in2(7);
% Calculate global tec (assume always exists and is good)
  gtec1(n) = gtec1(n) + da(nl)*in1(7);
  gtec2(n) = gtec2(n) + da(nl)*in2(7);
end	% nl=1:nlat
end	% nh=1:nt
gtec1(n) = gtec1(n)/area5deg;
gtec2(n) = gtec2(n)/area5deg;
  fprintf(OUTG,'%4.4d %3.3d %2.2d %2.2d %5.2f %5.2f %5.2f %5.2f \n',in1(1),nda,in1(2),in1(3),glon(lon),gtec1(n),glon(lon)+5,gtec2(n));
end	% n=1:ndays


% Calculate the median/ave TEC each day for each longitude
for n=1:ndays
nda = nda0 + n;
% Get name for output file of median TEC values for each glon for each day
outfil=sprintf('long_%3.3d_06%3.3d_tec_igs_2p5deg.dat',glon(lon),nda)
OUT = fopen(outfil,'w');
for nh=1:nt
 for nl=1:nlat
    avtec(lon,nh,nl) = 0.5*(tec1(lon,n,nh,nl)+tec2(lon,n,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
  if nl<33	% allow for minus sign for glats -90 to -5
  fprintf(OUT,'%2.2d %4.1f %5.1f %5.2f \n',(nh-1)*2,lat(nl),glon(lon)+2.5,avtec(lon,nh,nl));
  else
  fprintf(OUT,'%2.2d  %4.1f %5.1f %5.2f \n',(nh-1)*2,lat(nl),glon(lon)+2.5,avtec(lon,nh,nl));
  end
end	% if iprint==1
 end	% nl=1:nlat
end	% nh=1:nt   
end	% n=1:ndays
   
end	% lon=1:8




