clear all;

% Dec 2011  Program to calculate the 2 medians of the median TEC values from MIT GPS satellites
%        from 07325-08020 for 8 longitudes for the Emery et al., SEA13 Mar 2012 talk 
% 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];

% IRI TEC 20,200km ONLY in LT not UT! correct it in median! LT=0:23, UT=LT-glon/15=LT+corrh;
corrh = [-2,-6,-9,-12,-13,-17,-19,-23];
%ich = 1;
ich = 0;
%  12 Dec 2011:  IRI TEC now in UT and up to 500 km (low) and 1500 km so have to read 2 files
%ird2 = 1;
ird2 = 0;

% Only TIEGCM W05 dres TIMED lb has 0 UT meaning 0 UT.  Others are interpolated 
%   so 0 UT means 0030 UT.

%  The model longitudes are 25 deg wide like COSMIC - not like MIT GPS TEC
% lon_LLL_1_IRI(CTIPE,TIE-GCM)_2007.325-2008.020.txt where LLL=glon-10

% day 42 (Jan 1, 2008) is bad for 0 UT dres W05 TIMED lb
ndays = 61;
%ndays = 40
n36 = 36;
beglat36=-90:5:85;
beglat71=-87.5:2.5:87.5;
% TIEGCM-WT are dres with 72 lats
n36 = 72

% if5deg=1,0 if do,not use 5deg grid out output (can only be 0 if n36=72 aslo, so output will be 71 lats for IGS)
if5deg = 0
if (n36==36 & if5deg==0)
 if5deg==1
end
if (if5deg==0)
 nl36 = 71;
 beglat = beglat71;
else
 nl36 = 36;
 beglat = beglat36;
end


for lon=1:8
glon1 = glon(lon)-10;
glon2 = glon(lon)+15;
if (n36==72)
% NOT QUITE CORRECT (2.5deg), but close enough in integer form
 glon1 = glon(lon)+1;
 glon2 = glon(lon)+4;	
end

% Get name for output file of median TEC values for each glon
%outfil1=sprintf('long_%3.3d_07325_08020_1_IRI_1500km_hi.dat',glon(lon))
%outfil2=sprintf('long_%3.3d_07325_08020_1_IRI_1500km_lo.dat',glon(lon))
%outfile=sprintf('long_%3.3d_07325_08020_1_IRI.dat',glon(lon))
%outfile=sprintf('long_%3.3d_07325_08020_1_CTIPE.dat',glon(lon))
%outfil1=sprintf('long_%3.3d_07325_08020_1_CTIPE_hi.dat',glon(lon))
%outfil2=sprintf('long_%3.3d_07325_08020_1_CTIPE_lo.dat',glon(lon))
%outfile=sprintf('long_%3.3d_07325_08020_1_TIE-GCM.dat',glon(lon))
%outfil1=sprintf('long_%3.3d_07325_08020_1_TIE-GCM_hi.dat',glon(lon))
%outfil2=sprintf('long_%3.3d_07325_08020_1_TIE-GCM_lo.dat',glon(lon))
if (if5deg==1)
outfile=sprintf('long_%3.3d_07325_08020_2_TIE-GCM_alld.dat',glon(lon))
outfil1=sprintf('long_%3.3d_07325_08020_2_TIE-GCM_hid.dat',glon(lon))
outfil2=sprintf('long_%3.3d_07325_08020_2_TIE-GCM_lod.dat',glon(lon))
else
% Have double the number of lats (36*2=72) - or 71 lats to correspond with IGS
outfile=sprintf('long_%3.3d_07325_08020_2_TIE-GCM_all_2p5deg.dat',glon(lon))
outfil1=sprintf('long_%3.3d_07325_08020_2_TIE-GCM_hi_2p5deg.dat',glon(lon))
outfil2=sprintf('long_%3.3d_07325_08020_2_TIE-GCM_lo_2p5deg.dat',glon(lon))
end
OUT = fopen(outfile,'w');
OUT1 = fopen(outfil1,'w');
OUT2 = fopen(outfil2,'w');

