clear all;

% Mar 2012  Program to calculate the medians of the SAMI3 TEC, N(h)mF2 values 
%        from 07325-08020 for 8 longitudes for the Emery et al., SWW April 2012 talk 
% SAMI3 from Sarah McDonald has 37 values every 5 deg from -90 to +90 (ave 2 together to get 36 glats -87.5 to +87.5!)
%   year month day UTh (0.5 to 23.5 UT) glat glon TEC(TECU) NmF2(e+5 m-3) hmF2(km) where missing are 0 at hi lats
%     (ignore NmF2 if hmF2=0; have TEC=0 when hmF2=0!)
% DO NOT Read dailyTEC (ave, median, and std): year mo da glat glon av med std(daily TEC) with missing=0. at hi lats but calculate it!

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

% 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 = 24;
nlat = 37;
dlat = 5.;
dlon = 5.;
% calculate daily 'global' tec for each 5 deg longitude slice

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


for lon=1:8

% Get name for output file of median TEC, NmF2, hmF2 values for each glon
outfil=sprintf('long_%3.3d_07325_08020_Ne_sami3.dat',glon(lon))
outfil1=sprintf('long_%3.3d_07325_08020_tec_sami3_hi.dat',glon(lon))
outfil2=sprintf('long_%3.3d_07325_08020_tec_sami3_lo.dat',glon(lon))
OUT = fopen(outfil,'w');
OUT1 = fopen(outfil1,'w');
OUT2 = fopen(outfil2,'w');

% Get names of SAMI3 files
 file1=sprintf('samitec_2007_lon%3.3d.txt',glon(lon))
 file2=sprintf('samitec_2008_lon%3.3d.txt',glon(lon))

% Open files
IGS1 = fopen(file1,'r');
IGS2 = fopen(file2,'r');
%Skip first header line
 line = fgetl(IGS1);
 line = fgetl(IGS2);

% Calculate the 'global' TEC for each of 8 glons
outglob=sprintf('global_tec_sami3_07325_08020_lon%3.3d.dat',glon(lon))
OUTG = fopen(outglob,'w');

for n=1:ndays
nda = nda0 + n;
if (nda>365)
 nda = nda - 365;
end
gtec(n) = 0.;
area(n) = 0.;
gnmf2(n) = 0.;
ghmf2(n) = 0.;

for nh=1:nt
 if (nda>nda0)
  line = fgetl(IGS1);
  in = sscanf(line,'%f');
 else
  line = fgetl(IGS2);
  in = sscanf(line,'%f');
 end
% Save southern most glat
  tecs=in(7);
  nmf2s = in(8);
  hmf2s = in(9);

 for nl=1:nlat-1
 if (nda>nda0)
  line = fgetl(IGS1);
  in = sscanf(line,'%f');
 else
  line = fgetl(IGS2);
  in = sscanf(line,'%f');
 end

% year month day_of_month uthr glat glon TEC NmF2 hmF2
% Accumulate in arrays the values of TEC
 if (in(7)>0.1 & tecs>0.1)
   tec(lon,n,nh,nl) = 0.5*(tecs+in(7));
   nmf2(lon,n,nh,nl) = 0.5*(nmf2s+in(8));
   hmf2(lon,n,nh,nl) = 0.5*(hmf2s+in(9));
 else
   tec(lon,n,nh,nl) = 0.;
   nmf2(lon,n,nh,nl) = 0.;
   hmf2(lon,n,nh,nl) = 0.;
   if (in(7)>0.1)
     tec(lon,n,nh,nl) = in(7);
     nmf2(lon,n,nh,nl) = in(8);
     hmf2(lon,n,nh,nl) = in(9);
   end
   if (tecs>0.1)
     tec(lon,n,nh,nl) = tecs;
     nmf2(lon,n,nh,nl) = nmf2s;
     hmf2(lon,n,nh,nl) = hmf2s;
   end
  end
  if (tec(lon,n,nh,nl)>0.01)
% Calculate global tec 
   gtec(n) = gtec(n) + da(nl)*tec(lon,n,nh,nl);
   gnmf2(n) = gnmf2(n) + da(nl)*nmf2(lon,n,nh,nl);
   ghmf2(n) = ghmf2(n) + da(nl)*hmf2(lon,n,nh,nl);
   area(n) = area(n) + da(nl);
  end
