clear all;
close all;

% Nov 2012  Program to calculate the median of values of sat alt (~350km) and 
%  400km density from 28 Sep 2012 Doornbos official (Luehr) CHAMP values
%    into 3-deg lat bins and 1-h LT bins for various days
%    for the Emery et al., CCMC talk and Poster (w values +/-13.1 deg for C/NOFS)

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

fsize = 12;


% /hao/aim3/emery/cetic_climatology/champ_density accel.... files
%  400km den from MSIS as:  den(400) = den(z)*MSIS(400)/MSIS(z)
% Missing from Jan 6 >=0322UT and Jan 8 <1200UT) (missing when year=#008!)
% Missing also Dec 22 at 12:39:20 and 9 Jan 13:55:00-40
% 8640 lines/day with 10 sec data and missing with #007 or #008 for year
%  Had to replace #007 or #008 with 9999 to get '%d' OK for smaller gaps


% Name 'global' den 400km output file for daily global means (from equal wt ascent/descent)
 outglob=sprintf('champ_density/global_LTden_doornbos_07325_08020.datb')
OUTG = fopen(outglob,'w');


% CHAMP densities missing for days 2008 6-8

% Midpoints 24h and 60 lats
lt24 = 0.5:1:23.5;
lat = -88.5:3:88.5;

% 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,9999]
ndahi = [325,326,327,328,329,330,345,346,351,352,353,354,355,356,357, 370,371,372,373,379,380,381,382,383,384,9999]
% 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,9999]
ndalo = [332,333,334,335,336,337,338,339,340,341,342,343,349,350,358,359,360,362,363,364, 366,367,368,369,376,9999]


% Zero out LT,glat arrays for accumulations
for i=1:24
 for j=1:60
  for k=1:3
% 2=400km, md=MSIS, j8=Jacchia-Bowman 2008, hd=HASDM AF model, un,vn,wn(hwm)=winds(HWM07 model)
   num(i,j,k) = 0;
   den(i,j,k) = 0.;
   denmd(i,j,k) = 0.;
   denj8(i,j,k) = 0.;
   denhd(i,j,k) = 0.;
   un(i,j,k) = 0.;
   vn(i,j,k) = 0.;
   wn(i,j,k) = 0.;
   unhwm(i,j,k) = 0.;
   vnhwm(i,j,k) = 0.;
   wnhwm(i,j,k) = 0.;
   den2(i,j,k) = 0.;
   denmd2(i,j,k) = 0.;
   denj82(i,j,k) = 0.;
   denhd2(i,j,k) = 0.;
   ht(i,j,k) = 0.;
  end
 end
end

neqlt = 24;
xeqlt = 0;

nt = 3;
nhi = 1;
nlo = 1;

nht = 5.e+7;
xht = 0.;
if (nt<3)
 ndt = 25;
else
 ndt = 61;
end

for nda=325:385
mt = 0;
if (nda==ndahi(nhi))
 mt = 1;
end
if (nda==ndalo(nlo))
 mt = 2;
end

for k=1:2
% numad = number of ascent (8-3LT) k=1 and descent (20-15LT) k=2 times per day (d..=ad)
 numad(k) = 0;
 dden(k) = 0.;
 ddenmd(k) = 0.;
 ddenj8(k) = 0.;
 ddenhd(k) = 0.;
 dun(k) = 0.;
 dvn(k) = 0.;
 dwn(k) = 0.;
 dunhwm(k) = 0.;
 dvnhwm(k) = 0.;
 dwnhwm(k) = 0.;
 dden2(k) = 0.;
 ddenmd2(k) = 0.;
 ddenj82(k) = 0.;
 ddenhd2(k) = 0.;
end

% Doornbos files
 if (nda==325)
  infil=sprintf('/hao/aim3/emery/cetic_climatology/champ_density/acceldrag_datarequest_2007-11-21_31d.txt')
  INP=fopen(infil,'r');
  for n=1:21
   linhd = fgetl(INP);
  end
 end
 if (nda==356)
  infil=sprintf('/hao/aim3/emery/cetic_climatology/champ_density/acceldrag_datarequest_2007-12-22_10d.txt');
  INP=fopen(infil,'r');
  for n=1:21
   linhd = fgetl(INP);
  end
 end
 if (nda==366)
  infil=sprintf('/hao/aim3/emery/cetic_climatology/champ_density/acceldrag_datarequest_2008-01-01_20d.txt');
  INP=fopen(infil,'r');
  for n=1:21
   linhd = fgetl(INP);
  end
 end
 jda = nda
  if (nda>365)
   iyr = 2008;
   jda = nda-365;
  else
   iyr = 2007;
   jda = nda;
  end

