%%
%% 11/12 bae Put 8 glons together (zonal ave and also keep separate),
%%   average over LT for 61 days, and average +/-13.1deg (4 bins -10 
%%   to +10glat, then 5th bin is av of 2 bins from +/-(10-15glat) )
%% BA Emery, November 2012 /aim3/emery/cetic_climatology/cnofs8lon.m
clear all; close all




% if5deg = 1,0 for 5 deg or 2.5 deg dlat (only possible for IGS and TIE-WT)
if5deg = 1;

iftec = 0;
tecx = 35;
% ifmit=1,0 for mit or IGS TEC
ifmit = 1;
nt = 24;
if (iftec==1 & ifmit==0)
 nt = 12;
end
% if12 = 0 for 61-day climatology, =1,2 for 25 day climatology (hi,lo)
if12 = 0;
ifhmf2 = 0;

if iftec==1
 vartit = 'TEC'
 if (ifmit==1)
  gpstit='MIT'
 else
  gpstit='IGS'
 end
else
 gpstit='COSMIC'
 if ifhmf2==0
  vartit = 'NmF2'
 else
  vartit = 'hmF2'
 end
end


% Data
 loncutsE=[25,90,140,175,200,250,285,345,385]; % specify climatology longitudes to plot
 diflon=[40,65,50,35,25,50,35,60,40] 	% cannot do 360/8=45deg!
for j=1:9
 glon(j) = loncutsE(j)-diflon(j)/2.;
end
dlat = 5;		% 36*5deg=180 deg from 90S-90N	
nl1 = 17;
nl2 = 20;	
lat5deg=-90:dlat:90;
lat2p5deg = -87.5:2.5:90.;
	

min120 = 0:2:24;
min60=0:1:24;
dut = 1;

if (iftec==1 & ifmit==0)
 mingps=min120;
 dut = 2;
 if (if5deg==0)
  nl1 = 33;
  nl2 = nl1 + 9;
  dlat = 2.5;
  latgps=lat2p5deg;
 else
  latgps = lat5deg;
 end
else
 mingps=min60;	
 latgps = lat5deg;
end

nlat = length(latgps) - 1

tit = sprintf('%s %s +/-13.1glat %1.1dh LT vs glon',gpstit,vartit,dut)
if (iftec==0 & ifhmf2==0)
 tit = sprintf('COSMIC log10(NmF2) +/-13.1glat 1h LT vs glon');
end
outlow = sprintf('%s_%s_%1.1dh_LT_8glon_cnofs',gpstit,vartit,dut)
outlat = sprintf('%s_%s_%1.1dh_LT_8glon_glat',gpstit,vartit,dut)
outzlo = sprintf('%s_%s_%1.1dh_LT_zav_cnofs',gpstit,vartit,dut)
outzav = sprintf('%s_%s_%1.1dh_LT_zav_glat',gpstit,vartit,dut)
OUTLO = fopen(outlow,'w');
%OUTLA = fopen(outlat,'w');
OUTZL = fopen(outzlo,'w');
%OUTZA = fopen(outzav,'w');
fprintf(OUTZL,'%s %s %f deg dlat %f hour dut \n LT %s\n',gpstit,vartit,dlat,dut,gpstit);
fprintf(OUTLO,'%s %s %f deg dlat %f hour dut \n',gpstit,vartit,dlat,dut);
fprintf(OUTLO,'%d %d %d %d %d %d %d %d %d E glon \n LT 8glon %s\n',loncutsE,gpstit);


% Data
for it=1:24
 for nl=1:nlat
  for j=1:9
   nlt(it,nl,j) = 0;
   arrmed(it,nl,j) = 0;
  end
 end
end


for j=1:8

if iftec==1
if ifmit==1

 if if12==0
   DTECnam = sprintf('/hao/aim3/emery/cetic_climatology/mit_gps/long_%3.3d_07325_08020_tec_mit_all.dat',loncutsE(j))
 else
  if if12==1
   DTECnam = sprintf('/hao/aim3/emery/cetic_climatology/mit_gps/long_%3.3d_07325_08020_tec_mit_hi.dat',loncutsE(j))
  else
   DTECnam = sprintf('/hao/aim3/emery/cetic_climatology/mit_gps/long_%3.3d_07325_08020_tec_mit_lo.dat',loncutsE(j))
  end
 end
% ut_beg glat glon1 glon2 tecmed tecmean tecstd wttec averrtec numdays (max 61)
DTECin = fopen(DTECnam,'r');