% Save southern most glat
  tecs=in(7);
  nmf2s = in(8);
  hmf2s = in(9);
end	% nl=1:nlat-1
end	% nh=1:nt
gtec(n) = gtec(n)/area(n);
gnmf2(n) = gnmf2(n)/area(n);
ghmf2(n) = ghmf2(n)/area(n);


  fprintf(OUTG,'%4.4d %3.3d %3.3d %5.3f %5.2f %8.3e %6.2f\n',in(1),nda,glon(lon),area(n)/area5deg,gtec(n),gnmf2(n),ghmf2(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-1
  clear tecnd;
  clear tecav;
  clear tecer;
  clear nmf2nd;
  clear nmf2av;
  clear nmf2er;
  clear hmf2nd;
  clear hmf2av;
  clear hmf2er;
  num = 0;
  for n=1:ndays
   if (n12==3 & tec(lon,n,nh,nl)>0.01)
    num = num + 1;
    tecnd(num) = tec(lon,n,nh,nl);
    nmf2nd(num) = nmf2(lon,n,nh,nl);
    hmf2nd(num) = hmf2(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 & tec(lon,n,nh,nl)>0.01)
    num = num + 1;
    tecnd(num) = tec(lon,n,nh,nl);
    nmf2nd(num) = nmf2(lon,n,nh,nl);
    hmf2nd(num) = hmf2(lon,n,nh,nl);
    end
   end	% nn=1:25
  end	% n12<3

  end	% n=1:ndays
  medtec(lon,nh,nl) = -1.;
  mednmf2(lon,nh,nl) = -1.;
  medhmf2(lon,nh,nl) = -1.;
  av = -1.;
  stdev = -1.;
  avn = -1.;
  stdn = -1.;
  avh = -1.;
  stdh = -1.;
  if (num > 0)
   medtec(lon,nh,nl) = median(tecnd);
   mednmf2(lon,nh,nl) = median(nmf2nd);
   medhmf2(lon,nh,nl) = median(hmf2nd);
   av = 0.;
   avn = 0.;
   avh = 0.;
   for k=1:num
    av = av + tecnd(k);
    avn = avn + nmf2nd(k);
    avh = avh + hmf2nd(k);
   end
   av = av/num;
   avn = avn/num;
   avh = avh/num;
   stdev = 0.;
   stdn = 0.;
   stdh = 0.;
   for k=1:num
    stdev = stdev + (av-tecnd(k))^2;
    stdn = stdn + (avn-nmf2nd(k))^2;
    stdh = stdh + (avh-hmf2nd(k))^2;
   end
    stdev = sqrt(stdev/num);
    stdn = sqrt(stdn/num);
    stdh = sqrt(stdh/num);
  end 

if (n12==1)
% Print out median tec over ndays
  fprintf(OUT1,'%2.2d %5.1f %3.3d %3.3d %5.2f %5.2f %5.2f %8.3e %8.3e %8.3e %6.2f %6.2f %6.2f %2.2d\n',nh-1,lat(nl),glon(lon),glon(lon)+5,medtec(lon,nh,nl),av,stdev,mednmf2(lon,nh,nl),avn,stdn,medhmf2(lon,nh,nl),avh,stdh,num);
end
if (n12==2)
% Print out median tec over ndays
  fprintf(OUT2,'%2.2d %5.1f %3.3d %3.3d %5.2f %5.2f %5.2f %8.3e %8.3e %8.3e %6.2f %6.2f %6.2f %2.2d\n',nh-1,lat(nl),glon(lon),glon(lon)+5,medtec(lon,nh,nl),av,stdev,mednmf2(lon,nh,nl),avn,stdn,medhmf2(lon,nh,nl),avh,stdh,num);
end
if (n12==3)
% Print out median tec over ndays
  fprintf(OUT,'%2.2d %5.1f %3.3d %3.3d %5.2f %5.2f %5.2f %8.3e %8.3e %8.3e %6.2f %6.2f %6.2f %2.2d\n',nh-1,lat(nl),glon(lon),glon(lon)+5,medtec(lon,nh,nl),av,stdev,mednmf2(lon,nh,nl),avn,stdn,medhmf2(lon,nh,nl),avh,stdh,num);
end
 end	% nl=1:nlat
end	% nh=1:nt 
end	% n12=1:3  
   
end	% lon=1:8