for nut=1:24
for n10s=1:360
   line = fgetl(INP);
% Skip big gaps of # lines or read 9999 for date instead of #007 or #008 when missing
ifskp = 0;
if (nda>370 & nda<374)
 numskp = (nut-1)*360 + n10s;
 if (nda==371 & numskp>=1212)
  ifskp = 1;
 end
 if (nda==372)
  ifskp = 1;
 end
 if (nda==373 & numskp<4321)
  ifskp = 1;
 end
end
if (ifskp==0)
   ymdh = sscanf(line,'%d');
   year = ymdh(1);
   month = -ymdh(2);
   day = -ymdh(3);
   utsec = ymdh(4);
  if (year==2007 | year==2008) 
   input = sscanf(line,'%*27c %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e');
% 18 inputs: alt(m),glon,glat,slt,msis400,msisalt,jb08alt,hasdmalt,denalt,un,vn,wn,scyd,scyi,scyhwm07,unhwm,vnhwm,wnhwm
 nlt = fix(input(4)) + 1;
% nad = 1 for ascent (8.3 to 2.8LT) and 2 for descent (20.4 to 14.8LT) 
 nad = 2;
if (input(4)<12.)
 nad = 1;
end
 nlat = fix((input(3)+90)/3.) + 1;

% Accumulate in arrays the density values at 400km (2) or at ~348km
 num(nlt,nlat,nt) = num(nlt,nlat,nt) + 1;
 un(nlt,nlat,nt) = un(nlt,nlat,nt) + input(10);
 vn(nlt,nlat,nt) = vn(nlt,nlat,nt) + input(11);
 wn(nlt,nlat,nt) = wn(nlt,nlat,nt) + input(12);
 unhwm(nlt,nlat,nt) = unhwm(nlt,nlat,nt) + input(16);
 vnhwm(nlt,nlat,nt) = vnhwm(nlt,nlat,nt) + input(17);
 wnhwm(nlt,nlat,nt) = wnhwm(nlt,nlat,nt) + input(18);
 den(nlt,nlat,nt) = den(nlt,nlat,nt) + input(9);
 denmd(nlt,nlat,nt) = denmd(nlt,nlat,nt) + input(6);
 denj8(nlt,nlat,nt) = denj8(nlt,nlat,nt) + input(7);
 denhd(nlt,nlat,nt) = denhd(nlt,nlat,nt) + input(8);
 denmd2(nlt,nlat,nt) = denmd2(nlt,nlat,nt) + input(5);
 den2(nlt,nlat,nt) = den2(nlt,nlat,nt) + input(9)*input(5)/input(6);
 denj82(nlt,nlat,nt) = denj82(nlt,nlat,nt) + input(7)*input(5)/input(6);
 denhd2(nlt,nlat,nt) = denhd2(nlt,nlat,nt) + input(8)*input(5)/input(6);
 ht(nlt,nlat,nt) = ht(nlt,nlat,nt) + input(1);
 if (mt>0)
  ntsv = nt;
  nt = mt;
 num(nlt,nlat,nt) = num(nlt,nlat,nt) + 1;
 un(nlt,nlat,nt) = un(nlt,nlat,nt) + input(10);
 vn(nlt,nlat,nt) = vn(nlt,nlat,nt) + input(11);
 wn(nlt,nlat,nt) = wn(nlt,nlat,nt) + input(12);
 unhwm(nlt,nlat,nt) = unhwm(nlt,nlat,nt) + input(16);
 vnhwm(nlt,nlat,nt) = vnhwm(nlt,nlat,nt) + input(17);
 wnhwm(nlt,nlat,nt) = wnhwm(nlt,nlat,nt) + input(18);
 den(nlt,nlat,nt) = den(nlt,nlat,nt) + input(9);
 denmd(nlt,nlat,nt) = denmd(nlt,nlat,nt) + input(6);
 denj8(nlt,nlat,nt) = denj8(nlt,nlat,nt) + input(7);
 denhd(nlt,nlat,nt) = denhd(nlt,nlat,nt) + input(8);
 denmd2(nlt,nlat,nt) = denmd2(nlt,nlat,nt) + input(5);
 den2(nlt,nlat,nt) = den2(nlt,nlat,nt) + input(9)*input(5)/input(6);
 denj82(nlt,nlat,nt) = denj82(nlt,nlat,nt) + input(7)*input(5)/input(6);
 denhd2(nlt,nlat,nt) = denhd2(nlt,nlat,nt) + input(8)*input(5)/input(6);
 ht(nlt,nlat,nt) = ht(nlt,nlat,nt) + input(1);
 nt = ntsv;
 end	% if (mt>0)
 numad(nad) = numad(nad) + 1;
 dun(nad) = dun(nad) + input(10);
 dvn(nad) = dvn(nad) + input(11);
 dwn(nad) = dwn(nad) + input(12);
 dunhwm(nad) = dunhwm(nad) + input(16);
 dvnhwm(nad) = dvnhwm(nad) + input(17);
 dwnhwm(nad) = dwnhwm(nad) + input(18);
 dden(nad) = dden(nad) + input(9);
 ddenmd(nad) = ddenmd(nad) + input(6);
 ddenhd(nad) = ddenhd(nad) + input(8);
 ddenj8(nad) = ddenj8(nad) + input(7);
 dden2(nad) = dden2(nad) + input(9)*input(5)/input(6);
 ddenhd2(nad) = ddenhd2(nad) + input(8)*input(5)/input(6);
 ddenj82(nad) = ddenj82(nad) + input(7)*input(5)/input(6);
 ddenmd2(nad) = ddenmd2(nad) + input(5);
 if (input(1)<nht)
  nht = input(1);
 end
 if (input(1)>xht)
  xht = input(1);
 end