% Get name of 61-day output file (24h*61d*36glats=52704 lines each)
%infile=sprintf('ccmc_output/lon_%3.3d_1_IRI_2007.325-2008.020.txt',glon1)
%infilt=sprintf('ccmc_output/IRI_TEC_UT/TEC_lon_%3.3d_1_IRI_2007.325-2008.020.1500.txt',glon1)
%infile=sprintf('ccmc_output/lon_%3.3d_1_CTIPE_2007.325-2008.020.txt',glon1)
%infile=sprintf('ccmc_output/lon_%3.3d_1_TIE-GCM_2007.325-2008.020.txt',glon1)
%%infile=sprintf('ccmc_output/lon_%3.3d_2_TIE-GCM_2007.325-2008.026c.txt',glon(lon))
%infile=sprintf('ccmc_output/lon_%3.3d_07325_08026_2_TIE-GCMc.dat',glon(lon))
infile=sprintf('ccmc_output/lon_%3.3d_2_TIE-GCM_2007.325-2008.026d.txt',glon(lon))

% Open file
INP = fopen(infile,'r');
if (ird2==1)
 INPT = fopen(infilt,'r');
end

for n=1:ndays
for mh=1:24
nh = mh;
% Go from LT to UT
%if ich==1
% nh = mh + corrh(lon);
% if nh<1
%  nh = nh + 24;
% end
%end
for nlat=1:n36
 line = fgetl(INP);
% year daynum beghr begmin begsec beg_glat beg_glon TEC(TECU) NmF2(cm-5?) hmF2(km)
% dres TIE adds wiav(m/s) on end
 input = sscanf(line,'%f');
if (ird2==1)
 lint = fgetl(INPT);
 inputt = sscanf(lint,'%f');
end

if (n36==36)
% Accumulate in arrays the values of TEC, NmF2, hmF2
tec(lon,n,nh,nlat) = input(8);
nmf2(lon,n,nh,nlat) = input(9);
hmf2(lon,n,nh,nlat) = input(10);
 if (ird2==1)
  tec(lon,n,nh,nlat) = inputt(8);
 end
else
% 2_TIE-GCM is dres so have to average 2 pts for single 5deg lat bin
tecb(nh,nlat) = input(8);
% Must correct nmf2 in 2_TIE-GCM by 1.e+6 before c.txt files
%nmf2b(nh,nlat) = input(9)*1.e+6;
nmf2b(nh,nlat) = input(9);
hmf2b(nh,nlat) = input(10);
wiavb(nh,nlat) = input(11);
wiavlo(nh,nlat) = input(12);
uiavlo(nh,nlat) = input(13);
end
% Model TEC was never negative
if (input(8) < 0.)
 input
 stop
end
% TIEGCM NmF2 was sometimes negative, but often hmF2 in low-lats <100km so skip
%if (input(9) < 0. | input(10)<100. | input(10)>800.)
%input
%end
end	% nlat=1:36
end	% nh=1:24

% Now combine 2 dres lat bins together (means have 55 to 60 from 53.75 to 58.75)
% Therefore, do NOT combine, but pick an exact latitude like 57.5 (ie, nl+1)!
if (n36==72)
if (if5deg==1)
for nh=1:24
for nl=1:2:72
nlat = (nl+1)/2;
%tec(lon,n,nh,nlat) = 0.5*(tecb(nh,nl)+tecb(nh,nl+1));
%wiav(lon,n,nh,nlat) = 0.5*(wiavb(nh,nl)+wiavb(nh,nl+1));
%if (hmf2b(nh,nl)>100. & hmf2b(nh,nl) <800. & hmf2b(nh,nl+1)>100. & hmf2b(nh,nl+1)<800.)
% nmf2(lon,n,nh,nlat) = 0.5*(nmf2b(nh,nl)+nmf2b(nh,nl+1));
% hmf2(lon,n,nh,nlat) = 0.5*(hmf2b(nh,nl)+hmf2b(nh,nl+1));
%else
%  nmf2(lon,n,nh,nlat) = -1.;
%  hmf2(lon,n,nh,nlat) = 0.; 
% if (hmf2b(nh,nl)>100. & hmf2b(nh,nl) <800.)
%  nmf2(lon,n,nh,nlat) = nmf2b(nh,nl);
%  hmf2(lon,n,nh,nlat) = hmf2b(nh,nl);
% end
% if (hmf2b(nh,nl+1)>100. & hmf2b(nh,nl+1) <800.)
%  nmf2(lon,n,nh,nlat) = nmf2b(nh,nl+1);
%  hmf2(lon,n,nh,nlat) = hmf2b(nh,nl+1);
% end
%end	% if,else
% Choose nl+1 latitude (center of 5deg bin in data)
tec(lon,n,nh,nlat) = tecb(nh,nl+1);
wiav(lon,n,nh,nlat) = wiavb(nh,nl+1);
wiavl(lon,n,nh,nlat) = wiavlo(nh,nl+1);
uiav(lon,n,nh,nlat) = uiavlo(nh,nl+1);
if (hmf2b(nh,nl+1)>100. & hmf2b(nh,nl+1) <800.)
  nmf2(lon,n,nh,nlat) = nmf2b(nh,nl+1);
  hmf2(lon,n,nh,nlat) = hmf2b(nh,nl+1);
