clear all;

% May 2012  Program to calculate the median of values of 400km density from Luehr's
%    CHAMP into 3-deg lat bins and 1-h LT bins for various days
%    for the Emery et al., June 2012 MSU talk and CEDAR Poster

% 332-354km for Luehr instead of 337-368km for France/U CO/Sutton

fsize = 12;

%ifall = 1,0 if get all 61 days or hi/lo days
ifall = 1
% Density in /hao/aim3/stans/data/champ/2007(8).acc at 337-368km,400 km in 6,7 cols (MSIS 8,9 cols)
% Copied over to /hao/aim3/emery/cetic_climatology/champ_density as 07325_08020_luehr
% Date in Julian fractional days (2881-2941, missing from Jan 5-9 2926.5-2930.5)
% ~8640 lines/day (skip 2617919 lines in 2007.acc and read 354240 lines; read 138234 lines in 2008.acc)
% Found 07356 (Dec 22) missing also! 
n13 = 6;
%n13 = 7;



if (ifall==1)
% Name 'global' den 400km output file for daily global means (from equal wt ascent/descent)
 outglob=sprintf('champ_density/global_sat_400km_LTden_luehr_07325_08020_b.dat')
OUTG = fopen(outglob,'w');
end	% ifall==1

if (n13>6)
 kmdantit = sprintf('at %dkm 07325-08020',alt);
else
 kmdantit = sprintf('~337-368km 07325-08020');
end
hitit = sprintf('Kp>=2 Vsw>=500km/s');
lotit = sprintf('Kp<=1 Vsw<=450km/s');

nt1 = 1;
nt2 = 2;
if (ifall==1)
 nt1 = 3;
 nt2 = 3;
end

% CHAMP densities missing for days 08_006-8

lt24 = 0:1:23
%lat = -87:3:87
lat = -90:3:87

% 25 daynumbers for Kp>=2 and Vsw>=500km/s: 07325-330, 345-6, 351-7, 08005-8, 14-19
ndahi = [325,326,327,328,329,330,345,346,351,352,353,354,355,356,357, 5,6,7,8,14,15,16,17,18,19]
% 25 daynumbers for Kp=<1 and Vsw=<450km/s: 07332-343, 349-350, 358-360, 362-4, 08001-4, 11
ndalo = [332,333,334,335,336,337,338,339,340,341,342,343,349,350,358,359,360,362,363,364, 1,2,3,4,11]


% Zero out LT,glat arrays for accumulations
for i=1:24
 for j=1:60
  for k=1:3
   num(i,j,k) = 0;
   den(i,j,k) = 0.;
   denmd(i,j,k) = 0.;
   den2(i,j,k) = 0.;
   denmd2(i,j,k) = 0.;
   ht(i,j,k) = 0.;
  end
 end
end

neqlt = 24;
xeqlt = 0;

for nt=nt1:nt2
nht = 5000.;
xht = 0.;
if (nt<3)
 ndt = 25;
else
 ndt = 61;
%if (n13==12)
%   ndt = 61+50;
% end
end
% Luehr file
 denfile=sprintf('/hao/aim3/emery/cetic_climatology/champ_density/07325_08020_luehr');
 datden = dlmread(denfile,'',0,0);
 nll=length(datden);

for k=1:2
 numad(k) = 0;
 dden(k) = 0.;
 ddenmd(k) = 0.;
 dden2(k) = 0.;
 ddenmd2(k) = 0.;
end
ndap = 325;
for n=1:nll
 nda = fix(datden(n,1)) - 2881 + 325;
% Calculate global ascent/descent values at end of each day
 if nda>ndap
  if (ndap>365)
   iyr = 2008;
   jda = ndap-365;
  else
   iyr = 2007;
   jda = ndap;
  end
  fprintf(OUTG,'%4.4d %3.3d %4.4d %12.4e %12.4e %12.4e %12.4e %4.4d %12.4e %12.4e %12.4e %12.4e \n',iyr,jda,numad(1),dden(1)/numad(1),ddenmd(1)/numad(1),dden2(1)/numad(1),ddenmd2(1)/numad(1),numad(2),dden(2)/numad(2),ddenmd(2)/numad(2),dden2(2)/numad(2),ddenmd2(2)/numad(2));