end	%   if (year==2007 | year==2008)   (otherwise is set to be 999)
end	%   if (ifskp==0)	(for big gap Jan 6-8)
end	% for nlt=1:24
end	% for n10s=1:360

  fprintf(OUTG,'%4.4d %3.3d %4.4d %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %12.4e %12.4e %12.4e %12.4e %12.4e %12.4e %12.4e %12.4e %4.4d %12.4e %12.4e %12.4e %12.4e %12.4e %12.4e %12.4e %12.4e \n', ...
iyr,jda,numad(1),dun(1)/numad(1),dvn(1)/numad(1),dwn(1)/numad(1),dunhwm(1)/numad(1),dvnhwm(1)/numad(1),dwnhwm(1)/numad(1),dden(1)/numad(1),ddenmd(1)/numad(1), ...
ddenj8(1)/numad(1),ddenhd(1)/numad(1),dden2(1)/numad(1),ddenmd2(1)/numad(1),ddenj82(1)/numad(1),ddenhd2(1)/numad(1),numad(2),dden(2)/numad(2),ddenmd(2)/numad(2), ...
ddenj8(2)/numad(2),ddenhd(2)/numad(2),dden2(2)/numad(2),ddenmd2(2)/numad(2),ddenj82(2)/numad(2),ddenhd2(2)/numad(2));

if (mt==1)
 nhi = nhi + 1;
end
if (mt==2)
 nlo = nlo + 1;
end

end	% for nda=325:386
nht,xht




% Print out 2 files of den over 25 days 
for nt=1:3
 if nt==1
  outfile = sprintf('champ_density/hi_den_lat_vs_lt_accelb')
 end
 if nt==2
  outfile = sprintf('champ_density/lo_den_lat_vs_lt_accelb')
 end
 if nt==3
  outfile = sprintf('champ_density/all_den_lat_vs_lt_accelb')
 end
OUT = fopen(outfile,'w');
for nh=1:24
 for nl=1:60
  if (num(nh,nl,nt)>0)
    un(nh,nl,nt) = un(nh,nl,nt)/num(nh,nl,nt);
    vn(nh,nl,nt) = vn(nh,nl,nt)/num(nh,nl,nt);
    wn(nh,nl,nt) = wn(nh,nl,nt)/num(nh,nl,nt);
    unhwm(nh,nl,nt) = unhwm(nh,nl,nt)/num(nh,nl,nt);
    vnhwm(nh,nl,nt) = vnhwm(nh,nl,nt)/num(nh,nl,nt);
    wnhwm(nh,nl,nt) = wnhwm(nh,nl,nt)/num(nh,nl,nt);
    den(nh,nl,nt) = den(nh,nl,nt)/num(nh,nl,nt);
    denmd(nh,nl,nt) = denmd(nh,nl,nt)/num(nh,nl,nt);
    denj8(nh,nl,nt) = denj8(nh,nl,nt)/num(nh,nl,nt);
    denhd(nh,nl,nt) = denhd(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);
    denj82(nh,nl,nt) = denj82(nh,nl,nt)/num(nh,nl,nt);
    denhd2(nh,nl,nt) = denhd2(nh,nl,nt)/num(nh,nl,nt);
  end

  fprintf(OUT,'%5.1f %6.1f %d %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %12.4e %12.4e %12.4e %12.4e %12.4e %12.4e %12.4e %12.4e \n', ...
   lt24(nh),lat(nl),num(nh,nl,nt),un(nh,nl,nt),vn(nh,nl,nt),wn(nh,nl,nt), ...
   unhwm(nh,nl,nt),vnhwm(nh,nl,nt),wnhwm(nh,nl,nt),den(nh,nl,nt),denmd(nh,nl,nt), ...
   denj8(nh,nl,nt),denhd(nh,nl,nt),den2(nh,nl,nt),denmd2(nh,nl,nt),denj82(nh,nl,nt),denhd2(nh,nl,nt));