else
  nmf2(lon,n,nh,nlat) = -1.;
  hmf2(lon,n,nh,nlat) = 0.; 
 if (hmf2b(nh,nl)>100. & hmf2b(nh,nl) <800.)
  nmf2(lon,n,nh,nlat) = nmf2b(nh,nl);
  hmf2(lon,n,nh,nlat) = hmf2b(nh,nl);
 end
 if (nl<71)
  if (hmf2b(nh,nl+2)>100. & hmf2b(nh,nl+2) <800.)
  nmf2(lon,n,nh,nlat) = nmf2b(nh,nl+2);
  hmf2(lon,n,nh,nlat) = hmf2b(nh,nl+2);
  end
 end
end		% if,else
end	% nl=1:2:72
end	% nh=1:24
else
% Average 2 lats together to go from 72 to 71 degrees for IGS comparisons
for nh=1:24
for nl=1:71
nlat = nl;
% Choose nl+1 latitude (center of 5deg bin in data)
tec(lon,n,nh,nlat) = 0.5*(tecb(nh,nl)+tecb(nh,nl+1));
wiav(lon,n,nh,nlat) = 0.5*(wiavb(nh,nl)+wiavb(nh,nl+1));
wiavl(lon,n,nh,nlat) = 0.5*(wiavlo(nh,nl)+wiavlo(nh,nl+1));
uiav(lon,n,nh,nlat) = 0.5*(uiavlo(nh,nl)+uiavlo(nh,nl+1));
if (hmf2b(nh,nl)>100. & hmf2b(nh,nl+1)>100. & hmf2b(nh,nl) <800. & hmf2b(nh,nl+1) <800.)
  nmf2(lon,n,nh,nlat) = 0.5*(nmf2b(nh,nl)+nmf2b(nh,nl+1));
  hmf2(lon,n,nh,nlat) = 0.5*(hmf2b(nh,nl)+hmf2b(nh,nl+1));
else
  nmf2(lon,n,nh,nlat) = -1.;
  hmf2(lon,n,nh,nlat) = 0.; 
 if (hmf2b(nh,nl)>100. & hmf2b(nh,nl) <800.)
  nmf2(lon,n,nh,nlat) = nmf2b(nh,nl);
  hmf2(lon,n,nh,nlat) = hmf2b(nh,nl);
 end
 if (hmf2b(nh,nl+1)>100. & hmf2b(nh,nl+1) <800.)
  nmf2(lon,n,nh,nlat) = nmf2b(nh,nl+1);
  hmf2(lon,n,nh,nlat) = hmf2b(nh,nl+1);
 end 
end		% if,else
end	% nl=1:71
end	% nh=1:24
end	% if5deg==1
end	% if n36==72
end	% n=1:ndays

% Calculate the median TEC over n days for each longitude
for n12=1:3
for nh=1:24
 for nl=1:nl36
  clear tecnd;
  clear nmf2nd;
  clear hmf2nd;
  clear wind;
  clear windlo;
  clear uie;
  num = 0;
  m = 0;
  for n=1:ndays
% revise 0 UT n=42 (1 Jan 2008) for dres
% revise 0 UT n=42 (1 Jan 2008) for dres
   if n36==72 & nh==1 & n==42
    tec(lon,n,nh,nl) = tec(lon,n+1,nh,nl);
    hmf2(lon,n,nh,nl) = 0.;
    wiav(lon,n,nh,nl) = wiav(lon,n+1,nh,nl);
    wiavl(lon,n,nh,nl) = wiavl(lon,n+1,nh,nl);
    uiav(lon,n,nh,nl) = uiav(lon,n+1,nh,nl);
   end
   if (n12==3)
    m = m + 1;
    tecnd(m) = tec(lon,n,nh,nl);
   if hmf2(lon,n,nh,nl)>100.
