clear all;
close all;

% Nov 2012  Program to calculate the ave of older legacy CHAMP winds from 
%  Kathrin Haeusler (using calc of spherical earth so alt ~5km lower)
%    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)
%  Revise to separate days of various Ap and F107 levels
% Kp1=0-1Ap=0-4(5), Kp2=2-4-Ap=7-22(27); F1=<75,F2=75-95,F3=95-125,F4=125-175,F5>175
%  10 cases:  Kp1F1, Kp2F1, Kp1F2, Kp2F2, Kp1F3, Kp2F3, Kp1F4, Kp2F4, Kp1F5, Kp2F5
%  4 seasons of 131 days:  NHwin=290-055; NHspr=15-145; NHsum=107-237; NHfal=199-329
fsize = 12;

% MJD = model Julian Day: First=0517=1 June 2001 (01152) and Last=3285=08364
% MJD 001: 2001=366 2002=731 2003=1096, 2004=1461, 2005=1827, 2006=2192, 2007=2557, 2008=2922
%                 sum01(136-179), fal01(136-236), win01(200-228)
% spr02(119-139), sum02(153-191), fal02(165-188), win02(121-165)
% spr03(118-139), sum03(119-133), fal03(113-150), win03(104-150)
% spr04(102-111), sum04(100-122), fal04(104-122), win04(92-105)
% spr05(87-102),  sum05(87-102),  fal05(75-88),   win05(75-88)
% spr06(75-90),   sum06(78-90),   fal06(74-84),   win06(74-84)
% spr07(71-81),   sum07(68-71),   fal07(67-74),   win07(67-72)
% spr08(69-72),   sum08(68-71),   fal08(67-68),   win08(67-68)
% spr09(68-72),   sum09(69-71),   fal09(69-72),   win09(72-83)
% spr10(75-83)

mjd01=[366 731 1096 1461 1827 2192 2557 2922];
apskp1=[10593 10958 11323 11688 12054 12419 12784 13149];
longd1=[15 107 199 290];	% for ndt=131 days for full LT orbit
shortd1=[35 127 219 310];	% for ndt=91 days for shorter season

% 1 Jan 2007=2557
iyear0 = 2007;
ndpy = 365;
ndt = 61;
ndt = 131;
% mjd=3156-3285 (2922=1 Jan 2008) or 235-364 but start at 290 (mjd=3211-85)
iyear0 = 2008;
ndt = 75;
% 1 Jan 2006=2192, mjd=2481-2611 (06290-07055), missing 2517, 2536, 2542-6, 2552-3, 2561, 2579, 2585, 2587-8, 2597, 2599-2601
%iyear0 = 2006;
%ndt = 131;
% 1 jan 2005 = 1827, mjd=2116-2246 (05290-06055), missing 2147-9, 2154-6, 2165, 2191, 2209, 2217, 2240
iyear0 = 2005;
ndt = 131;

miss05=[2147 2148 2149 2154 2155 2156 2165 2191 2209 2217 2240]
miss06=[2484 2517 2536 2542 2543 2544 2545 2546 2552 2553 2561 2579 2585 2587 2588 2597 2599 2600 2601];



% Get daily av Kp(5), Ap(8), Vsw(37) from eln_tot_72_11_combnx7_daily (daynum=2)
filap = sprintf('/hao/polar1/emery/aurora_i/xuw/combined_output_daily/eln_tot_72_11_combnx7_daily')
INAP = fopen(filap,'r');
if (iyear0==2007)
 nskp = 13073;	% to get to 07290
 if (ndt==61)
  nskp = nskp + 35;
 end
end
if (iyear0==2008)
% mjd=3156-3285 (2922=1 Jan 2008) or 235-364 but start at 290 (mjd=3211-85)
 nskp = 13073+365;
end
if (iyear0==2006)
 nskp=13073-365;
end
if (iyear0==2005)
 nskp=13073-2*365;
end
% Skip over lines
for n=1:nskp
     rdhd = fgetl(INAP);
end