for k=1:2
 numad(k) = 0;
 dden(k) = 0.;
 ddenmd(k) = 0.;
 dden2(k) = 0.;
 ddenmd2(k) = 0.;
end
 end	% nda>ndap
% nad = 1 for ascent (8.3 to 2.8LT) and 2 for descent (20.4 to 14.8LT) 
 nad = 2;
if (datden(n,2)<12.)
 nad = 1;
end
 nlt = fix(datden(n,2)) + 1;
 nlat = fix((datden(n,4)+90)/3.) + 1;

% Accumulate in arrays the density values at 400km or at 345km
 num(nlt,nlat,nt) = num(nlt,nlat,nt) + 1;
 den(nlt,nlat,nt) = den(nlt,nlat,nt) + datden(n,6);
 den2(nlt,nlat,nt) = den2(nlt,nlat,nt) + datden(n,7);
 denmd(nlt,nlat,nt) = denmd(nlt,nlat,nt) + datden(n,8);
 denmd2(nlt,nlat,nt) = denmd2(nlt,nlat,nt) + datden(n,9);
 ht(nlt,nlat,nt) = ht(nlt,nlat,nt) + datden(n,3);
 numad(nad) = numad(nad) + 1;
 dden(nad) = dden(nad) + datden(n,6);
 dden2(nad) = dden2(nad) + datden(n,7);
 ddenmd(nad) = ddenmd(nad) + datden(n,8);
 ddenmd2(nad) = ddenmd2(nad) + datden(n,9);
 if (datden(n,3)<nht)
  nht = datden(n,3);
 end
 if (datden(n,3)>xht)
  xht = datden(n,3);
 end
 ndap = nda;
end	% n=1:nl
nht,xht

% Calculate global ascent/descent values at end of each day
   iyr = 2008;
   jda = ndap-365;
  fprintf(OUTG,'%4.4d %3.3d %4.4d %12.4e %12.4e %12.4e %12.4e %4.4d %12.4e %12.4e %12.4e %12.4e \n',iyr,jda,numad(1),dden(1)/numad(1),ddenmd(1)/numad(1),dden2(1)/numad(1),ddenmd2(1)/numad(1),numad(2),dden(2)/numad(2),ddenmd(2)/numad(2),dden2(2)/numad(2),ddenmd2(2)/numad(2));

end	% nt=nt1:nt2


% Print out 2 files of den over 25 days 
for nt=nt1:nt2
  outfile = sprintf('champ_density/all_den_lat_vs_lt_luehr')
 if nt==1
  outfile = sprintf('champ_density/hi_den_%dkm_lat_vs_lt_luehr',alt)
 end
 if nt==2
  outfile = sprintf('champ_density/lo_den_%dkm_lat_vs_lt_luehr',alt)
 end
OUT = fopen(outfile,'w');
for nh=1:24
 for nl=1:60
  if (num(nh,nl,nt)>0)
    den(nh,nl,nt) = den(nh,nl,nt)/num(nh,nl,nt);
    denmd(nh,nl,nt) = denmd(nh,nl,nt)/num(nh,nl,nt);
    den2(nh,nl,nt) = den2(nh,nl,nt)/num(nh,nl,nt);
    denmd2(nh,nl,nt) = denmd2(nh,nl,nt)/num(nh,nl,nt);
  end
  if (nl>1)
  fprintf(OUT,'%5.1f %2.2d %d %12.4e %12.4e %12.4e %12.4e \n', ...
   lt24(nh)+0.5,lat(nl),num(nh,nl,nt),den(nh,nl,nt),denmd(nh,nl,nt), ...
   den2(nh,nl,nt),denmd2(nh,nl,nt));
  end
