clear all;
close all;

% Dec 2012  Program to read the 2 LT +/-13.1glat older legacy CHAMP winds from 
%  Kathrin Haeusler found in kh_winds_latvlt.m to plot as glon vs LT in 131d Dec solstice
% (function of daily Kp?)


fsize = 12;
unx=300.;
unn=-300.;
dnx=2.0e-12;

ndt = 131;

iyear0 = 2007;
%iyear0 = 2008;
%ndt = 75;
%iyear0 = 2006;
%ndt = 131;
%iyear0 = 2005;
%ndt = 131;

for ih=1:24
 for jl=1:24
  numt(jl,ih) = 0;
  unt(jl,ih) = 0.;
  d400t(jl,ih) = 0.;
 end
end

 out4y=sprintf('champ_density/lowlat_LTvsglon_haeusler_4y.dat')
 OUT = fopen(out4y,'w');

for iy=2005:2008

iyear0 = iy


if (iyear0==2005)
 outglob=sprintf('champ_density/global_LTden_haeusler_05_06.dat')
 wrap = sprintf('champ_density/daily_glon_2LTs_05_06')
end
if (iyear0==2006)
 outglob=sprintf('champ_density/global_LTden_haeusler_06_07.dat')
% Save +/-13.1glat w SLT av every day ht,Un,Den400 (av Kp,Ap,Vsw) plus SLT vs glon
 wrap = sprintf('champ_density/daily_glon_2LTs_06_07')
end
if (iyear0==2007)
% 2846-2976 MJD (290-420) vs 2881-2941 (325-385)
 outglob=sprintf('champ_density/global_LTden_haeusler_07_08.dat2')
% Save +/-13.1glat w SLT av every day ht,Un,Den400 (av Kp,Ap,Vsw) plus SLT vs glon
 wrap = sprintf('champ_density/daily_glon_2LTs_07_08')
end
if (iyear0==2008)
 ndt = 75;
 outglob=sprintf('champ_density/global_LTden_haeusler_08.dat')
 wrap = sprintf('champ_density/daily_glon_2LTs_08')
end
 ING = fopen(wrap,'r');
 linhed = fgetl(ING);
%  fprintf(OUTAP,'iyr,jda,f107,kpav,apav,vswav,numd(1),stld(1),htd(1),numd(2),stld(2),htd(2) then 4 lines Un1,2 and Den4001,2 for glon=0:15:355 bins \n');
datday = dlmread(outglob,'',0,0);
size(datday)


for ih=1:24
 for jl=1:24
  num(jl,ih) = 0;
  un(jl,ih) = 0.;
  d400(jl,ih) = 0.;
 end
end
for nd=1:ndt
 if (datday(nd,3)>0)
% Read 5 lines of day, Un(LT1,2) and d400(LT1,2)
  lineh = fgetl(ING);
  datah = sscanf(lineh,'%f');
  ih1 = fix(datah(8)) + 1;
  ih2 = fix(datah(11)) + 1;
  linu1 = fgetl(ING);
  datu1 = sscanf(linu1,'%f');
  linu2 = fgetl(ING);
  datu2 = sscanf(linu2,'%f');
  lind1 = fgetl(ING);
  datd1 = sscanf(lind1,'%f');
  lind2 = fgetl(ING);
  datd2 = sscanf(lind2,'%f');
  for jl=1:24
   if (isnan(datu1(jl))==0)
    num(jl,ih1) = num(jl,ih1) + 1;
    un(jl,ih1) = un(jl,ih1) + datu1(jl);
    d400(jl,ih1) = d400(jl,ih1) + datd1(jl);
    numt(jl,ih1) = numt(jl,ih1) + 1;
    unt(jl,ih1) = unt(jl,ih1) + datu1(jl);
    d400t(jl,ih1) = d400t(jl,ih1) + datd1(jl);
   end
   if (isnan(datu2(jl))==0)
    num(jl,ih2) = num(jl,ih2) + 1;
    un(jl,ih2) = un(jl,ih2) + datu2(jl);
    d400(jl,ih2) = d400(jl,ih2) + datd2(jl);
    numt(jl,ih2) = numt(jl,ih2) + 1;
    unt(jl,ih2) = unt(jl,ih2) + datu2(jl);
    d400t(jl,ih2) = d400t(jl,ih2) + datd2(jl);
   end
  end	%   for jl=1:24  
 end	% if (datday(nd,3)>0)
end	% for nd=1:ndt
for ih=1:24
 for jl=1:24
  if num(jl,ih)>0;
   un(jl,ih) = un(jl,ih)/num(jl,ih);
   d400(jl,ih) = d400(jl,ih)/num(jl,ih);
  else
   un(jl,ih) = NaN;
   d400(jl,ih) = NaN;
  end
 end
end

% Midpoints 24h and 60 lats
lt24 = 0.5:1:23.5;
lon=[7.5:15:352.5];
nlon = length(lon);