% 1_TIE-GCM bad hmnf2 at: lon=80, lat=+5; lon=140 lat=-5 to +10; lon=165 lat=-5 to +5 and 80-85N(13UT);
%   lon=190 lat=-10 to 0; lon=240 lat=-15 to -5
    num = num + 1;
    nmf2nd(num) = nmf2(lon,n,nh,nl);
    hmf2nd(num) = hmf2(lon,n,nh,nl);
   end 
   else	% if n12==3  
   for nn=1:25
    n25 = nd1hi(nn);
    if (n12==2)
     n25 = nd2lo(nn);
    end
  if (n==n25)
   m = m + 1;
   tecnd(m) = tec(lon,n,nh,nl);
   if hmf2(lon,n,nh,nl)>100.
% 1_TIE-GCM bad hmnf2 at: lon=80, lat=+5; lon=140 lat=-5 to +10; lon=165 lat=-5 to +5 and 80-85N(13UT);
%   lon=190 lat=-10 to 0; lon=240 lat=-15 to -5
    num = num + 1;
    nmf2nd(num) = nmf2(lon,n,nh,nl);
    hmf2nd(num) = hmf2(lon,n,nh,nl);
   end 
  end		% if n==n25
  end		% nn=1:25 
  end		% n12~=3  
  end		% n=1:ndays
   medtec(lon,nh,nl) = median(tecnd);
   mednmf2(lon,nh,nl) = median(nmf2nd);
   medhmf2(lon,nh,nl) = median(hmf2nd);
   avetec(lon,nh,nl) = mean(tecnd);
   avenmf2(lon,nh,nl) = mean(nmf2nd);
   avehmf2(lon,nh,nl) = mean(hmf2nd);
   stdtec = 0.;
   stdnmf2 = 0.;
   stdhmf2 = 0.;
   for n=1:m
    stdtec = stdtec + (avetec(lon,nh,nl)-tecnd(n))^2;
   end
   for n=1:num
    stdnmf2 = stdnmf2 + (avenmf2(lon,nh,nl)-nmf2nd(n))^2;
    stdhmf2 = stdhmf2 + (avehmf2(lon,nh,nl)-hmf2nd(n))^2;
   end
    stdtec = sqrt(stdtec/m);
    stdnmf2 = sqrt(stdnmf2/num);
    stdhmf2 = sqrt(stdhmf2/num);
if (n36==72)
  mw = 0;
  for n=1:ndays
   if (n12==3)
    mw = mw + 1;
    wind(mw) = wiav(lon,n,nh,nl);
    windlo(mw) = wiavl(lon,n,nh,nl);
    uie(mw) = uiav(lon,n,nh,nl);
   else	% if n12==3 
   for nn=1:25
    n25 = nd1hi(nn);
    if (n12==2)
     n25 = nd2lo(nn);
    end
  if (n==n25)
   mw = mw + 1;
   wind(mw) = wiav(lon,n,nh,nl);
    windlo(mw) = wiavl(lon,n,nh,nl);
    uie(mw) = uiav(lon,n,nh,nl);
  end		% n==n25
  end		% nn=1:25
  end		% n12~=3
  end		% n=1:ndays
   medwi(lon,nh,nl) = median(wind);
   avewi(lon,nh,nl) = mean(wind);   
   stdwi = 0.;
   medwil(lon,nh,nl) = median(windlo);
   avewil(lon,nh,nl) = mean(windlo);   
   stdwil = 0.;
   meduil(lon,nh,nl) = median(uie);
   aveuil(lon,nh,nl) = mean(uie);   
   stduil = 0.;
  for n=1:mw
   stdwi = stdwi + (avewi(lon,nh,nl)-wind(n))^2;
   stdwil = stdwil + (avewil(lon,nh,nl)-windlo(n))^2;
   stduil = stduil + (aveuil(lon,nh,nl)-uie(n))^2;
  end 
   stdwi = sqrt(stdwi/mw); 
   stdwil = sqrt(stdwil/mw); 
   stduil = sqrt(stduil/mw); 
end	% if n36=72

