clear all;

% 11/11  Program to calculate the median of the median TEC values from MIT GPS satellites
%        from 07325-08020 for 8 longitudes for the Emery et al., Fall AGU 2011 poster 
% 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

ndays = 61;
%ndays = 3;

for lon=1:8

% Get name for output file of median TEC values for each glon
outfile=sprintf('long_%3.3d_07325_08020_tec_mit.dat',glon(lon))
OUT = fopen(outfile,'w');

for n=1:ndays
maxnum = 24*36;
for nh=1:24
 for nl=1:36
  tecmed(lon,n,nh,nl)=-99.;
  tecave(lon,n,nh,nl)=-99.;
  tecerr(lon,n,nh,nl)=-99.;
  numtec(lon,n,nh,nl)=0;
 end
end

% Missing days: Jan 1-3, 2008 (n=42-44)
%if (n<42 | n>44)

% Get name of MIT GPS TEC file
if n<11 
 file=sprintf('mit_gps/long_%d_2007_11_%d.dat',glon(lon),n+20)
end
if (n>10 & n<42)
 file=sprintf('mit_gps/long_%d_2007_12_%d.dat',glon(lon),n-10)
end
if n>41
 file=sprintf('mit_gps/long_%d_2008_1_%d.dat',glon(lon),n-41)
end

% Open file
GPS = fopen(file,'r');

% Read 24h daily file, skipping some glats (-90 to +85)
iend = 0;
glatend = -95.;
icount = 0;

while (icount<maxnum+1) && iend==0
 icount = icount + 1;
 linegps = fgetl(GPS);
% year month day_of_month beghr begmin beg_glat beg_glon end_glon TEC:mean median wted_err std_dev #pts
 ingps = sscanf(linegps,'%f');
% Save last glat for UT=22 assuming this is also the last glat for UT=23
 if ingps(4) == 22
   glatend = ingps(6);
% Exceptions
% 25-30E
   if lon==1 && n==8
    glatend = 75;
   end
   if lon==1 && n==44
    glatend = 75;
   end
   if lon==1 && n>=24 && n<=41
    glatend = 80;
   end
% 90-95E
   if lon==2 && n>=13 && n<=17
    glatend = 70;
   end
   if lon==2 && n>=18 && n<=36
    glatend = 75;
   end
   if lon==2 && n<11
    glatend = 80;
   end
   if lon==2 && n==53
    glatend = 80;
   end
   if lon==2 && n>56
    glatend = 80;
   end
   if lon==2 && n>=39 && n<=41
    glatend = 85;
   end
% 140-145E
   if lon==3 && n<4
    glatend = 60;
   end
   if lon==3 && n>=54 && n<=59
    glatend = 65;
   end
   if lon==3 && n>=4 && n<=5
    glatend = 70;
   end
   if lon==3 && n==41
    glatend = 70;
   end
   if lon==3 && n>=60 && n<=61
    glatend = 70;
   end
   if lon==3 && n>=6 && n<=25
    glatend = 75;
   end
    if lon==3 && n==40
    glatend = 80;
   end
  if lon==3 && n==42
    glatend = 80;
   end
   if lon==3 && n>=44 && n<=53
    glatend = 80;
   end
% 175-180E
   if lon==4 && n<4
    glatend = 70;
   end
   if lon==4 && n>=42 && n<=51
    glatend = 70;
   end
   if lon==4 && n>=53 && n<=56
    glatend = 70;
   end
   if lon==4 && n>=4 && n<=25
    glatend = 75;
   end
   if lon==4 && n>=57 && n<=61
    glatend = 75;
   end
% 200-205E
   if lon==5 && n>=14 && n<=20
    glatend = 75;
   end
   if lon==5 && n>=43 && n<=59
    glatend = 75;
   end
   if lon==5 && n>=1 && n<=13
    glatend = 80;
   end
   if lon==5 && n>=37 && n<=41
    glatend = 85;
   end
% 250-255E
   if lon==6 && n==11
    glatend = 75;
   end
   if lon==6 && n>=4 && n<=6
    glatend = 80;
   end
   if lon==6 && n>=45 && n<=52
    glatend = 80;
   end
   if lon==6 && n>=54 && n<=56
    glatend = 80;
   end
   if lon==6 && n==60
    glatend = 85;
   end
% 285-290E
   if lon==7 && n>=28 && n<=29
    glatend = 80;
   end
   if lon==7 && n>=43 && n<=44
    glatend = 85;
   end
% 345-350E
   if lon==8 && n>=41 && n<=42
    glatend = 75;
   end
   if lon==8 && n==40
    glatend = 80;
   end
   if lon==8 && n>=34 && n<=38
    glatend = 80;
   end
 end
% Accumulate in arrays the values of median TEC ingps(10)
nlat = fix((ingps(6)+95.1)/5.);
tecmed(lon,n,ingps(4)+1,nlat) = ingps(10);
tecave(lon,n,ingps(4)+1,nlat) = ingps(9);
tecerr(lon,n,ingps(4)+1,nlat) = ingps(11);
numtec(lon,n,ingps(4)+1,nlat) = ingps(12);
% Find if this is the last line of the file and end the 'while'
 if ingps(4) == 23 && ingps(6) == glatend
   iend = 1;
 end
end	% while
%end	% missing days n=42-44
end	% n=1:ndays

% Calculate the median TEC over n days for each longitude
for nh=1:24
 for nl=1:36
  clear tecnd;
  clear tecav;
  clear tecer;
  num = 0;
  for n=1:ndays
   if (tecmed(lon,n,nh,nl)>-1.)
    num = num + 1;
    tecnd(num) = tecmed(lon,n,nh,nl);
    tecav(num) = tecave(lon,n,nh,nl)/tecerr(lon,n,nh,nl)^2;
    tecer(num) = 1./tecerr(lon,n,nh,nl)^2;
   end
  end
  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.;
   wtav = 0.;
   errinv = 0.;
   err = 0.;
   for k=1:num
    av = av + tecav(k)/tecer(k);
    err = err + sqrt(1/tecer(k));
    wtav = wtav + tecav(k);
    errinv = errinv + tecer(k);
   end
   av = av/num;
   err = err/num;
   wtavtec(lon,nh,nl) = wtav/errinv;
   avertec(lon,nh,nl) = sqrt(1/errinv);
   stdev = 0.;
   for k=1:num
    stdev = stdev + (av-tecav(k)/tecer(k))^2;
   end
    stdev = sqrt(stdev/num);
  end 
% Print out median tec over ndays
  if nl<19	% allow for minus sign for glats -90 to -5
  fprintf(OUT,'%2.2d %2.2d %3.3d %3.3d %5.2f %5.2f %5.2f %5.2f %5.2f %2.2d\n',nh-1,nl*5-95,ingps(7),ingps(8),medtec(lon,nh,nl),av,stdev,wtavtec(lon,nh,nl),avertec(lon,nh,nl),num);
  else
  fprintf(OUT,'%2.2d  %2.2d %3.3d %3.3d %5.2f %5.2f %5.2f %5.2f %5.2f %2.2d\n',nh-1,nl*5-95,ingps(7),ingps(8),medtec(lon,nh,nl),av,stdev,wtavtec(lon,nh,nl),avertec(lon,nh,nl),num);
  end
 end	% nl=1:36
end	% nh=1:24   
   
end	% lon=1:8