for nh=1:nt	% UT=1:23 start of hourly bins
     slt = mingps(nh) + loncutsE(j)/15.;
     if (slt>=24.)
      slt = slt - 24.;
     end
     it = fix(slt/dut) + 1;	
 for nl=1:nlat
 linegps = fgetl(DTECin);
 ingps = sscanf(linegps,'%f');
  numtec = ingps(10);
  if (ingps(10) > 0)
   if (ingps(5) > 0.) 
    tecarrmed = ingps(5);
    nlt(it,nl,j) = nlt(it,nl,j) + 1;
    arrmed(it,nl,j) = arrmed(it,nl,j) + tecarrmed;
   else
    numtec = 0;
   end
  end

 end	% for 1:nlat
end	% for nh=1:nt

else
% Read IGS TEC every 2 hours (deleted every other glat to have 5 deg glat)

 if if12==0
  if (if5deg==1)
   DTECnam = sprintf('/hao/aim3/emery/cetic_climatology/igs_tec/long_%3.3d_07325_08020_tec_igs.dat',loncutsE(j))
  else
   DTECnam = sprintf('/hao/aim3/emery/cetic_climatology/igs_tec/long_%3.3d_07325_08020_tec_igs_2p5deg.dat',loncutsE(j))
  end
 else
  if if12==1
   if (if5deg==1)
   DTECnam = sprintf('/hao/aim3/emery/cetic_climatology/igs_tec/long_%3.3d_07325_08020_tec_igs_hi.dat',loncutsE(j))
   else
   DTECnam = sprintf('/hao/aim3/emery/cetic_climatology/igs_tec/long_%3.3d_07325_08020_tec_igs_hi_2p5deg.dat',loncutsE(j))
   end
  else
   if (if5deg==1)
   DTECnam = sprintf('/hao/aim3/emery/cetic_climatology/igs_tec/long_%3.3d_07325_08020_tec_igs_lo.dat',loncutsE(j))
   else
   DTECnam = sprintf('/hao/aim3/emery/cetic_climatology/igs_tec/long_%3.3d_07325_08020_tec_igs_lo_2p5deg.dat',loncutsE(j))
   end
  end
 end
% utmid glatmid glonmid tecmed tecmean tecstd numdays (max 61)
DTECin = fopen(DTECnam,'r');

for nh=1:nt	% UT=1:23 start of hourly bins
     slt = mingps(nh) + loncutsE(j)/15.;
     if (slt>=24.)
      slt = slt - 24.;
     end
     it = fix(slt/dut) + 1;	
 for nl=1:nlat
 linegps = fgetl(DTECin);
 ingps = sscanf(linegps,'%f');
  numtec = ingps(7);
    tecarrmed = ingps(4);
    nlt(it,nl,j) = nlt(it,nl,j) + 1;
    arrmed(it,nl,j) = arrmed(it,nl,j) + tecarrmed;
 end	% for 1:nlat
end	% for nh=1:nt

end		% end of iftec==1
else		
 if ifhmf2==0
DCOSMICnam = sprintf('/hao/aim3/emery/cetic_climatology/usu_cosmic/COSMIC_Medians_07325_08020/data_07325_08020/NmF2_%3.3d_%3.3d_07325_08020.dat',loncutsE(j)-10,loncutsE(j)+15)
 else
DCOSMICnam = sprintf('/hao/aim3/emery/cetic_climatology/usu_cosmic/COSMIC_Medians_07325_08020/data_07325_08020/hmF2_%3.3d_%3.3d_07325_08020.dat',loncutsE(j)-10,loncutsE(j)+15)
 end
DCOSMICin = fopen(DCOSMICnam,'r');

for nh=1:nt	% UT=1:23 start of hourly bins
     slt = mingps(nh) + loncutsE(j)/15.;
     if (slt>=24.)
      slt = slt - 24.;
     end
     it = fix(slt/dut) + 1;	
 for nl=1:nlat
 linecosmic = fgetl(DCOSMICin);
 incosmic = sscanf(linecosmic,'%f');
  numcosmic = incosmic(10);
  if (incosmic(10) > 0) 
   cosmicarrmed = incosmic(8);
    nlt(it,nl,j) = nlt(it,nl,j) + 1;
    arrmed(it,nl,j) = arrmed(it,nl,j) + cosmicarrmed;
  end
 end	% for 1:nlat
end	% for nh=1:nt

end


    end				% end of plotting 8 longitudes

% Divide sum by nlt to get average
for it=1:nt
 for j=1:8
 nltlow(it,j) = 0;
 arrlow(it,j) = 0.;
  for nl=1:nlat
   if (nlt(it,nl,j)>0)
    arrmed(it,nl,j) = arrmed(it,nl,j)/nlt(it,nl,j);
    if (nl>=nl1 & nl<=nl2)
     nltlow(it,j) = nltlow(it,j) + 1;
     arrlow(it,j) = arrlow(it,j) + arrmed(it,nl,j);
    end
   end
  end	% for nl=1:nlat
  if (if5deg==1)
   if (nlt(it,nl1-1,j)>0 | nlt(it,nl2+1,j)>0)
    facs = 0.5;
    facn = 0.5;
    if (nlt(it,nl1-1,j)>0 & nlt(it,nl2+1,j)==0)
     facs = 1.;
     facn = 0.;
    end
    if (nlt(it,nl1-1,j)==0 & nlt(it,nl2+1,j)>0)
     facs = 0.;
     facn = 1.;
    end
    nltlow(it,j) = nltlow(it,j) + 1;
    arrlow(it,j) = arrlow(it,j) + facs*arrmed(it,nl1-1,j) + facn*arrmed(it,nl2+1,j);
   end
  end	% if (if5deg==1)