if (n12==1)
if (n36==36)
% Print out median tec over ndays
  fprintf(OUT1,'%2.2d %5.1f %3.3d %3.3d %5.2f %5.2f %5.2f %7.0f %7.0f %7.0f %5.1f %5.1f %5.1f %2.2d %2.2d\n',nh-1,beglat(nl),glon1,glon2,medtec(lon,nh,nl),avetec(lon,nh,nl),stdtec,mednmf2(lon,nh,nl),avenmf2(lon,nh,nl),stdnmf2,medhmf2(lon,nh,nl),avehmf2(lon,nh,nl),stdhmf2,ndays,num);
else
  fprintf(OUT1,'%2.2d %5.1f %3.3d %3.3d %5.2f %5.2f %5.2f %7.0f %7.0f %7.0f %5.1f %5.1f %5.1f %2.2d %2.2d %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f\n',nh-1,beglat(nl),glon1,glon2,medtec(lon,nh,nl),avetec(lon,nh,nl),stdtec,mednmf2(lon,nh,nl),avenmf2(lon,nh,nl),stdnmf2,medhmf2(lon,nh,nl),avehmf2(lon,nh,nl),stdhmf2,ndays,num,medwi(lon,nh,nl),avewi(lon,nh,nl),stdwi,medwil(lon,nh,nl),avewil(lon,nh,nl),stdwil,meduil(lon,nh,nl),aveuil(lon,nh,nl),stduil);
end
end	% if n12==1
if (n12==2)
if (n36==36)
% Print out median tec over ndays
  fprintf(OUT2,'%2.2d %5.1f %3.3d %3.3d %5.2f %5.2f %5.2f %7.0f %7.0f %7.0f %5.1f %5.1f %5.1f %2.2d %2.2d\n',nh-1,beglat(nl),glon1,glon2,medtec(lon,nh,nl),avetec(lon,nh,nl),stdtec,mednmf2(lon,nh,nl),avenmf2(lon,nh,nl),stdnmf2,medhmf2(lon,nh,nl),avehmf2(lon,nh,nl),stdhmf2,ndays,num);
else
  fprintf(OUT2,'%2.2d %5.1f %3.3d %3.3d %5.2f %5.2f %5.2f %7.0f %7.0f %7.0f %5.1f %5.1f %5.1f %2.2d %2.2d %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f\n',nh-1,beglat(nl),glon1,glon2,medtec(lon,nh,nl),avetec(lon,nh,nl),stdtec,mednmf2(lon,nh,nl),avenmf2(lon,nh,nl),stdnmf2,medhmf2(lon,nh,nl),avehmf2(lon,nh,nl),stdhmf2,ndays,num,medwi(lon,nh,nl),avewi(lon,nh,nl),stdwi,medwil(lon,nh,nl),avewil(lon,nh,nl),stdwil,meduil(lon,nh,nl),aveuil(lon,nh,nl),stduil);
end
end	% if n12==2
if (n12==3)
if (n36==36)
% Print out median tec over ndays
  fprintf(OUT,'%2.2d %2.2d %3.3d %3.3d %5.2f %5.2f %5.2f %7.0f %7.0f %7.0f %5.1f %5.1f %5.1f %2.2d %2.2d\n',nh-1,beglat(nl),glon1,glon2,medtec(lon,nh,nl),avetec(lon,nh,nl),stdtec,mednmf2(lon,nh,nl),avenmf2(lon,nh,nl),stdnmf2,medhmf2(lon,nh,nl),avehmf2(lon,nh,nl),stdhmf2,ndays,num);
else
  fprintf(OUT,'%2.2d %5.1f %3.3d %3.3d %5.2f %5.2f %5.2f %7.0f %7.0f %7.0f %5.1f %5.1f %5.1f %2.2d %2.2d %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f\n',nh-1,beglat(nl),glon1,glon2,medtec(lon,nh,nl),avetec(lon,nh,nl),stdtec,mednmf2(lon,nh,nl),avenmf2(lon,nh,nl),stdnmf2,medhmf2(lon,nh,nl),avehmf2(lon,nh,nl),stdhmf2,ndays,num,medwi(lon,nh,nl),avewi(lon,nh,nl),stdwi,medwil(lon,nh,nl),avewil(lon,nh,nl),stdwil,meduil(lon,nh,nl),aveuil(lon,nh,nl),stduil);
end
end	% if n12==3
 end	% nl=1:nl36
end	% nh=1:24   
end	% n12=1:2
   
end	% lon=1:8