% Plot Low-lat CHAMP densities and winds as glon vs LT
titu = sprintf('%d Dec Low-lat CHAMP Haeusler Un (m/s)',iyear0);

figure
hold on;
            pcolor(lon,lt24,un')
            title([titu],'FontSize',fsize,'Fontweight','bold');
            caxis([unn unx]);
            shading flat
	    grid;
            xlabel('Longitude','FontSize',fsize,'Fontweight','bold');
            ylabel('LT','FontSize',fsize,'Fontweight','bold');
	    set (gca,'ylim',[0,24],'Ytick',([0:6:24]),'YTickLabel',{'00:00','06:00','12:00','18:00','24:00'},'xlim',[0,360]);
    hc=colorbar('horiz');

    hold off;	

titd = sprintf('%d Dec Low-lat CHAMP Haeusler Neutral Density 400km',iyear0);

figure
hold on;
            pcolor(lon,lt24,d400')
            title([titd],'FontSize',fsize,'Fontweight','bold');
            caxis([0 dnx]);
            shading flat
	    grid;
            xlabel('Longitude','FontSize',fsize,'Fontweight','bold');
            ylabel('LT','FontSize',fsize,'Fontweight','bold');
	    set (gca,'ylim',[0,24],'Ytick',([0:6:24]),'YTickLabel',{'00:00','06:00','12:00','18:00','24:00'},'xlim',[0,360]);
    hc=colorbar('horiz');

    hold off;	

end	% for iy=2005:2008


  fprintf(OUT,'Un vs LT for all glon and for glon=7.5:15:352.5 Dec solstice CHAMP 131d for 2005-2008 \n LT all_glon_Un(even wt; wted w #) Un(7.5E to 352.5E) \n');

for ih=1:24
 numlt(ih) = 0;
 unlt(ih) = 0.;
 uneven(ih) = 0.;
 for jl=1:24
  if numt(jl,ih)>0;
   unt(jl,ih) = unt(jl,ih)/numt(jl,ih);
   d400t(jl,ih) = d400t(jl,ih)/numt(jl,ih);
  else
   unt(jl,ih) = NaN;
   d400t(jl,ih) = NaN;
  end
  numlt(ih) = numlt(ih) + numt(jl,ih);
  unlt(ih) = unlt(ih) + numt(jl,ih)*unt(jl,ih);
  uneven(ih) = uneven(ih) + unt(jl,ih);
 end
 unlt(ih) = unlt(ih)/numlt(ih);
 uneven(ih) = uneven(ih)/24.;
  fprintf(OUT,'%5.1f %8.2f %8.2f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f \n',lt24(ih),uneven(ih),unlt(ih),unt(:,ih));
end

% Plot Low-lat CHAMP densities and winds as glon vs LT
titu = sprintf('2005-2008 Dec Low-lat CHAMP Haeusler Un (m/s)');

figure
hold on;
            pcolor(lon,lt24,unt')
            title([titu],'FontSize',fsize,'Fontweight','bold');
            caxis([unn unx]);
            shading flat
	    grid;
            xlabel('Longitude','FontSize',fsize,'Fontweight','bold');
            ylabel('LT','FontSize',fsize,'Fontweight','bold');
	    set (gca,'ylim',[0,24],'Ytick',([0:6:24]),'YTickLabel',{'00:00','06:00','12:00','18:00','24:00'},'xlim',[0,360]);
    hc=colorbar('horiz');

    hold off;	

titd = sprintf('2005-2008 Dec Low-lat CHAMP Haeusler Neutral Density 400km');

figure
hold on;
            pcolor(lon,lt24,d400t')
            title([titd],'FontSize',fsize,'Fontweight','bold');
            caxis([0 dnx]);
            shading flat
	    grid;
            xlabel('Longitude','FontSize',fsize,'Fontweight','bold');
            ylabel('LT','FontSize',fsize,'Fontweight','bold');
	    set (gca,'ylim',[0,24],'Ytick',([0:6:24]),'YTickLabel',{'00:00','06:00','12:00','18:00','24:00'},'xlim',[0,360]);
    hc=colorbar('horiz');

    hold off;	

zero = [0.,0.];
xtot = [0,24];

figure
hold on;
title('Four year 2005-8 131-d December solstice CHAMP Zonal Wind','Fontweight','bold');
plot (xtot,zero,'--k');
plot(lt24,uneven,'*-r','LineWidth',2);
plot(lt24,unlt,'--b');
text(10.,0.9*unx,'even wts','Color','r','FontSize',fsize,'Fontweight','bold');
text(10.,0.8*unx,'wted w # per glon','Color','b','FontSize',fsize,'Fontweight','bold');
grid;
            xlabel('LT','FontSize',fsize,'Fontweight','bold');
            ylabel('Un (m/s)','FontSize',fsize,'Fontweight','bold');
	    set (gca,'xlim',[0,24],'Xtick',([0:6:24]),'XTickLabel',{'00:00','06:00','12:00','18:00','24:00'},'ylim',[unn,unx]);
hold off;