% Set 0's to NaN for plots
  if (num(nh,nl,nt)==0)
   den(nh,nl,nt) = NaN;
   denmd(nh,nl,nt) = NaN;
   den2(nh,nl,nt) = NaN;
   denmd2(nh,nl,nt) = NaN;
  else
   den(nh,nl,nt) = den(nh,nl,nt)*1.e+13;
   denmd(nh,nl,nt) = denmd(nh,nl,nt)*1.e+13;
   den2(nh,nl,nt) = den2(nh,nl,nt)*1.e+13;
   denmd2(nh,nl,nt) = denmd2(nh,nl,nt)*1.e+13;
  end
 end	% nl=1:60
end	% nh=1:24   

end	% nt=nt1:nt2

if (nt2>nt1)
tit1num = [hitit,' CHAMP number in bins ',kmdantit];
tit1ave = [hitit,' CHAMP average density*1.e+13 in kg/m3 ',kmdantit];
tit1med = [hitit,' CHAMP median density*1.e+13 in kg/m3 ',kmdantit];
tit2num = [lotit,' CHAMP number in bins ',kmdantit];
tit2ave = [lotit,' CHAMP average density*1.e+13 in kg/m3 ',kmdantit];
tit2med = [lotit,' CHAMP median density*1.e+13 in kg/m3 ',kmdantit];


figure
% position [left bottom width height]
% Set size to be 8.5x11 inches
%set(gcf,'position',[142    99   647   822])
%H = gcf;
clf

