clear all;

% Mar 2012  Program to calculate the 2 medians of the median TEC values from IGS GPS satellites
%        from 07325-08020 for 8 longitudes for the Emery et al., SEA13 Mar 2012 talk 
%  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];
  
%  07325-08020 is 2007_11_21(to30), 2007_12_1(to31), 2008_1_1(to20) for file names - cat first!
%  Watch for missing latitudes in GPS files

% There are 25 days in each median for #1 (Kp>=2, Vsw>=500) and #2 (Kp<=1, Vsw<=450)
% 25 HSS days: 07325-330, 345-346, 351-357, 370-373, 379-384
% 25 Low Kp days: 07332-343, 349-350, 358-360, 362-364, 366-369, 376
nd1hi = [1,2,3,4,5,6,21,22,27,28,29,30,31,32,33,46,47,48,49,55,56,57,58,59,60];
nd2lo = [8,9,10,11,12,13,14,15,16,17,18,19,25,26,34,35,36,38,39,40,42,43,44,45,52];

ndays = 61;
%ndays = 3;
nda0 = 324;

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 for output file of median TEC values for each glon
%outfil=sprintf('long_%3.3d_07325_08020_tec_igs.dat',glon(lon))
%outfil1=sprintf('long_%3.3d_07325_08020_tec_igs_hi.dat',glon(lon))
%outfil2=sprintf('long_%3.3d_07325_08020_tec_igs_lo.dat',glon(lon))
outfil=sprintf('long_%3.3d_07325_08020_tec_igs_2p5deg.dat',glon(lon))
outfil1=sprintf('long_%3.3d_07325_08020_tec_igs_hi_2p5deg.dat',glon(lon))
outfil2=sprintf('long_%3.3d_07325_08020_tec_igs_lo_2p5deg.dat',glon(lon))
OUT = fopen(outfil,'w');
OUT1 = fopen(outfil1,'w');
OUT2 = fopen(outfil2,'w');

% Get name of IGS GPS TEC file
 file1=sprintf('igs_tec/igstec_07325_08020_lon%3.3d.txt',glon(lon))
 file2=sprintf('igs_tec/igstec_07325_08020_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_07325_08020_lon%3.3d-%3.3d.dat',glon(lon),glon(lon)+5)
OUTG = fopen(outglob,'w');

for n=1:ndays
nda = nda0 + n
if (nda>365)
 nda = nda - 365
end
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 TEC over n days for each longitude
for n12=1:3
for nh=1:nt
 for nl=1:nlat
  clear tecnd;
  clear tecav;
  clear tecer;
  num = 0;
  for n=1:ndays
   if (n12==3)
    num = num + 1;
    tecnd(num) = 0.5*(tec1(lon,n,nh,nl)+tec2(lon,n,nh,nl));
   end
   if (n12<3)
   for nn=1:25
    n25 = nd1hi(nn);
    if (n12==2)
     n25 = nd2lo(nn);
    end
    if (n==n25)
    num = num + 1;
    tecnd(num) = 0.5*(tec1(lon,n,nh,nl)+tec2(lon,n,nh,nl));
    end
   end	% nn=1:25
  end	% n12<3

  end	% n=1:ndays
  medtec(lon,nh,nl) = -1.;
  wtavtec(lon,nh,nl) = -1.;
  avertec(lon,nh,nl) = -1.;
  av = -1.;
  stdev = -1.;
  if (num > 0)
   medtec(lon,nh,nl) = median(tecnd);
   av = 0.;
   for k=1:num
    av = av + tecnd(k);
   end
   av = av/num;
   stdev = 0.;
   for k=1:num
    stdev = stdev + (av-tecnd(k))^2;
   end
    stdev = sqrt(stdev/num);
  end 
% 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
if (n12==1)
% Print out median tec over ndays
  if nl<33	% allow for minus sign for glats -90 to -5
  fprintf(OUT1,'%2.2d %4.1f %5.1f %5.2f %5.2f %5.2f %2.2d\n',(nh-1)*2,lat(nl),glon(lon)+2.5,medtec(lon,nh,nl),av,stdev,num);
  else
  fprintf(OUT1,'%2.2d  %4.1f %5.1f %5.2f %5.2f %5.2f %2.2d\n',(nh-1)*2,lat(nl),glon(lon)+2.5,medtec(lon,nh,nl),av,stdev,num);
  end
end
if (n12==2)
% Print out median tec over ndays
  if nl<33	% allow for minus sign for glats -90 to -5
  fprintf(OUT2,'%2.2d %4.1f %5.1f %5.2f %5.2f %5.2f %2.2d\n',(nh-1)*2,lat(nl),glon(lon)+2.5,medtec(lon,nh,nl),av,stdev,num);
  else
  fprintf(OUT2,'%2.2d  %4.1f %5.1f %5.2f %5.2f %5.2f %2.2d\n',(nh-1)*2,lat(nl),glon(lon)+2.5,medtec(lon,nh,nl),av,stdev,num);
  end
end
if (n12==3)
% 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 %5.2f %5.2f %2.2d\n',(nh-1)*2,lat(nl),glon(lon)+2.5,medtec(lon,nh,nl),av,stdev,num);
  else
  fprintf(OUT,'%2.2d  %4.1f %5.1f %5.2f %5.2f %5.2f %2.2d\n',(nh-1)*2,lat(nl),glon(lon)+2.5,medtec(lon,nh,nl),av,stdev,num);
  end
end
end	% if iprint==1
 end	% nl=1:nlat
end	% nh=1:nt 
end	% n12=1:3  
   
end	% lon=1:8