if (ndt==131 & iyear0==2005)
 outglob=sprintf('champ_density/global_LTden_haeusler_05_06.dat')
 wrap = sprintf('champ_density/daily_glon_2LTs_05_06')
end
if (ndt==131 & 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 (ndt==61)
% Name 'global' den 400km output file for daily global means (from equal wt ascent/descent)
outglob=sprintf('champ_density/global_LTden_haeusler_07325_08020.dat')
 wrap = sprintf('champ_density/daily_glon_2LTs_07325_08020b')
end
if (ndt==131 & 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 (ndt==75)
 outglob=sprintf('champ_density/global_LTden_haeusler_08.dat')
 wrap = sprintf('champ_density/daily_glon_2LTs_08')
end
 OUTG = fopen(outglob,'w');
 OUTAP = fopen(wrap,'w');
  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');
% CHAMP densities missing for days 2008 Jan 6-8 and for Dec 22 (356) 2007

% Midpoints 24h and 60 lats
lt24 = 0.5:1:23.5;
lat = -88.5:3:88.5;
lon=[0:15:355];
nlon = length(lon);

% 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, 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, 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.;
   un(i,j,k) = 0.;
   vn(i,j,k) = 0.;
   wn(i,j,k) = 0.;
   ht(i,j,k) = 0.;
  end
 end
end
for nd=1:ndt
for k=1:2
% numad = number of ascent (8-3LT) k=1 and descent (20-15LT) k=2 times per day (d..=ad)
 numad(nd,k) = 0;
 dden(nd,k) = 0.;
 dun(nd,k) = 0.;
 dvn(nd,k) = 0.;
 dwn(nd,k) = 0.;
end
end

neqlt = 24;
xeqlt = 0;

nt = 3;
nhi = 1;
nlo = 1;
mt = 0;

nht = 5.e+7;
xht = 0.;

% Kathrin's file
%  datkh=dlmread('/hao/aim3/emery/cetic_climatology/champ_density/khwinds/kh_CHAMP_07325-08020.alldat','',0,0);
%nll = length(datkh)	% = 482692
%for n=1:nll
%if (n==275)
%datkh(n,:)
%end
%rmjd = datkh(n,1);
%mjd = fix(rmjd);
% Problem with fix(rmjd)=2880 instead of 2881 at n=275 - break into days!

if (ndt==61)
 nda1 = 325;
else
 nda1 = 290;
end
nda2 = nda1+ndt-1

for nda=nda1:nda2
iyr = iyear0;
jda = nda
if (nda>ndpy)
 iyr = iyear0+1;
 jda = nda - ndpy;
end
% Read indices file
    rdap = fgetl(INAP);
    datap = sscanf(rdap,'%f');
    kpav = datap(5)
    apav = datap(8);
    vswav = datap(37);
    f107 = datap(10);
% Prep arrays for daily values
  undlon(1:nlon,1:2) = 0.;
  d400dlon(1:nlon,1:2) = 0.;
  numdlon(1:nlon,1:2) = 0.;
  numd2 = 0;
  numd(1:2) = 0;
  stld(1:2) = 0.;
  htd(1:2) = 0.;
clear datkh;
ifskp = 0;
% Missing 2860, 2874
if (iyear0==2007)
if (nda==304 | nda==318)
 ifskp=1;
end
if (nda==356)
 ifskp=1;
end
if (nda>370 & nda<374)
 ifskp=1;
end
mjd = 2881 + nda-325;
end	% if iyear0==2007
if (iyear0==2008)
 mjd = nda + 2921;
end
if (iyear0==2006)
 mjd = nda + 2191;
 for nm=1:19
  if (mjd==miss06(nm))
    ifskp=1;
  end
 end
end
if (iyear0==2005)
 mjd = nda + 1826;
 for nm=1:11
  if (mjd==miss05(nm))
    ifskp=1;
  end
 end
end

if (ifskp==0)

nd = nda - nda1 + 1;
mt = 0;
if (nda==ndahi(nhi))
 mt = 1;
end
if (nda==ndalo(nlo))
 mt = 2;
end

%if (ndt==131 & iyear0==2005)
%filin=sprintf('/hao/aim3/emery/cetic_climatology/champ_density/khwinds/CHAMP_Barbara/CHAMP_winter_05:06/CHAMP_MJD_%d.txt',mjd);
%end
if (ndt==61)
filin=sprintf('/hao/aim3/emery/cetic_climatology/champ_density/khwinds/CHAMP_07325-08020/%d.dat',mjd);
else
filin=sprintf('/hao/aim3/emery/cetic_climatology/champ_density/khwinds/MJD/CHAMP_MJD_%d.txt',mjd);
end
datkh = dlmread(filin,'',0,0);
nll = length(datkh)
for n=1:nll

% 13 inputs: mjd,slt,mlt,glat,glat360,glon,mlat,mlon,alt,den,Vn,Un,Wn
 nlt = fix(datkh(n,2)) + 1;

% Choose first LT from first of each day - can sort out later
 
% Accumulate in daily arrays of glon at 2 LTs
  if (datkh(n,4)<13.1 & datkh(n,4)>-13.1)
   numd2 = numd2 + 1;
   if (numd2==1)
% Get first value of day
    stld(1) = datkh(n,2);
    htd(1) = datkh(n,9) - 6371.;
   end
   diflt = abs(stld(1)-datkh(n,2));
   if (diflt<1. | diflt>22.)
    mad = 1;
   else
    mad = 2;
% Get near last value of day
    stld(2) = datkh(n,2);
    htd(2) = datkh(n,9) - 6371.;
   end
   glon = datkh(n,6);
   if (glon<0.)
    glon = glon + 360.;
   end
   jl = fix(glon/15.) + 1;
   numdlon(jl,mad) = numdlon(jl,mad) + 1;
   numd(mad) = numd(mad) + 1;
   undlon(jl,mad) = undlon(jl,mad) + datkh(n,12);
   d400dlon(jl,mad) = d400dlon(jl,mad) + datkh(n,10);
  end	%  if (datkh(n,4)<13.1 & datkh(n,4)>-13.1)

 nlat = fix((datkh(n,4)+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) + datkh(n,12);
 vn(nlt,nlat,nt) = vn(nlt,nlat,nt) + datkh(n,11);
 wn(nlt,nlat,nt) = wn(nlt,nlat,nt) + datkh(n,13);
 den(nlt,nlat,nt) = den(nlt,nlat,nt) + datkh(n,10);
 if (mt>0)
  ntsv = nt;
  nt = mt;
 num(nlt,nlat,nt) = num(nlt,nlat,nt) + 1;
 un(nlt,nlat,nt) = un(nlt,nlat,nt) + datkh(n,12);
 vn(nlt,nlat,nt) = vn(nlt,nlat,nt) + datkh(n,11);
 wn(nlt,nlat,nt) = wn(nlt,nlat,nt) + datkh(n,13);
 den(nlt,nlat,nt) = den(nlt,nlat,nt) + datkh(n,10);
 nt = ntsv;
 end	% if (mt>0)

% Choose nad old way
% nad = 1 for ascent (8.3 to 2.8LT) and 2 for descent (20.4 to 14.8LT) 
 nad = 2;
if (datkh(n,2)<12.)
 nad = 1;
end

 numad(nd,nad) = numad(nd,nad) + 1;
 dun(nd,nad) = dun(nd,nad) + datkh(n,12);
 dvn(nd,nad) = dvn(nd,nad) + datkh(n,11);
 dwn(nd,nad) = dwn(nd,nad) + datkh(n,13);
 dden(nd,nad) = dden(nd,nad) + datkh(n,10);


end	% for n=1,nll
if (mt==1)
 nhi = nhi + 1;
end
if (mt==2)
 nlo = nlo + 1;
end
% Print-out of daily arrays of glon at 2 LTs 
  for mad=1:2
   for jl=1:nlon
    undlon(jl,mad) = undlon(jl,mad)/numdlon(jl,mad);
    d400dlon(jl,mad) = d400dlon(jl,mad)/numdlon(jl,mad);
   end
  end
  fprintf(OUTAP,'%4.4d %3.3d %8.2f %8.2f %8.2f %8.2f %d %8.2f %8.2f %d %8.2f %8.2f \n', iyr,jda,f107,kpav,apav,vswav,numd(1),stld(1),htd(1),numd(2),stld(2),htd(2));
  fprintf(OUTAP,'%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',undlon(:,1));
  fprintf(OUTAP,'%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',undlon(:,2));
  fprintf(OUTAP,'%e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e \n',d400dlon(:,1));
  fprintf(OUTAP,'%e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e \n',d400dlon(:,2));
end	% if (ifskp==0)
end	% for nda=nda1:nda2

nd = 0;
for nda=nda1:nda2
 nd=nd+1;
 iyr = iyear0;
 jda = nda;
  if (nda>ndpy)
   iyr = iyear0+1;
   jda = nda-ndpy;
  end
  fprintf(OUTG,'%4.4d %3.3d %4.4d %8.2f %8.2f %8.2f %12.4e %4.4d %8.2f %8.2f %8.2f %12.4e \n', iyr,jda,numad(nd,1),dun(nd,1)/numad(nd,1),dvn(nd,1)/numad(nd,1),dwn(nd,1)/numad(nd,1),dden(nd,1)/numad(nd,1),numad(nd,2),dun(nd,2)/numad(nd,2),dvn(nd,2)/numad(nd,2),dwn(nd,2)/numad(nd,2),dden(nd,2)/numad(nd,2));
end	% for nda=nda1:nda2


% Print out 2 files of den over 25 days 
for nh=1:24
 for ilat=1:9
    num12(nh,ilat) = 0;
    un12(nh,ilat) = 0.;
    d400(nh,ilat) = 0.;
 end
end
%for nt=1:3
for nt=3:3
 if nt==1
  outfile = sprintf('champ_density/hi_den_lat_vs_lt_kh2')
 end
 if nt==2
  outfile = sprintf('champ_density/lo_den_lat_vs_lt_kh2')
 end
 if nt==3
  outfile = sprintf('champ_density/all_den_lat_vs_lt_kh_%d',iyear0)
%  outcnofs = sprintf('champ_winds_cnofs_vs_lt_%d',iyear0)
  outcnofs = sprintf('champ_winds_9glat_vs_lt_290-055_%d',iyear0)
  OUTC = fopen(outcnofs,'w');
  fprintf (OUTC,'SLT 9(# Un(m/s) NDen(kg/m3)) for glat=-60:15:60 \n');
 end
OUT = fopen(outfile,'w');
% vn and wn are components of the cross-track obs in the +N and +up directions (are small so ignore them)
for nh=1:24

 for nl=1:60
  if (num(nh,nl,nt)>0)

if (nt==3)
% 60 3-deg glats from -88.5:3:88.5 (collect 15deg -60:15:60)
   if (nl>=27 & nl<=34) 
% Average over nl=27-34 (+/-12glat) for C/NOFS
    ilat = 5;
    num12(nh,ilat) = num12(nh,ilat) + num(nh,nl,nt);
    un12(nh,ilat) = un12(nh,ilat) + un(nh,nl,nt);
    d400(nh,ilat) = d400(nh,ilat) + den(nh,nl,nt);
   end
   if (nl>=23 & nl<=28) 
    ilat = 4;
    num12(nh,ilat) = num12(nh,ilat) + num(nh,nl,nt);
    un12(nh,ilat) = un12(nh,ilat) + un(nh,nl,nt);
    d400(nh,ilat) = d400(nh,ilat) + den(nh,nl,nt);
   end
   if (nl>=18 & nl<=23) 
    ilat = 3;
    num12(nh,ilat) = num12(nh,ilat) + num(nh,nl,nt);
    un12(nh,ilat) = un12(nh,ilat) + un(nh,nl,nt);
    d400(nh,ilat) = d400(nh,ilat) + den(nh,nl,nt);
   end
   if (nl>=13 & nl<=18) 
    ilat = 2;
    num12(nh,ilat) = num12(nh,ilat) + num(nh,nl,nt);
    un12(nh,ilat) = un12(nh,ilat) + un(nh,nl,nt);
    d400(nh,ilat) = d400(nh,ilat) + den(nh,nl,nt);
   end
   if (nl>=8 & nl<=13) 
    ilat = 1;
    num12(nh,ilat) = num12(nh,ilat) + num(nh,nl,nt);
    un12(nh,ilat) = un12(nh,ilat) + un(nh,nl,nt);
    d400(nh,ilat) = d400(nh,ilat) + den(nh,nl,nt);
   end
   if (nl>=33 & nl<=38) 
    ilat = 6;
    num12(nh,ilat) = num12(nh,ilat) + num(nh,nl,nt);
    un12(nh,ilat) = un12(nh,ilat) + un(nh,nl,nt);
    d400(nh,ilat) = d400(nh,ilat) + den(nh,nl,nt);
   end
   if (nl>=38 & nl<=43) 
    ilat = 7;
    num12(nh,ilat) = num12(nh,ilat) + num(nh,nl,nt);
    un12(nh,ilat) = un12(nh,ilat) + un(nh,nl,nt);
    d400(nh,ilat) = d400(nh,ilat) + den(nh,nl,nt);
   end
   if (nl>=43 & nl<=48) 
    ilat = 8;
    num12(nh,ilat) = num12(nh,ilat) + num(nh,nl,nt);
    un12(nh,ilat) = un12(nh,ilat) + un(nh,nl,nt);
    d400(nh,ilat) = d400(nh,ilat) + den(nh,nl,nt);
   end
   if (nl>=48 & nl<=53) 
    ilat = 9;
    num12(nh,ilat) = num12(nh,ilat) + num(nh,nl,nt);
    un12(nh,ilat) = un12(nh,ilat) + un(nh,nl,nt);
    d400(nh,ilat) = d400(nh,ilat) + den(nh,nl,nt);
   end
end	% if (nt==3)

    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);
    den(nh,nl,nt) = den(nh,nl,nt)/num(nh,nl,nt);

  end	% if (num(nh,nl,nt)>0)

  fprintf(OUT,'%5.1f %6.1f %d %8.2f %8.2f %8.2f %12.4e \n', ...
   lt24(nh),lat(nl),num(nh,nl,nt),un(nh,nl,nt),vn(nh,nl,nt),wn(nh,nl,nt),den(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;
   den(nh,nl,nt) = NaN;
  else
   den(nh,nl,nt) = den(nh,nl,nt)*1.e+13;
  end
 end	% nl=1:60

% Print out C/NOFS info as func of LT only
 if (nt==3)
  for ilat=1:9
  if (num12(nh,ilat)>0);
    un12(nh,ilat) = un12(nh,ilat)/num12(nh,ilat);
    d400(nh,ilat) = d400(nh,ilat)/num12(nh,ilat);
  end
  end
   fprintf(OUTC,'%5.1f %d %8.2f %12.4e %d %8.2f %12.4e %d %8.2f %12.4e %d %8.2f %12.4e %d %8.2f %12.4e %d %8.2f %12.4e %d %8.2f %12.4e %d %8.2f %12.4e %d %8.2f %12.4e \n',lt24(nh),num12(nh,1),un12(nh,1),d400(nh,1),num12(nh,2),un12(nh,2),d400(nh,2),num12(nh,3),un12(nh,3),d400(nh,3),num12(nh,4),un12(nh,4),d400(nh,4),num12(nh,5),un12(nh,5),d400(nh,5),num12(nh,6),un12(nh,6),d400(nh,6),num12(nh,7),un12(nh,7),d400(nh,7),num12(nh,8),un12(nh,8),d400(nh,8),num12(nh,9),un12(nh,9),d400(nh,9));
 end	% if (nt==3)
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=300.;
unn=-300.;
vnx=20.;
vnn=-vnx;
wnx=6.;
wnn=-wnx;

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

% Plot lo and hi Kp days together
if (ndt==131 & iyear0==2005)
 kmsattit = sprintf('at ~348km 05290-06055');
 km400tit = sprintf('at 400km 05290-06055');
end
if (ndt==131 & iyear0==2006)
 kmsattit = sprintf('at ~348km 06290-07055');
 km400tit = sprintf('at 400km 06290-07055');
end
if (ndt==61)
 kmsattit = sprintf('at ~348km 07325-08020');
 km400tit = sprintf('at 400km 07325-08020');
end
if (ndt==131 & iyear0==2007)
 kmsattit = sprintf('at ~348km 07290-08055');
 km400tit = sprintf('at 400km 07290-08055');
end
if (ndt==75)
 kmsattit = sprintf('at ~348km 08290-0864');
 km400tit = sprintf('at 400km 08290-08364');
end
hitit = sprintf('Kp>=2 Vsw>=500km/s');
lotit = sprintf('Kp<=1 Vsw<=450km/s');

tit1ave = [hitit,' CHAMP Haeusler density*1.e+13 in kg/m3 ',kmsattit];
tit1un = [hitit,' KH CHAMP Un(+E) m/s ',kmsattit];
tit2un = [lotit,' KH CHAMP Un(+E) m/s ',kmsattit];
tit1vn = [hitit,' KH CHAMP Vn(+N) m/s ',kmsattit];
tit2vn = [lotit,' KH CHAMP Vn(+N) m/s ',kmsattit];
tit1wn = [hitit,' KH CHAMP Wn(+up) m/s ',kmsattit];
tit2wn = [lotit,' KH CHAMP Wn(+up) m/s ',kmsattit];


figure
clf
subplot(2,1,1)
            pcolor(lt24,lat,un(:,:,1)')
            title([tit1un]);
            caxis([unn unx]);
            shading flat
	    grid;
            ylabel('Latitude');
            xlabel('LT');
	    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');
subplot(2,1,2)
            pcolor(lt24,lat,un(:,:,2)')
            title([tit2un]);
            caxis([unn unx]);
            shading flat
	    grid;
            xlabel('LT');
            ylabel('Latitude');
	    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
clf
subplot(2,1,1)
            pcolor(lt24,lat,vn(:,:,1)')
            title([tit1vn]);
            caxis([vnn vnx]);
            shading flat
	    grid;
            ylabel('Latitude');
	    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');
subplot(2,1,2)
            pcolor(lt24,lat,vn(:,:,2)')
            title([tit2vn]);
            caxis([vnn vnx]);
            shading flat
	    grid;
            xlabel('LT');
            ylabel('Latitude');
	    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
clf
subplot(2,1,1)
            pcolor(lt24,lat,wn(:,:,1)')
            title([tit1wn]);
            caxis([wnn wnx]);
            shading flat
	    grid;
            ylabel('Latitude');
	    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');
subplot(2,1,2)
            pcolor(lt24,lat,wn(:,:,2)')
            title([tit2wn]);
            caxis([wnn wnx]);
            shading flat
	    grid;
            xlabel('LT');
            ylabel('Latitude');
	    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
clf
subplot(1,1,1)
            pcolor(lt24,lat,den(:,:,1)')
            title([tit1ave]);
            caxis([nden xden]);
            shading flat
	    grid;
            ylabel('Latitude');
	    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 CHAMP densities and winds with HWM 
titavu = ['All CHAMP Haeusler av Un (m/s) ',kmsattit];
titavv = ['All CHAMP Haeusler av Vn (m/s) ',kmsattit];
titavw = ['All CHAMP Haeusler av Wn (m/s) ',kmsattit];

figure
hold on;
            pcolor(lt24,lat,un(:,:,3)')
            title([titavu],'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
hold on;
            pcolor(lt24,lat,vn(:,:,3)')
            title([titavv]);
            caxis([vnn vnx]);
            shading flat
	    grid;
            ylabel('Latitude');
	    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,wn(:,:,3)')
            title([titavw]);
            caxis([wnn wnx]);
            shading flat
	    grid;
            ylabel('Latitude');
	    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;	


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

titave = ['All CHAMP Haeusler av density*1.e+13 in kg/m3 ',km400tit];
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;	