subplot(2,1,1)
            pcolor(lt24,lat,den(:,:,1)')
            title([tit1ave]);
            caxis([nden xden]);
            ylabel('Latitude');
	    set (gca,'xlim',[0,24],'Xtick',([0:6:24]),'XTickLabel',{'00:00','06:00','12:00','18:00','24:00'});

    hc=colorbar('horiz');

subplot(2,1,2)
            pcolor(lt24,lat,den(:,:,2)')
            title([tit2ave]);
            caxis([nden xden]);
            shading flat
            xlabel('LT');
            ylabel('Latitude');
	    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;	


figure
% position [left bottom width height]
% Set size to be 8.5x11 inches
%set(gcf,'position',[142    99   647   822])
%H = gcf;
clf

subplot(2,1,1)
            pcolor(lt24,lat,denmd(:,:,1)')
            title([tit1med]);
            caxis([nden xden]);
            shading flat
            xlabel('LT');
            ylabel('Latitude');
	    set (gca,'xlim',[0,24],'Xtick',([0:6:24]),'XTickLabel',{'00:00','06:00','12:00','18:00','24:00'});

    hc=colorbar('horiz');

subplot(2,1,2)
            pcolor(lt24,lat,denmd(:,:,2)')
            title([tit2med]);
            caxis([nden xden]);
            shading flat
            xlabel('LT');
            ylabel('Latitude');
	    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;	

figure
% position [left bottom width height]
% Set size to be 8.5x11 inches
%set(gcf,'position',[142    99   647   822])
%H = gcf;
clf

subplot(2,1,1)
            pcolor(lt24,lat,num(:,:,1)')
            title([tit1num]);
            caxis([0 150]);
            shading flat
            xlabel('LT');
            ylabel('Latitude');
	    set (gca,'xlim',[0,24],'Xtick',([0:6:24]),'XTickLabel',{'00:00','06:00','12:00','18:00','24:00'});

    hc=colorbar('horiz');

subplot(2,1,2)
            pcolor(lt24,lat,num(:,:,2)')
            title([tit2num]);
            caxis([0 150]);
            shading flat
            xlabel('LT');
            ylabel('Latitude');
	    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;	

else

% Plot All CHAMP densities
figure
% position [left bottom width height]
% Set size to be 8.5x11 inches
%set(gcf,'position',[142    99   647   822])
%H = gcf;
clf

alt = 345;
nden = 7.5;
xden = 50;

titave = ['All CHAMP av density*1.e+13 in kg/m3 at ~332-354km (Luehr) 07325-08020'];
titavmsis = ['All MSIS00 av density*1.e+13 in kg/m3 at ~332-354km (Luehr) 07325-08020'];
subplot(2,1,1)

            pcolor(lt24,lat,den(:,:,3)')
            title([titave]);
            caxis([nden xden]);
            shading flat
            xlabel('LT');
            ylabel('Latitude');
	    set (gca,'xlim',[0,24],'Xtick',([0:6:24]),'XTickLabel',{'00:00','06:00','12:00','18:00','24:00'});
    hc=colorbar('horiz');

subplot(2,1,2)
            pcolor(lt24,lat,denmd(:,:,3)')
            title([titavmsis]);
            caxis([nden xden]);
            shading flat
            xlabel('LT');
            ylabel('Latitude');
	    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;	

figure
% position [left bottom width height]
% Set size to be 8.5x11 inches
%set(gcf,'position',[142    99   647   822])
%H = gcf;
clf
%subplot('position',[0.10,0.10,0.87,0.85])
hold on;
            pcolor(lt24,lat,den(:,:,3)')
            title([titave],'FontSize',fsize,'Fontweight','bold');
            caxis([nden xden]);
            shading flat
	    grid;
            xlabel('LT','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'},'ylim',[-90,90]);
    hc=colorbar('horiz');

    hold off;	

figure
hold on;
            pcolor(lt24,lat,denmd(:,:,3)')
            title([titavmsis],'FontSize',fsize,'Fontweight','bold');
            caxis([nden xden]);
            shading flat
	    grid;
            xlabel('LT','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'},'ylim',[-90,90]);
    hc=colorbar('horiz');

    hold off;	

% Plot All 400km  densities
figure
% position [left bottom width height]
% Set size to be 8.5x11 inches
%set(gcf,'position',[142    99   647   822])
%H = gcf;
clf

 alt = 400;
 nden = 0;
 xden = 15;

titave = ['All CHAMP av density*1.e+13 in kg/m3 at 400km (Luehr) 07325-08020'];
titavmsis = ['All MSIS00 av density*1.e+13 in kg/m3 at 400km (Luehr) 07325-08020'];
subplot(2,1,1)

            pcolor(lt24,lat,den2(:,:,3)')
            title([titave]);
            caxis([nden xden]);
            shading flat
            xlabel('LT');
            ylabel('Latitude');
	    set (gca,'xlim',[0,24],'Xtick',([0:6:24]),'XTickLabel',{'00:00','06:00','12:00','18:00','24:00'});
    hc=colorbar('horiz');

subplot(2,1,2)
            pcolor(lt24,lat,denmd2(:,:,3)')
            title([titavmsis]);
            caxis([nden xden]);
            shading flat
            xlabel('LT');
            ylabel('Latitude');
	    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;	

figure
% position [left bottom width height]
% Set size to be 8.5x11 inches
%set(gcf,'position',[142    99   647   822])
%H = gcf;
clf
%subplot('position',[0.10,0.10,0.87,0.85])
hold on;
            pcolor(lt24,lat,den2(:,:,3)')
            title([titave],'FontSize',fsize,'Fontweight','bold');
            caxis([nden xden]);
            shading flat
	    grid;
            xlabel('LT','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'},'ylim',[-90,90]);
    hc=colorbar('horiz');

    hold off;	


figure
% position [left bottom width height]
% Set size to be 8.5x11 inches
%set(gcf,'position',[142    99   647   822])
%H = gcf;
clf
%subplot('position',[0.10,0.10,0.87,0.85])
hold on;
            pcolor(lt24,lat,denmd2(:,:,3)')
            title([titavmsis],'FontSize',fsize,'Fontweight','bold');
            caxis([nden xden]);
            shading flat
	    grid;
            xlabel('LT','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'},'ylim',[-90,90]);
    hc=colorbar('horiz');

    hold off;	

	
end	% for nt=3:3 or all 61-3=58 days