% Set 0's to NaN for plots
  if (num(nh,nl,nt)==0)
   un(nh,nl,nt) = NaN;
   vn(nh,nl,nt) = NaN;
   wn(nh,nl,nt) = NaN;
   unhwm(nh,nl,nt) = NaN;
   vnhwm(nh,nl,nt) = NaN;
   wnhwm(nh,nl,nt) = NaN;
   den(nh,nl,nt) = NaN;
   denmd(nh,nl,nt) = NaN;
   denj8(nh,nl,nt) = NaN;
   denhd(nh,nl,nt) = NaN;
   den2(nh,nl,nt) = NaN;
   denmd2(nh,nl,nt) = NaN;
   denj82(nh,nl,nt) = NaN;
   denhd2(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;
   denj8(nh,nl,nt) = denj8(nh,nl,nt)*1.e+13;
   denhd(nh,nl,nt) = denhd(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;
   denj82(nh,nl,nt) = denj82(nh,nl,nt)*1.e+13;
   denhd2(nh,nl,nt) = denhd2(nh,nl,nt)*1.e+13;
  end
 end	% nl=1:60
end	% nh=1:24   

end	% nt=nt1:nt2

% Cannot go from mid-pt to one extra edge for plotting - bad dimensions!
%  But still only see 23 hours, not 24??? on plots?? is because is edging?
%lt24 = 0:1:24;
%lat = -90:3:90;

unx=180.;
unn=-unx;
vnx=30.;
vnn=-vnx;
wnx=10.;
wnn=-wnx;

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

% Plot lo and hi Kp days together
kmsattit = sprintf('at ~348km 07325-08020');
km400tit = sprintf('at 400km 07325-08020');
hitit = sprintf('Kp>=2 Vsw>=500km/s');
lotit = sprintf('Kp<=1 Vsw<=450km/s');

tit1ave = [hitit,' CHAMP Doornbos density*1.e+13 in kg/m3 ',km400tit];
tit1med = [hitit,' MSIS density*1.e+13 in kg/m3 ',km400tit];
tit2ave = [lotit,' CHAMP Doornbos density*1.e+13 in kg/m3 ',km400tit];
tit2med = [lotit,' MSIS density*1.e+13 in kg/m3 ',km400tit];
tit1un = [hitit,' CHAMP Un(+E) m/s ',kmsattit];
tit2un = [lotit,' CHAMP Un(+E) m/s ',kmsattit];
tit1vn = [hitit,' CHAMP Vn(+N) m/s ',kmsattit];
tit2vn = [lotit,' CHAMP Vn(+N) m/s ',kmsattit];
tit1wn = [hitit,' CHAMP Wn(+up) m/s ',kmsattit];
tit2wn = [lotit,' CHAMP Wn(+up) m/s ',kmsattit];
tit1unh = [hitit,' HWM07 Un(+E) m/s ',kmsattit];
tit2unh = [lotit,' HWM07 Un(+E) m/s ',kmsattit];
tit1vnh = [hitit,' HWM07 Vn(+N) m/s ',kmsattit];
tit2vnh = [lotit,' HWM07 Vn(+N) m/s ',kmsattit];
tit1wnh = [hitit,' HWM07 Wn(+up) m/s ',kmsattit];
tit2wnh = [lotit,' HWM07 Wn(+up) m/s ',kmsattit];


figure
clf
subplot(2,1,1)
            pcolor(lt24,lat,un(:,:,1)')
            title([tit1un]);
            caxis([unn unx]);
            shading flat
            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,un(:,:,2)')
            title([tit2un]);
            caxis([unn unx]);
            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
clf
subplot(2,1,1)
            pcolor(lt24,lat,unhwm(:,:,1)')
            title([tit1unh]);
            caxis([unn unx]);
            shading flat
            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,unhwm(:,:,2)')
            title([tit2unh]);
            caxis([unn unx]);
            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
clf
subplot(2,1,1)
            pcolor(lt24,lat,vn(:,:,1)')
            title([tit1vn]);
            caxis([vnn vnx]);
            shading flat
            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,vn(:,:,2)')
            title([tit2vn]);
            caxis([vnn vnx]);
            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
clf
subplot(2,1,1)
            pcolor(lt24,lat,vnhwm(:,:,1)')
            title([tit1vnh]);
            caxis([vnn vnx]);
            shading flat
            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,vnhwm(:,:,2)')
            title([tit2vnh]);
            caxis([vnn vnx]);
            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
clf
subplot(2,1,1)
            pcolor(lt24,lat,wn(:,:,1)')
            title([tit1wn]);
            caxis([wnn wnx]);
            shading flat
            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,wn(:,:,2)')
            title([tit2wn]);
            caxis([wnn wnx]);
            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
clf
subplot(2,1,1)
            pcolor(lt24,lat,wnhwm(:,:,1)')
            title([tit1wnh]);
            caxis([wnn wnx]);
            shading flat
            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,wnhwm(:,:,2)')
            title([tit2wnh]);
            caxis([wnn wnx]);
            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
clf
subplot(2,1,1)
            pcolor(lt24,lat,den2(:,:,1)')
            title([tit1ave]);
            caxis([nden xden]);
            shading flat
            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,den2(:,:,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
clf
subplot(2,1,1)
            pcolor(lt24,lat,denmd2(:,:,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,denmd2(:,:,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;	


% Plot All CHAMP densities and winds with HWM 
titavu = ['All CHAMP av Un(m/s) at ~336-368km (Doornbos) 07325-08020'];
titavuh = ['All av Un(m/s) from HWM07 at ~336-368km 07325-08020'];
titavv = ['All CHAMP av Vn(m/s) at ~336-368km (Doornbos) 07325-08020'];
titavvh = ['All av Vn(m/s) from HWM07 at ~336-368km 07325-08020'];
titavw = ['All CHAMP av Wn(m/s) at ~336-368km (Doornbos) 07325-08020'];
titavwh = ['All av Wn(m/s) from HWM07 at ~336-368km 07325-08020'];
unx=300.;
unn=-unx;
vnx=30.;
vnn=-vnx;
wnx=10.;
wnn=-wnx;

% Plot HWM07 Un
figure
hold on;
            pcolor(lt24,lat,unhwm(:,:,3)')
            title([titavuh],'FontSize',fsize,'Fontweight','bold');
            caxis([unn unx]);
            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
subplot(2,1,1)
            pcolor(lt24,lat,un(:,:,3)')
            title([titavu]);
            caxis([unn unx]);
            shading flat
            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,unhwm(:,:,3)')
            title([titavuh]);
            caxis([unn unx]);
            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
subplot(2,1,1)
            pcolor(lt24,lat,vn(:,:,3)')
            title([titavv]);
            caxis([vnn vnx]);
            shading flat
            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,vnhwm(:,:,3)')
            title([titavvh]);
            caxis([vnn vnx]);
            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
subplot(2,1,1)
            pcolor(lt24,lat,wn(:,:,3)')
            title([titavw]);
            caxis([wnn wnx]);
            shading flat
            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,wnhwm(:,:,3)')
            title([titavwh]);
            caxis([wnn wnx]);
            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

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

titave = ['All CHAMP av density*1.e+13 in kg/m3 at ~336-368km (Doornbos) 07325-08020'];
titavmsis = ['All MSIS00 av density*1.e+13 in kg/m3 at ~336-368km 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 (Doornbos) 07325-08020'];
titavmsis = ['All MSIS00 av density*1.e+13 in kg/m3 at 400km 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;
	
figure
titave = ['All HASDM av density*1.e+13 in kg/m3 at 400km 07325-08020'];
titavmsis = ['All JB-2008 av density*1.e+13 in kg/m3 at 400km 07325-08020'];
subplot(2,1,1)
            pcolor(lt24,lat,denhd2(:,:,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,denj82(:,:,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;	
	