% Divide sum by nltlow to get average
  if (nltlow(it,j)>0)
   arrlow(it,j) = arrlow(it,j)/nltlow(it,j);
  end
 end	% for j=1:8
% Print out low-lats for 8 glons
fprintf (OUTLO,'%5.1f %d %6.2f %d %6.2f %d %6.2f %d %6.2f %d %6.2f %d %6.2f %d %6.2f %d %6.2f\n',mingps(it)+dut/2.,nltlow(it,1),arrlow(it,1),nltlow(it,2),arrlow(it,2),nltlow(it,3),arrlow(it,3), ...
  nltlow(it,4),arrlow(it,4),nltlow(it,5),arrlow(it,5),nltlow(it,6),arrlow(it,6),nltlow(it,7),arrlow(it,7),nltlow(it,8),arrlow(it,8));
% Average over 8 glons
  nltzav(it) = 0;
  arrzav(it) = 0.;
  for j=1:8
   if (nltlow(it,j)>0)
    nltzav(it) = nltzav(it) + 1;
    arrzav(it) = arrzav(it) + arrlow(it,j);
   end
  end	% for j=1:8
  if (nltzav(it)>0)
   arrzav(it) = arrzav(it)/nltzav(it);
  end
% Print out low-lats for zonal ave
fprintf (OUTZL,'%5.1f %d %6.2f \n',mingps(it)+dut/2.,nltzav(it),arrzav(it));
% zonal ave for all nlat
 for nl=1:nlat
  allzav(it,nl) = 0.;
  numall(it,nl) = 0;
  for j=1:8
   if (nlt(it,nl,j)>0)
    numall(it,nl) = numall(it,nl) + 1;
    allzav(it,nl) = allzav(it,nl) + arrmed(it,nl,j);
   end
  end
% Divide sum by numall
 if (numall(it,nl)>0)
  allzav(it,nl) = allzav(it,nl)/numall(it,nl);
 else
  allzav(it,nl) = NaN;
 end
 end	% for nl=1:nlat
end	% for it=1:24

fsize = 12;

% Plot out low-lat LT vs glon
% Set 9th j to be same as first for plots
 j = 9;
 for it=1:nt
  arrlow(it,j) = arrlow(it,1);
 end
for j=1:9
 arrlow(nt+1,j) = arrlow(1,j);
end

figure
if (iftec==1)
            pcolor(mingps,glon,arrlow')
            caxis([0 tecx]);		% 
else
 if (ifhmf2==1)
            pcolor(mingps,glon,arrlow')
            caxis([200 400]);
 else
            pcolor(mingps,glon,log10(arrlow)')
            caxis([4.5 6.5]); 
 end
end
           title([tit],'FontSize',15,'FontWeight','bold');
            shading flat
            xlabel('LT (hours)','FontSize',fsize,'FontWeight','bold');
            ylabel('Longitude','FontSize',fsize,'FontWeight','bold');
	    set (gca,'xlim',[0,24],'Xtick',([0:6:24]),'XTickLabel',{'00:00','06:00','12:00','18:00','24:00'});
    hc=colorbar('horiz');
hold off;

% Plot zonal av vs glat
% Set nlat+1 to be same as first for plots
 for it=1:nt
  allzav(it,nlat+1) = allzav(it,1);
 end
for j=1:nlat+1
 allzav(nt+1,j) = allzav(1,j);
end
% Plot out
figure
if (iftec==1)
            pcolor(mingps,latgps,allzav')
            caxis([0 35]);		% for mit and igs
else
 if (ifhmf2==1)
            pcolor(mingps,latgps,allzav')
            caxis([200 400]);
 else
            pcolor(mingps,latgps,log10(allzav)')
            caxis([3.5 6.5]); 
 end
end
           title([tit],'FontSize',15,'FontWeight','bold');
            shading flat
            xlabel('LT (hours)','FontSize',fsize,'FontWeight','bold');
            ylabel('Latitude','FontSize',fsize,'FontWeight','bold');
	    set (gca,'xlim',[0,24],'Xtick',([0:6:24]),'XTickLabel',{'00:00','06:00','12:00','18:00','24:00'});
    hc=colorbar('horiz');
hold off;
