%%
%% This program is written for comparing the TIEGCM and SON ISR data
%% BA Emery, April 2011
clear all; close all

% notie=1 if have no TIE values to read (wd and wded have no TN,UN,VN)
notie = 0
% madskp = 1 means skip get_mad for obs in all directions
madskp = 0
% if1m=1 plots 1-min Sonde Zenith data instead of 20-min medians
if1m=0
% iflp=1 plots long-pulse data
iflp = 1

 year1 = 2005;
 year1 = 2006;
if (year1 < 2005.5)
 filerms = sprintf('rms_may05_son_tie');
 day1 = 135;
 da1 = 135;
 day2 = 136;
 da2 = 135+ 20./24.;
% Both SP and LP data have same number of times (so can be combined easily)
 n15m = 207;
% ExB file has 1 less in 20 min in May05, and Dec05 last is missing so could be deleted
 n20m = 48;
else
 filerms = sprintf('rms_dec06_son_tie');
 day1 = 348;
 da1 = 348 + 12./24.;
 day2 = 350;
 da2 = 350 + 0./24.;
 n15m = 192;
 n20m = 71;
end
begut = [da1,da1];
endut = [da2,da2];

nf1 = 1;
it1nf1 = 1;
nf2 = 1;

 iday0 = day1 - 1;
% 20 min saves
 ntmpda = 24*3;
lenda = day2 - day1;
nobs = ntmpda*lenda + 1
% Subtract 1 time so start at 0020 UT on first day
nobs = nobs - 1

% Number of 20 min points to look at the median absolute deviation (mad)
numhr = 4;
npmad = numhr*3+1;
% midm1=6 points in 4h of 13 20min obs
midm1 = (npmad-1)/2;
% Set the min number of nmad points small because is sparse data 
%   Therefore, for sparse data, may only find mad of combined data (and set madskp=1)
minnmad = midm1;

% Choose number in 20-min bin above which bin is plotted in contourf or used at ht300
binnum = 1.01;
binnum = 0.01;
% Look at number for ht300 (n300=#24)
if (iflp==0)
% kindat=5521 short pulse data
nalt = 32;
n300 = 24;
ht300 = 243;
else
% kindat=5501 long pulse data zenith
nalt = 21;
n300 = 8;
ht300 = 298;
% 10 alts from 141-497 km
nalt = 10
n300 = 6;
% kindat=5501 long pulse data 75 deg elev
%nalt = 22;
%n300 = 9;
%ht300 = 309;
end
ht300str = sprintf ('%d km',ht300);


z300 = [ht300,ht300];
zero = [0.,0.];
missing = -32767;

% Look at altitude ht300 km or range of +/-23 km 
ht1 = (ht300-23.)*1.e+5;
ht2 = (ht300+23.)*1.e+5;

purple = [0.48, 0.06, 0.89];
dark_green = [0 0.5 0];
red_brown = [0.6 0.2 0];
yel_grn = [0.75 0.75 0];

% Read IMF file (TIEGCM run used hourly IMF) of 15-min aves every 5 min (t=0 from t=-5 to -20m)
 if (year1 < 2005.5)
   Fname=['/hao/aim1/tgcm/data/imf_OMNI_2005001-2005365.nc'];
 else
   Fname=['/hao/aim1/tgcm/data/imf_OMNI_2006001-2006365.nc'];
 end
  nc=netcdf(Fname,'nowrite');
  yrdayfracimf=nc{'date'}(:);
  yrbximf=nc{'bx'}(:);
  yrbyimf=nc{'by'}(:);
  yrbzimf=nc{'bz'}(:);
  yrvswimf=nc{'swvel'}(:);
  close(nc)
  clear nc 
  ntimf=length(yrdayfracimf);
  nt1 = (day1-1)*288 + 1
  nt2 = (day2-1)*288 + 1
  ntot = nt2-nt1 + 1
  for n=nt1:nt2
    dayfracimf(n-nt1+1) = yrdayfracimf(n)-year1*1000;
    bximf(n-nt1+1) = yrbximf(n);
    byimf(n-nt1+1) = yrbyimf(n);
    bzimf(n-nt1+1) = yrbzimf(n);
    btmag(n-nt1+1) = sqrt(yrbximf(n)^2+yrbyimf(n)^2+yrbzimf(n)^2);
    vsw(n-nt1+1) = yrvswimf(n);
  end
for id=1:lenda
 ni=288*id
 dayfracimf(ni),vsw(ni)
end


% Sondrestrom coordinates (66.99N, 50.95W, j,i=32,27)  36 lats from -87.5 to +87.5, 72 lons from -180 to +175(?)
jj = 32;
ii = 27;

% Sonde 298 km Bmag is 0.4884gauss vs apex 0.4865gauss, I=79.995, D=-31.867
deginc = 79.995;
degdec = -31.867;

% Have SON obs start at 0020 UT for the 20 min values to be the same as the TIEGCM

for n=1:nobs
 for j=1:nalt
  wiobs_20min(j,n) = NaN;
  tiobs_20min(j,n) = NaN;
  teobs_20min(j,n) = NaN;
  neobs_20min(j,n) = NaN;
  wiobsv_20min(j,n) = NaN;
  tiobsv_20min(j,n) = NaN;
  teobsv_20min(j,n) = NaN;
  neobsv_20min(j,n) = NaN;
 end
 wi_obs_nan_20min(n) = NaN;
 ti_obs_nan_20min(n) = NaN;
 te_obs_nan_20min(n) = NaN;
 ne_obs_nan_20min(n) = NaN;
 nmax_obs_nan_20min(n) = NaN;
 hmax_obs_nan_20min(n) = NaN;
 tec_obs_nan_20min(n) = NaN;
 vpe_obs_nan_20min(n) = NaN;
 vpn_obs_nan_20min(n) = NaN;
% Start at 0000 UT
%time_obs_20min(n) = day1 + (n-1)/ntmpda;
% Start at 0020 UT (20 min is .0139 frac_days)
 time_obs_20min(n) = day1 + n/ntmpda;
end

if (year1 < 2005.5)
if (iflp==0)
 if (if1m==1)
 obs1fil = ['/hao/aim2/emery/cetic11/son050514_15_nmax'];
 obs2fil = ['/hao/aim2/emery/cetic11/son050514_15_profiles'];
 end
 obs1m20 = ['/hao/aim2/emery/cetic11/son050514_15_20m_nmax'];
 obs2m20 = ['/hao/aim2/emery/cetic11/son050514_15_20m_profiles'];
 obs3m20 = ['/hao/aim2/emery/cetic11/son_epec_050514_15_20m_ExB'];
else
 if (if1m==1)
 obs1fil = ['/hao/aim2/emery/cetic11/son_lp050514_15_nmax'];
 obs2fil = ['/hao/aim2/emery/cetic11/son_lp050514_15_profiles'];
 end
% LP90
 obs1m20 = ['/hao/aim2/emery/cetic11/son_lp050514_15_20m_nmax'];
 obs2m20 = ['/hao/aim2/emery/cetic11/son_lp050514_15_20m_profiles'];
 obs3m20 = ['/hao/aim2/emery/cetic11/son_epec_050514_15_20m_ExB'];
% LP75
%obs1m20 = ['/hao/aim2/emery/cetic11/son_lp75_050514_15_20m_nmax'];
%obs2m20 = ['/hao/aim2/emery/cetic11/son_lp75_050514_15_20m_profiles'];
% LP75+Z
%obs1m20 = ['/hao/aim2/emery/cetic11/son_lpc_050514_15_20m_nmax'];
%obs2m20 = ['/hao/aim2/emery/cetic11/son_lpc_050514_15_20m_profiles'];
end
else
 obs1m20 = ['/hao/aim2/emery/cetic11/son_lp061214_15_20m_nmax'];
 obs2m20 = ['/hao/aim2/emery/cetic11/son_lp061214_15_20m_profiles'];
 obs3m20 = ['/hao/aim2/emery/cetic11/son_epec_061214_15_20m_ExB'];
end


if (if1m==1)
n = 0;
%%  Read SON obs and assign each time to closest 20min model time
OBS1=fopen(obs1fil,'r');
OBS2=fopen(obs2fil,'r');
% Skip over first line
 linobs1 = fgetl(OBS1);
 linobs2 = fgetl(OBS2);

for m=1:n15m
% Read 1D numbers
 linobs1 = fgetl(OBS1);
 obs1dat = sscanf(linobs1,'%f %f %f %f %f %f %f %f %f %f');
 yrobs = obs1dat(1);
 daynumobs = obs1dat(2);
% Do not save if daynumobs < day1 or > day2-1
if (daynumobs >= day1 & daynumobs <= day2-1) 
 n = n + 1;
 utobs = obs1dat(4);
 timeobs(n) = daynumobs + utobs/24.;
 hmaxobs(n) = NaN;
 nmaxobs(n) = NaN;
 tecobs(n) = NaN;
% if (obs1dat(7) > 99.)
 if (obs1dat(7) > 99. & obs1dat(8) > 0.281)
   hmaxobs(n) = obs1dat(7);
 end
% foF2=0.28 seems to be a minimum MHz value and is not valid - delete Nmax and hmax
%if (obs1dat(9) > 0.)
 if (obs1dat(8) > 0.281)
   nmaxobs(n) = log10(obs1dat(9));
 end
 if (obs1dat(10) > 0.)
   tecobs(n) = log10(obs1dat(10));
 end

% Read 2D profiles
 tiobs300(n) = NaN;
 teobs300(n) = NaN;
 neobs300(n) = NaN;
  for n9=1:nalt
   linobs2 = fgetl(OBS2);
   obs2dat = sscanf(linobs2,'%f %f %f %f %f %f %f %f %f %f');
   if (n == 1)
     altobs(n9) = obs2dat(4);
   end
   tiobs(n9,n) = NaN;
   teobs(n9,n) = NaN;
   neobs(n9,n) = NaN;
 if (obs2dat(5) > 0.)
  tiobs(n9,n) = obs2dat(5);
 end
 if (obs2dat(7) > 0.)
  teobs(n9,n) = obs2dat(7);
 end
 if (obs2dat(9) > 0.)
  neobs(n9,n) = obs2dat(9);
 end
if (n9 == n300)
 if (obs2dat(5) > 0.)
  tiobs300(n) = obs2dat(5);
 end
 if (obs2dat(7) > 0.)
  teobs300(n) = obs2dat(7);
 end
 if (obs2dat(9) > 0.)
  neobs300(n) = obs2dat(9);
 end

end	% end of if n9=n300
end	% end n9=1:nalt
end	% end of if daynumobs between day1 and day2-1
end	% end of for m=1:n15m
nm1 = n;
end 	% end of if1m==1



% Read 20 min observational values
OBS1M=fopen(obs1m20,'r');
OBS2M=fopen(obs2m20,'r');
OBS3M=fopen(obs3m20,'r');
% Skip over first line
 linobs1 = fgetl(OBS1M);
 lin2obs = fgetl(OBS2M);
 lin3obs = fgetl(OBS3M);
% Skip over data description 
if (year1 < 2005.5 & iflp==0)
% for n=1:48
%  linobs1 = fgetl(OBS1M);
% end
end

nvi = 6;
nti = 9;
nte = 12;
nne = 15;
day1,time_obs_20min(1)
for m=1:n20m
% Read n20m-1 1-D ExB numbers
if (m<n20m)
 linobs3 = fgetl(OBS3M);
 obs3dat = sscanf(linobs3,'%f %f %f %f %f %f %f %f %f %f');
end

% Read 1D numbers
 linobs1 = fgetl(OBS1M);
 obs1dat = sscanf(linobs1,'%f %f %f %f %f %f %f %f %f %f');
 yrobs = obs1dat(1);
 daynumobs = obs1dat(2);
 utobs = obs1dat(4);
 dafrac = daynumobs + utobs/24.;
% Find nearest 20 min time
 n = fix( ((dafrac-time_obs_20min(1))*1440.+9.99)/20. ) + 1;
% Do not save if daynumobs < day1(0020UT) or > day2-1
noskp=0;
if (daynumobs < day1 | daynumobs > day2-1 | n < 1 | n > nobs)
 noskp=1;
end
if (noskp==0)
 time_obs_20min(n) = daynumobs + utobs/24.;
 hmax_obs_nan_20min(n) = obs1dat(7);
 if (obs1dat(8) > 0.2)
  nmax_obs_nan_20min(n) = log10(obs1dat(9));
  tec_obs_nan_20min(n) = log10(obs1dat(10));
 end

% Assign ExB data from obs3dat
if (m<n20m & obs3dat(20)>-32760.)
 vpe_obs_nan_20min(n) = obs3dat(25);
 vpn_obs_nan_20min(n) = obs3dat(20);
end

% Read 2D profiles
  for n9=1:nalt
   lin2obs = fgetl(OBS2M);
%  obs2dat = sscanf(lin2obs,'%f %f %f %f %f %f %f %f %f %f');
%  obs2dat = sscanf(lin2obs,'%f %f %f %f %f %f %f %f %f %f %f %f %f');
   obs2dat = sscanf(lin2obs,'%f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f');
     altobs(n9) = obs2dat(4);
 if (obs2dat(nvi) > 0. & obs2dat(nvi-1) > binnum)
  wiobs_20min(n9,n) = obs2dat(nvi);
  wiobsv_20min(n9,n) = obs2dat(nvi+1);
 end
 if (obs2dat(nti) > 0. & obs2dat(nti-1) > binnum)
  tiobs_20min(n9,n) = obs2dat(nti);
  tiobsv_20min(n9,n) = obs2dat(nti+1);
 end
 if (obs2dat(nte) > 0. & obs2dat(nte-1) > binnum)
  teobs_20min(n9,n) = obs2dat(nte);
  teobsv_20min(n9,n) = obs2dat(nte+1);
 end
 if (obs2dat(nne) > 0. & obs2dat(nne-1) > binnum)
  neobs_20min(n9,n) = log10(obs2dat(nne));
  neobsv_20min(n9,n) = log10(obs2dat(nne+1));
 end
if (n9 == n300)
 if (obs2dat(nvi) > 0. & obs2dat(nvi-1) > binnum)
  wi_obs_nan_20min(n) = obs2dat(nvi);
 end
 if (obs2dat(nti) > 0. & obs2dat(nti-1) > binnum)
  ti_obs_nan_20min(n) = obs2dat(nti);
 end
 if (obs2dat(nte) > 0. & obs2dat(nte-1) > binnum)
  te_obs_nan_20min(n) = obs2dat(nte);
 end
 if (obs2dat(nne) > 0. & obs2dat(nne-1) > binnum)
  ne_obs_nan_20min(n) = log10(obs2dat(nne));
 end

end	% end of if n9=n300
end	% end n9=1:nalt
end	% end of if daynumobs between day1 and day2-1 (if noskp==0)
end	% end of for m=1:n20m
n,time_obs_20min(n)


% Calculate mad from tn_obs_20min (assuming obs start at 0020min - FALSE FOR 75elev!!)
nmad1 = 1;
nmad2 = min(nobs-npmad+1,fix((time_obs_20min(n)-time_obs_20min(1))*1440./20.) + 1);
minmad1 = max(1,nmad1);
maxmad2 = min(nobs-midm1-1,nmad2);
kob = 0;
nob = 0;
jvi = 0;
jti = 0;
jte = 0;
jne = 0;
for i=minmad1:maxmad2
 icen = i+midm1;
 [obswimed,obswiresid,obswimad] = get_mad(wi_obs_nan_20min,npmad,i,midm1);
 [obstimed,obstiresid,obstimad] = get_mad(ti_obs_nan_20min,npmad,i,midm1);
 [obstemed,obsteresid,obstemad] = get_mad(te_obs_nan_20min,npmad,i,midm1);
 [obsnemed,obsneresid,obsnemad] = get_mad(ne_obs_nan_20min,npmad,i,midm1);
 if ~isnan(obswimed)
   jvi = jvi + 1;
   timemadwi_obs(jvi) = time_obs_20min(icen);
   wimed_obs(jvi) = obswimed;
   wiresid_obs(jvi) = obswiresid;
   wimad_obs(jvi) = obswimad;
 end
 if ~isnan(obstimed)
   jti = jti + 1;
   timemadti_obs(jti) = time_obs_20min(icen);
   timed_obs(jti) = obstimed;
   tiresid_obs(jti) = obstiresid;
   timad_obs(jti) = obstimad;
 end
 if ~isnan(obstemed)
   jte = jte + 1;
   timemadte_obs(jte) = time_obs_20min(icen);
   temed_obs(jte) = obstemed;
   teresid_obs(jte) = obsteresid;
   temad_obs(jte) = obstemad;
 end
 if ~isnan(obsnemed)
   jne = jne + 1;
   timemadne_obs(jne) = time_obs_20min(icen);
   nemed_obs(jne) = obsnemed;
   neresid_obs(jne) = obsneresid;
   nemad_obs(jne) = obsnemad;
 end
 [obsnmaxmed,obsnmaxresid,obsnmaxmad] = get_mad(nmax_obs_nan_20min,npmad,i,midm1);
 [obshmaxmed,obshmaxresid,obshmaxmad] = get_mad(hmax_obs_nan_20min,npmad,i,midm1);
 if ~isnan(obshmaxmed)
   kob = kob + 1;
   timemadn_obs(kob) = time_obs_20min(icen);
   hmaxmed_obs(kob) = obshmaxmed;
   hmaxresid_obs(kob) = obshmaxresid;
   hmaxmad_obs(kob) = obshmaxmad;
   nmaxmed_obs(kob) = obsnmaxmed;
   nmaxresid_obs(kob) = obsnmaxresid;
   nmaxmad_obs(kob) = obsnmaxmad;
 end
 [obstecmed,obstecresid,obstecmad] = get_mad(tec_obs_nan_20min,npmad,i,midm1);
 if ~isnan(obstecmed)
   nob = nob + 1;
   timemadt_obs(nob) = time_obs_20min(icen);
   tecmed_obs(nob) = obstecmed;
   tecresid_obs(nob) = obstecresid;
   tecmad_obs(nob) = obstecmad;
 end
end
jvi,jti,jte,jne
% Calculate total mads:
mad_wiobs = 0.;
mad_tiobs = 0.;
mad_teobs = 0.;
mad_neobs = 0.;
jmad = 0;
if (jvi > 0)
 for i=1:jvi
  if ~isnan(wimad_obs(i))
   jmad = jmad + 1;
   mad_wiobs = mad_wiobs + wimad_obs(i);
  end
 end
 mad_wiobs = mad_wiobs/jmad;
end
nmad_wiobs = jmad;
jmad = 0;
if (jti > 0)
 for i=1:jti
  if ~isnan(timad_obs(i))
   jmad = jmad + 1;
   mad_tiobs = mad_tiobs + timad_obs(i);
  end
 end
 mad_tiobs = mad_tiobs/jmad;
end
nmad_tiobs = jmad;
jmad = 0;
if (jte > 0)
 for i=1:jte
  if ~isnan(temad_obs(i))
   jmad = jmad + 1;
   mad_teobs = mad_teobs + temad_obs(i);
  end
 end
 mad_teobs = mad_teobs/jmad;
end
nmad_teobs = jmad;
jmad = 0;
if (jne > 0)
 for i=1:jne
  if ~isnan(nemad_obs(i))
   jmad = jmad + 1;
   mad_neobs = mad_neobs + nemad_obs(i);
  end
 end
 mad_neobs = mad_neobs/jmad;
end
nmad_neobs = jmad;
mad_nmaxobs = 0.;
mad_hmaxobs = 0.;
if (kob > 0)
kmad = 0;
 for i=1:kob
  if ~isnan(nmaxmad_obs(i))
   kmad = kmad + 1;
   mad_hmaxobs = mad_hmaxobs + hmaxmad_obs(i);
   mad_nmaxobs = mad_nmaxobs + nmaxmad_obs(i);
  end
 end
   mad_hmaxobs = mad_hmaxobs/kmad;
   mad_nmaxobs = mad_nmaxobs/kmad;
end
mad_tecobs = 0.;
if (nob > 0)
nmad = 0;
 for i=1:nob
  if ~isnan(tecmad_obs(i))
   nmad = nmad + 1;
   mad_tecobs = mad_tecobs + tecmad_obs(i);
  end
 end
 mad_tecobs = mad_tecobs/nmad;
end
 
 
if (notie == 0)


fpath = '/hao/ftp/pub/emery/';

%%  Read TIEGCM .nc files starting at 0020 UT every 20 min (so itwd = nobs-1)


itwd = 0;
it1 = 1;
if (year1 < 2005.5)
    Fname=[fpath '/wd194/wd.s05135_001.nc'];
else
    Fname=[fpath '/wd194/wd.s06348_349_001.nc'];
end
    [nt,dayfrac,vpeav,vpnav,uiav,viav,wiav,tnav,tiav,teav,neav,nmax,hmax,tec,nz,ht,htne,htte,htti,htmid,httn,htui,htvi,htwi] = get_netein_nc(Fname,ht1,ht2,ii,jj,deginc,degdec);
 for n=1:nz
  htav_tie(n) = 0.;
 end
for it=it1:nt
 nmax_wd(itwd+it-it1+1) = log10(nmax(it));
 hmax_wd(itwd+it-it1+1) = hmax(it);
 tec_wd(itwd+it-it1+1) = log10(tec(it));
 tn_wd(itwd+it-it1+1) = tnav(it);
 timewd(itwd+it-it1+1) = dayfrac(it);
 vpe_wd(itwd+it-it1+1) = vpeav(it);
 vpn_wd(itwd+it-it1+1) = vpnav(it);
 ui_wd(itwd+it-it1+1) = uiav(it);
 vi_wd(itwd+it-it1+1) = viav(it);
 wi_wd(itwd+it-it1+1) = wiav(it);
 ti_wd(itwd+it-it1+1) = tiav(it);
 te_wd(itwd+it-it1+1) = teav(it);
 ne_wd(itwd+it-it1+1) = log10(neav(it));
 for n=1:nz
  ht_timewd(n,itwd+it-it1+1) = timewd(itwd+it-it1+1);
  htav_tie(n) = htav_tie(n) + ht(n,it);
  ht_tie(n,itwd+it-it1+1) = ht(n,it);
  ne_tie(n,itwd+it-it1+1) = log10(htne(n,it));
  if (n < nz)
   ui_tie(n,itwd+it-it1+1) = htui(n,it);
   vi_tie(n,itwd+it-it1+1) = htvi(n,it);
   wi_tie(n,itwd+it-it1+1) = htwi(n,it);
   te_tie(n,itwd+it-it1+1) = htte(n,it);
   ti_tie(n,itwd+it-it1+1) = htti(n,it);
   htmid_time(n,itwd+it-it1+1) = timewd(itwd+it-it1+1);
   ht_mid(n,itwd+it-it1+1) = htmid(n,it);
   tn_tie(n,itwd+it-it1+1) = httn(n,it);
  end
 end
end
nt1 = nt;
itwd = itwd + nt - (it1-1);
if (year1 > 2005.5)
    Fname=[fpath '/wd194/wd.s06348_349_002.nc'];
    [nt,dayfrac,vpeav,vpnav,uiav,viav,wiav,tnav,tiav,teav,neav,nmax,hmax,tec,nz,ht,htne,htte,htti,htmid,httn,htui,htvi,htwi] = get_netein_nc(Fname,ht1,ht2,ii,jj,deginc,degdec);
for it=it1:nt
 nmax_wd(itwd+it-it1+1) = log10(nmax(it));
 hmax_wd(itwd+it-it1+1) = hmax(it);
 tec_wd(itwd+it-it1+1) = log10(tec(it));
 tn_wd(itwd+it-it1+1) = tnav(it);
 timewd(itwd+it-it1+1) = dayfrac(it);
 vpe_wd(itwd+it-it1+1) = vpeav(it);
 vpn_wd(itwd+it-it1+1) = vpnav(it);
 ui_wd(itwd+it-it1+1) = uiav(it);
 vi_wd(itwd+it-it1+1) = viav(it);
 wi_wd(itwd+it-it1+1) = wiav(it);
 ti_wd(itwd+it-it1+1) = tiav(it);
 te_wd(itwd+it-it1+1) = teav(it);
 ne_wd(itwd+it-it1+1) = log10(neav(it));
 for n=1:nz
  ht_timewd(n,itwd+it-it1+1) = timewd(itwd+it-it1+1);
  htav_tie(n) = htav_tie(n) + ht(n,it);
  ht_tie(n,itwd+it-it1+1) = ht(n,it);
  ne_tie(n,itwd+it-it1+1) = log10(htne(n,it));
  if (n < nz)
   ui_tie(n,itwd+it-it1+1) = htui(n,it);
   vi_tie(n,itwd+it-it1+1) = htvi(n,it);
   wi_tie(n,itwd+it-it1+1) = htwi(n,it);
   te_tie(n,itwd+it-it1+1) = htte(n,it);
   ti_tie(n,itwd+it-it1+1) = htti(n,it);
   htmid_time(n,itwd+it-it1+1) = timewd(itwd+it-it1+1);
   ht_mid(n,itwd+it-it1+1) = htmid(n,it);
   tn_tie(n,itwd+it-it1+1) = httn(n,it);
  end
 end
end
for n=1:nz
  htav_tie(n) = htav_tie(n)/(nt1+nt);
end
itwd = itwd + nt - (it1-1);
end	% of year1>2005.5

%% calculate rms and print from obs every 20 min 


RMS = fopen(filerms,'w');
fprintf(RMS, 'RMS of 20-min TIEGCM values of hmF2, h=NmF2, TEC from avs every 20 min of Sondrestrom ISR obs from %6.2f km +/-20km \n',ht300);
fprintf(RMS, '%d  nobs    hmF2     NmF2      TEC \n',year1);
iend = min(itwd,nobs);
 rmshmaxwd = 0.;
 [rmshmaxwd,nhmaxrms] = get_rms(hmax_obs_nan_20min,hmax_wd,1,iend);
 rmsnmaxwd = 0.;
 [rmsnmaxwd,nnmaxrms] = get_rms(nmax_obs_nan_20min,nmax_wd,1,iend);
 rmstecwd = 0.;
 [rmstecwd,ntecrms] = get_rms(tec_obs_nan_20min,tec_wd,1,iend);
 rmsvpe300wd = 0.;
 [rmsvpe300wd,nvpe300rms] = get_rms(vpe_obs_nan_20min,vpe_wd,1,iend);
 rmsvpn300wd = 0.;
 [rmsvpn300wd,nvpn300rms] = get_rms(vpn_obs_nan_20min,vpn_wd,1,iend);
 rmswi300wd = 0.;
 [rmswi300wd,nwi300rms] = get_rms(wi_obs_nan_20min,wi_wd,1,iend);
 rmsti300wd = 0.;
 [rmsti300wd,nti300rms] = get_rms(ti_obs_nan_20min,ti_wd,1,iend);
 rmste300wd = 0.;
 [rmste300wd,nte300rms] = get_rms(te_obs_nan_20min,te_wd,1,iend);
 rmsne300wd = 0.;
 [rmsne300wd,nne300rms] = get_rms(ne_obs_nan_20min,ne_wd,1,iend);
fprintf(RMS, '%d-%d %4.0d %13.4e %3.0d %13.4e %3.0d %13.4e %3.0d %13.4e %3.0d %7.2f %3.0d %7.2f %3.0d %7.2f %3.0d \n \n', day1,day2,nobs,rmshmaxwd,nhmaxrms,rmsnmaxwd,nnmaxrms,rmstecwd,ntecrms, rmsne300wd,nne300rms,rmste300wd,nte300rms,rmsti300wd,nti300rms,rmswi300wd,nwi300rms);

% Calculate the median absolute deviation (mad) sliding each 20 min for intervals of npmad pts
%  Best if npmad is odd number.  Usually is over an integer number of hours numhr.  
%  The median, the residual, and the mad are saved for the midpt of npmad pts or at (npmad+1)/2
nmad = itwd-npmad+1
mad_nmax_wd = 0.;
mad_hmax_wd = 0.;
mad_tec_wd = 0.;
mad_wi_wd = 0.;
mad_tn_wd = 0.;
mad_ti_wd = 0.;
mad_te_wd = 0.;
mad_ne_wd = 0.;
for i=1:nmad
 icen = i+midm1;
 timemad(icen-midm1) = time_obs_20min(icen);
 [hmaxmed_wd(icen-midm1),hmaxresid_wd(icen-midm1),hmaxmad_wd(icen-midm1)] = get_mad(hmax_wd,npmad,i,midm1);
 [nmaxmed_wd(icen-midm1),nmaxresid_wd(icen-midm1),nmaxmad_wd(icen-midm1)] = get_mad(nmax_wd,npmad,i,midm1);
 [tecmed_wd(icen-midm1),tecresid_wd(icen-midm1),tecmad_wd(icen-midm1)] = get_mad(tec_wd,npmad,i,midm1);
 [tnmed_wd(icen-midm1),tnresid_wd(icen-midm1),tnmad_wd(icen-midm1)] = get_mad(tn_wd,npmad,i,midm1);
 [wimed_wd(icen-midm1),wiresid_wd(icen-midm1),wimad_wd(icen-midm1)] = get_mad(wi_wd,npmad,i,midm1);
 [timed_wd(icen-midm1),tiresid_wd(icen-midm1),timad_wd(icen-midm1)] = get_mad(ti_wd,npmad,i,midm1);
 [temed_wd(icen-midm1),teresid_wd(icen-midm1),temad_wd(icen-midm1)] = get_mad(te_wd,npmad,i,midm1);
 [nemed_wd(icen-midm1),neresid_wd(icen-midm1),nemad_wd(icen-midm1)] = get_mad(ne_wd,npmad,i,midm1);
% Calculate total mads:
 mad_hmax_wd = mad_hmax_wd + hmaxmad_wd(i);
 mad_nmax_wd = mad_nmax_wd + nmaxmad_wd(i);
 mad_tec_wd = mad_tec_wd + tecmad_wd(i);
 mad_wi_wd = mad_wi_wd + wimad_wd(i);
 mad_tn_wd = mad_tn_wd + tnmad_wd(i);
 mad_ti_wd = mad_ti_wd + timad_wd(i);
 mad_te_wd = mad_te_wd + temad_wd(i);
 mad_ne_wd = mad_ne_wd + nemad_wd(i);
end
mad_hmax_wd = mad_hmax_wd/nmad;
mad_nmax_wd = mad_nmax_wd/nmad;
mad_tec_wd = mad_tec_wd/nmad;
mad_tn_wd = mad_tn_wd/nmad;
mad_ti_wd = mad_ti_wd/nmad;
mad_te_wd = mad_te_wd/nmad;
mad_ne_wd = mad_ne_wd/nmad;
fprintf(RMS, '\n MAD (median absolute deviations) of 20-min Sondrestrom ISR obs and TIEGCM values \n');
fprintf(RMS, '%d     obs   wd1.94 hmax nmax tec at %6.2f km: ne te ti tn wi \n',year1,ht300);
fprintf(RMS, '%d-%d %4.0d %13.4e %13.4e %13.4e %13.4e %7.2f %7.2f %7.2f %7.2f \n \n', day1,day2,nmad,mad_hmax_wd,mad_nmax_wd,mad_tec_wd,mad_ne_wd,mad_te_wd,mad_ti_wd,mad_tn_wd,mad_wi_wd);

end	% Skip reading TIEGCM if notie=1


%% Make plots

% position [left bottom width height]
fsize = 12;
f9 = 10;
wr1c = 0.33;
wr1m = 0.06;
wr1v = 0.01;
pcwr = 0.10;

yi1 = -40;
yi2 = 40;
yi95 = yi1 + 0.92*(yi2-yi1);
yitot = [yi1,yi2];

 y1=501;
 y2 = 3499;
y95 = y1 + 0.92*(y2-y1);
y90 = y1 + 0.85*(y2-y1);
y85 = y1 + 0.78*(y2-y1);

yw1 = -99;
yw2 = 299;
% Add VpE,N to Wi
yw1 = -1599.;
yw2 = 1599.;
yw95 = yw1 + 0.92*(yw2-yw1);
yw90 = yw1 + 0.75*(yw2-yw1);

yh1 = 101;
yh2 = 599;

ye1 = .9e+9;
ye1 = 0.1;
ye2 = .39e+12;
% Change to log10
%ye1 = 9.01;
%ye2 = 12.99;
ye95 = ye1 + 0.92*(ye2-ye1);

yt1 = 0.;
yt2 = .99e+17;
% Change to log10
%yt1 = 15.01;
%yt2 = 17.99;
yt95 = yt1 + 0.92*(yt2-yt1);

alt1 = 100;
alt2 = 699;
altcb = 0.90*(alt2-alt1)+alt1;


 x1 = day1;
 x2 = day2;
if (year1 < 2005.5)
 if (iflp==0)
 daytit = sprintf('Sondrestrom SP (66.99N, 50.95W) 15 May 2005');
 plottit = sprintf('son_sp05135');
 plotht = sprintf('sonht_sp05135');
 plottie = sprintf('sontie_sp05135');
 plotmadt = sprintf('sonmadt_sp05135');
 plotmadn = sprintf('sonmadn_sp05135');
 if (notie==1)
  plottit = sprintf('son_data_sp05135');
 end
 else
 daytit = sprintf('Sondrestrom LP Z (66.99N, 50.95W) 15 May 2005');
 plottit = sprintf('son_lp05135');
 plotht = sprintf('sonht_lp05135');
 plothtv = sprintf('sonhtv_lp05135');
 plottie = sprintf('sontie_lp05135');
 plotmadt = sprintf('sonmadt_lp05135');
 plotmadn = sprintf('sonmadn_lp05135');
 if (notie==1)
  plottit = sprintf('son_data_lp05135');
 end
%daytit = sprintf('Sondrestrom LP elev=75 (66.99N, 50.95W) 15 May 2005');
%plottit = sprintf('son_lp75_05135');
%plotht = sprintf('sonht_lp75_05135');
%plothtv = sprintf('sonhtv_lp75_05135');
%plottie = sprintf('sontie_lp75_05135');
%plotmadt = sprintf('sonmadt_lp75_05135');
%plotmadn = sprintf('sonmadn_lp75_05135');
%if (notie==1)
% plottit = sprintf('son_data_lp75_05135');
%end
% daytit = sprintf('Sondrestrom LP 75d+Z (66.99N, 50.95W) 15 May 2005');
% plottit = sprintf('son_lpc_05135');
% plotht = sprintf('sonht_lpc_05135');
% plothtv = sprintf('sonhtv_lpc_05135');
% plottie = sprintf('sontie_lpc_05135');
% plotmadt = sprintf('sonmadt_lpc_05135');
% plotmadn = sprintf('sonmadn_lpc_05135');
% if (notie==1)
%  plottit = sprintf('son_data_lpc_05135');
% end
 end
else
 daytit = sprintf('Sondrestrom (66.99N, 50.95W) 14-15 Dec 2006');
 plottit = sprintf('son_06348_9');
 plotht = sprintf('sonht_06348_9');
 plothtv = sprintf('sonhtv_06348_9');
 plottie = sprintf('sontie_06348_9');
 plotmadt = sprintf('sonmadt_06348_9');
 plotmadn = sprintf('sonmadn_06348_9');
 if (notie==1)
  plottit = sprintf('son_data_06348_9');
 end
end

xtot = [x1,x2];

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

subplot('position',[0.08,0.750,0.87,0.210])	% 1/4-5
hold on;
title(daytit,'FontSize',fsize);
plot (dayfracimf,bximf,'k');
plot (dayfracimf,byimf,'g');
plot (dayfracimf,bzimf,'r');
plot (dayfracimf,vsw/100.,'c');
plot(xtot,zero,':k');
plot (begut,yitot,':k');
plot (endut,yitot,':k');
text (x1+0.01*(x2-x1),yi95,'Bx','Color','k','FontSize',f9);
text (x1+0.05*(x2-x1),yi95,'By','Color','g','FontSize',f9);
text (x1+0.09*(x2-x1),yi95,'Bz','Color','r','FontSize',f9);
text (x1+0.13*(x2-x1),yi95,'Vsw/100','Color','c','FontSize',f9);
ylabel('IMF (nT), Vsw/100 (km/s)', 'FontSize',fsize);
xlabel(' ','FontSize',fsize);
set (gca,'xlim',[x1,x2],'ylim',[yi1,yi2],'Xtick',([x1:0.25:x2]),'XTickLabel',{' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';},'FontSize',fsize,'box','on');

% Obs in order Tn, Ti, Te for colors/lines: k c r
subplot('position',[0.08,0.540,0.87,0.210])	% 2/4-5
hold on;
if (if1m==1)
 plot (timeobs,tiobs300,'Color','c');
 plot (timeobs,teobs300,'Color','r');
else
 plot (time_obs_20min,ti_obs_nan_20min,'Color','c','LineStyle','-');
 plot (time_obs_20min,te_obs_nan_20min,'Color','r','LineStyle','-');
end
 text (x1+0.01*(x2-x1),y95,'Tn','Color','k','FontSize',f9);
 text (x1+0.10*(x2-x1),y95,'Ti','Color','c','FontSize',f9);
 text (x1+0.19*(x2-x1),y95,'Te','Color','r','FontSize',f9);
if (notie ~=1)
text (x1+0.28*(x2-x1),y90,'RMS','Color','k','FontSize',f9);
rms = sprintf('%4.2f',rmste300wd);
text (x1+0.19*(x2-x1),y90,rms,'Color','r','FontSize',f9);
rms = sprintf('%4.2f',rmsti300wd);
text (x1+0.10*(x2-x1),y90,rms,'Color','c','FontSize',f9);
text (x1+0.28*(x2-x1),y95,'--- av TIEGCM ','Color','k','FontSize',f9);
text (x1+0.80*(x2-x1),y95,ht300str,'Color','k','FontSize',f9);
plot (timewd,tn_wd,'Color','k','LineStyle','--');
plot (timewd,ti_wd,'Color','c','LineStyle','--');
plot (timewd,te_wd,'Color','r','LineStyle','--');
end
 ylabel('T (K)','FontSize',fsize);
xlabel('Daynumber','FontSize',fsize);
set (gca,'xlim',[x1,x2],'ylim',[y1,y2],'Xtick',([x1:0.25:x2]),'XTickLabel',{' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';},'FontSize',fsize,'box','on');

subplot('position',[0.08,0.430,0.87,0.110])	% 3/6
hold on;
if (if1m==1)
 plot (timeobs,wiobs300,'Color','r');
else
 plot (time_obs_20min,wi_obs_nan_20min*10.,'Color','r');
 plot (time_obs_20min,vpe_obs_nan_20min,'Color','g');
 plot (time_obs_20min,vpn_obs_nan_20min,'Color','b');
 text (x1+0.05*(x2-x1),yw95,'Wi (x10)','Color','r','FontSize',f9);
 text (x1+0.16*(x2-x1),yw95,'VpE','Color','g','FontSize',f9);
 text (x1+0.26*(x2-x1),yw95,'VpN','Color','b','FontSize',f9);
end
plot(xtot,zero,':k');
if (notie ~= 1)
 rms = sprintf('RMS %4.2f',rmswi300wd);
 text (x1+0.01*(x2-x1),yw90,rms,'Color','r','FontSize',f9);
 plot (timewd,wi_wd*10.,'Color','r','LineStyle','--');
 rms = sprintf('%4.2f',rmsvpe300wd);
 text (x1+0.16*(x2-x1),yw90,rms,'Color','g','FontSize',f9);
 plot (timewd,vpe_wd,'Color','g','LineStyle','--');
 rms = sprintf('%4.2f',rmsvpn300wd);
 text (x1+0.26*(x2-x1),yw90,rms,'Color','b','FontSize',f9);
 plot (timewd,vpn_wd,'Color','b','LineStyle','--');
end
ylabel('Vi (m/s)','FontSize',fsize);
xlabel('Daynumber','FontSize',fsize);
set (gca,'xlim',[x1,x2],'ylim',[yw1,yw2],'Xtick',([x1:0.25:x2]),'XTickLabel',{' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';},'FontSize',fsize,'box','on');


subplot('position',[0.08,0.330,0.87,0.100])	% 4/6
hold on;
if (if1m==1)
 plot (timeobs,hmaxobs,'Color','r');
else
 plot (time_obs_20min,hmax_obs_nan_20min,'Color','r');
end
plot(xtot,z300,':k');
if (notie ~= 1)
 plot (timewd,hmax_wd,'Color','r','LineStyle','--');
end
ylabel('hmax (km)','FontSize',fsize);
xlabel('Daynumber','FontSize',fsize);
set (gca,'xlim',[x1,x2],'ylim',[yh1,yh2],'Xtick',([x1:0.25:x2]),'XTickLabel',{' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';},'FontSize',fsize,'box','on');

subplot('position',[0.08,0.190,0.87,0.140])	% 5/6
hold on;
if (if1m==1)
 plot (timeobs,nmaxobs,'Color','r');
 plot (timeobs,neobs300,'Color','k');
else
 plot (time_obs_20min,10.^ne_obs_nan_20min,'Color','k');
 plot (time_obs_20min,10.^nmax_obs_nan_20min,'Color','r');
%plot (time_obs_20min,ne_obs_nan_20min,'Color','k');
%plot (time_obs_20min,nmax_obs_nan_20min,'Color','r');
end
text (x1+0.01*(x2-x1),ye95,'Nmax','Color','r','FontSize',f9);
text (x1+0.10*(x2-x1),ye95,ht300str,'Color','k','FontSize',f9);
text (x1+0.55*(x2-x1),ye95,'Contact Anja Stromme for ISR data usage','Color','k','FontSize',f9);
if (notie ~= 1)
% plot (timewd,ne_wd,'Color','k','LineStyle','--');
% plot (timewd,nmax_wd,'r','LineStyle','--');
 plot (timewd,10.^ne_wd,'Color','k','LineStyle','--');
 plot (timewd,10.^nmax_wd,'r','LineStyle','--');
 text (x1+0.19*(x2-x1),ye95,'-- TIEGCM','Color','k','FontSize',f9);
end
%ylabel('log10(Ne m-3)','FontSize',fsize);
ylabel('Ne m-3','FontSize',fsize);
xlabel('Daynumber','FontSize',fsize);
set (gca,'xlim',[x1,x2],'ylim',[ye1,ye2],'Xtick',([x1:0.25:x2]),'XTickLabel',{' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';},'FontSize',fsize,'box','on');

% Markers: +,o,.,x,*,s,d,p,h,^,v,<,>
subplot('position',[0.08,0.050,0.87,0.140])	% 6/6
hold on;
if (if1m==1)
 plot (timeobs,tecobs,'Color','r');
else
% plot (time_obs_20min,tec_obs_nan_20min,'Color','r');
 plot (time_obs_20min,10.^tec_obs_nan_20min,'Color','r');
end
if (notie ~= 1)
text (x1+0.19*(x2-x1),yt95,'-- TIEGCM','Color','k','FontSize',f9);
% plot (timewd,tec_wd,'r','LineStyle','--');
 plot (timewd,10.^tec_wd,'r','LineStyle','--');
end
%ylabel('Log10(TEC m-2)','FontSize',fsize);
ylabel('TEC (m-2)','FontSize',fsize);
xlabel('Daynumber','FontSize',fsize);
set (gca,'xlim',[x1,x2],'Xtick',([x1:0.25:x2]),'ylim',[yt1,yt2],'FontSize',fsize,'box','on');
hold off;

%set(gcf,'color',[1 1 1],'inverthardcopy','off');
%set(gcf,'PaperPosition',[0.5, 0.5, 7.5, 10.5]);
set(gcf,'PaperPosition',[0.5, 0.5, 7.5, 10.5]);
set(gcf,'PaperPositionmode', 'auto');
%saveas(gcf,plottit,'tif')
saveas(gcf,plottit,'epsc')
%saveas(gcf,plottit,'jpg')


%% Plot contours of Ne, Te, Ti, Tn for obs (and std) and TIEGCM
% Define contour levels for Ne, Te, Ti-Tn
%wilev = [-300 -250 -200 -150 -100 -50 0 50 100 150 200 250 300 350 400 450 500];
wilev = [-100 -50 0 50 100 150 200 250 300 350 400];
tnlev = [200 350 500 650 800 950 1100 1250 1400 1550 1700 1850];
telev = [250 500 750 1000 1250 1500 1750 2000 2250 2500 2750 3000];
nelev = [9.0 9.33 9.67 10.0 10.33 10.67 11.0 11.33 11.67 12.0];

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

h1=subplot('position',[0.08,0.735,0.87,0.225]);	% 1/4
hold on;
title(daytit,'FontSize',fsize);
% min/max Dec06 6.6-13.7; 10-3000, 10-2079, 184-1011
if (if1m==1)
 contourf(timeobs,altobs,neobs,nelev);
 plot (timeobs,hmaxobs,'Color','r');
else
 contourf(time_obs_20min,altobs,neobs_20min,nelev);
 plot (time_obs_20min,hmax_obs_nan_20min,'Color','r');
end
[min1,max1] = caxis;
min1,max1
caxis([9.,12.]);
xlabel(' ','FontSize',fsize);
set (gca,'xlim',[x1,x2],'ylim',[alt1,alt2],'Xtick',([x1:0.25:x2]),'XTickLabel',{' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';},'FontSize',fsize,'box','on');
text(0.01*(x2-x1)+x2, altcb, 'lg(Ne) m-3','FontSize',fsize);
pos1=get(h1,'position');
set(h1,'position',pos1);
hp=colorbar;
pos=get(hp,'position');
pos(4)=0.19;
%pos(3)=0.5*pos(3)	% narrows the color-bar by factor of 2
set(hp,'position',pos);
pos1(3)=0.80;
set(gca,'position',pos1);


h1=subplot('position',[0.08,0.510,0.87,0.225]);	% 2/4
hold on;
if (if1m==1)
 contourf(timeobs,altobs,teobs,telev);
else
 contourf(time_obs_20min,altobs,teobs_20min,telev);
end
[min1,max1] = caxis;
min1,max1
caxis([250.,3000.]);
plot(xtot,z300,':k');
xlabel(' ','FontSize',fsize);
set (gca,'xlim',[x1,x2],'ylim',[alt1,alt2],'Xtick',([x1:0.25:x2]),'XTickLabel',{' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';},'FontSize',fsize,'box','on');
text (0.01*(x2-x1)+x2, altcb, 'Te(K)','FontSize',fsize);
pos1=get(h1,'position');
set(h1,'position',pos1)
hp=colorbar;
pos=get(hp,'position');
pos(4) = 0.19;
set(hp,'position',pos);
pos1(3)=0.80;
set(gca,'position',pos1);

h1=subplot('position',[0.08,0.285,0.87,0.225]);	% 3/4
hold on;
if (if1m==1)
 contourf(timeobs,altobs,tiobs,tnlev);
else
 contourf(time_obs_20min,altobs,tiobs_20min,tnlev);
end
[min1,max1] = caxis;
min1,max1
caxis([200.,1850.]);
plot(xtot,z300,':k');
xlabel(' ','FontSize',fsize);
set (gca,'xlim',[x1,x2],'ylim',[alt1,alt2],'Xtick',([x1:0.25:x2]),'XTickLabel',{' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';},'FontSize',fsize,'box','on');
text (0.01*(x2-x1)+x2, altcb, 'Ti(K)','FontSize',fsize);
pos1=get(h1,'position');
set(h1,'position',pos1)
hp=colorbar;
pos=get(hp,'position');
pos(4) = 0.19;
set(hp,'position',pos);
pos1(3)=0.80;
set(gca,'position',pos1);

h1=subplot('position',[0.08,0.06,0.87,0.225]);	% 4/4
hold on;
if (if1m==1)
 contourf(timeobs,altobs,wiobs,wilev);
else
 contourf(time_obs_20min,altobs,wiobs_20min,wilev);
end
[min1,max1] = caxis;
min1,max1
caxis([-100.,400.]);
plot(xtot,z300,':k');
xlabel(' ','FontSize',fsize);
set (gca,'xlim',[x1,x2],'ylim',[alt1,alt2],'Xtick',([x1:0.25:x2]),'FontSize',fsize,'box','on');
text(0.01*(x2-x1)+x2, altcb, 'Wi m/s','FontSize',fsize);
pos1=get(h1,'position');
set(h1,'position',pos1)
hp=colorbar;
pos=get(hp,'position');
pos(4) = 0.19;
set(hp,'position',pos);
pos1(3)=0.80;
set(gca,'position',pos1);

hold off;

set(gcf,'PaperPosition',[0.5, 0.5, 7.5, 10.5]);
set(gcf,'PaperPositionmode', 'auto');
saveas(gcf,plotht,'epsc')

% Plot the variance in the 20min medians
if (if1m~=1)
% Define contour levels for std Te, Ti
stdtei = [0 100 200 300 400 500 600 700 800 900 1000];
stdwi = [0 50 100 150 200 250 300 350 400 450 500];

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

h1=subplot('position',[0.08,0.735,0.87,0.225]);	% 1/4
hold on;
title(daytit,'FontSize',fsize);
% min/max Dec06 6.6-13.7; 10-3000, 10-2079, 184-1011
 contourf(time_obs_20min,altobs,neobsv_20min,nelev);
 plot (time_obs_20min,hmax_obs_nan_20min,'Color','r');
[min1,max1] = caxis;
min1,max1
caxis([9.,12.]);
xlabel(' ','FontSize',fsize);
set (gca,'xlim',[x1,x2],'ylim',[alt1,alt2],'Xtick',([x1:0.25:x2]),'XTickLabel',{' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';},'FontSize',fsize,'box','on');
text(0.01*(x2-x1)+x2, altcb, 'lgstdNe m-3','FontSize',fsize);
pos1=get(h1,'position');
set(h1,'position',pos1);
hp=colorbar;
pos=get(hp,'position');
pos(4)=0.19;
%pos(3)=0.5*pos(3)	% narrows the color-bar by factor of 2
set(hp,'position',pos);
pos1(3)=0.80;
set(gca,'position',pos1);


h1=subplot('position',[0.08,0.510,0.87,0.225]);	% 2/4
hold on;
 contourf(time_obs_20min,altobs,teobsv_20min,stdtei);
[min1,max1] = caxis;
min1,max1
caxis([0.,1000.]);
plot(xtot,z300,':k');
xlabel(' ','FontSize',fsize);
set (gca,'xlim',[x1,x2],'ylim',[alt1,alt2],'Xtick',([x1:0.25:x2]),'XTickLabel',{' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';},'FontSize',fsize,'box','on');
text (0.01*(x2-x1)+x2, altcb, 'stdTe(K)','FontSize',fsize);
pos1=get(h1,'position');
set(h1,'position',pos1)
hp=colorbar;
pos=get(hp,'position');
pos(4) = 0.19;
set(hp,'position',pos);
pos1(3)=0.80;
set(gca,'position',pos1);

h1=subplot('position',[0.08,0.285,0.87,0.225]);	% 3/4
hold on;
 contourf(time_obs_20min,altobs,tiobsv_20min,stdtei);
[min1,max1] = caxis;
min1,max1
caxis([0.,1000.]);
plot(xtot,z300,':k');
xlabel(' ','FontSize',fsize);
set (gca,'xlim',[x1,x2],'ylim',[alt1,alt2],'Xtick',([x1:0.25:x2]),'XTickLabel',{' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';},'FontSize',fsize,'box','on');
text (0.01*(x2-x1)+x2, altcb, 'stdTi(K)','FontSize',fsize);
pos1=get(h1,'position');
set(h1,'position',pos1)
hp=colorbar;
pos=get(hp,'position');
pos(4) = 0.19;
set(hp,'position',pos);
pos1(3)=0.80;
set(gca,'position',pos1);

h1=subplot('position',[0.08,0.06,0.87,0.225]);	% 4/4
hold on;
 contourf(time_obs_20min,altobs,wiobsv_20min,stdwi);
[min1,max1] = caxis;
min1,max1
caxis([0.,300.]);
plot(xtot,z300,':k');
xlabel(' ','FontSize',fsize);
set (gca,'xlim',[x1,x2],'ylim',[alt1,alt2],'Xtick',([x1:0.25:x2]),'FontSize',fsize,'box','on');
text(0.01*(x2-x1)+x2, altcb, 'stdVi m/s','FontSize',fsize);
pos1=get(h1,'position');
set(h1,'position',pos1)
hp=colorbar;
pos=get(hp,'position');
pos(4) = 0.19;
set(hp,'position',pos);
pos1(3)=0.80;
set(gca,'position',pos1);

hold off;

set(gcf,'PaperPosition',[0.5, 0.5, 7.5, 10.5]);
set(gcf,'PaperPositionmode', 'auto');
saveas(gcf,plothtv,'epsc')
end	% if if1m~=1

if (notie ~= 1)
%% Plot contours of Ne, Te, Ti, and Tn from TIEGCM
figure;
% Set size to be 8.5x11 inches
set(gcf,'position',[142    99   647   822])
H = gcf;

h1=subplot('position',[0.08,0.735,0.87,0.225]);	% 1/4
hold on;
title(daytit,'FontSize',fsize);
% min/max Dec06 6.6-13.7; 10-3000, 10-2079, 184-1011
%contourf(timewd,htav_tie,ne_tie,9);
contourf(ht_timewd,ht_tie,ne_tie,nelev);
[min1,max1] = caxis;
min1,max1
caxis([9.,12.]);
 plot (timewd,hmax_wd,'Color','r');
xlabel(' ','FontSize',fsize);
set (gca,'xlim',[x1,x2],'ylim',[alt1,alt2],'Xtick',([x1:0.25:x2]),'XTickLabel',{' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';},'FontSize',fsize,'box','on');
text(0.01*(x2-x1)+x2, altcb, 'lg(Ne) m-3','FontSize',fsize);
text (0.01*(x2-x1)+x2,0.05*(alt2-alt1)+alt2,'TIEGCM','FontSize',fsize)
pos1=get(h1,'position');
set(h1,'position',pos1);
hp=colorbar;
pos=get(hp,'position');
pos(4)=0.19;
%pos(3)=0.5*pos(3)	% narrows the color-bar by factor of 2
set(hp,'position',pos);
pos1(3)=0.80;
set(gca,'position',pos1);


h1=subplot('position',[0.08,0.510,0.87,0.225]);	% 2/4
hold on;
%contourf(ht_timewd,ht_tie,te_tie,13);
contourf(htmid_time,ht_mid,te_tie,telev);
[min1,max1] = caxis;
min1,max1
caxis([250.,3000.]);
plot(xtot,z300,':k');
xlabel(' ','FontSize',fsize);
set (gca,'xlim',[x1,x2],'ylim',[alt1,alt2],'Xtick',([x1:0.25:x2]),'XTickLabel',{' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';},'FontSize',fsize,'box','on');
text (0.01*(x2-x1)+x2, altcb, 'Te(K)','FontSize',fsize);
pos1=get(h1,'position');
set(h1,'position',pos1)
hp=colorbar;
pos=get(hp,'position');
pos(4) = 0.19;
set(hp,'position',pos);
pos1(3)=0.80;
set(gca,'position',pos1);

h1=subplot('position',[0.08,0.285,0.87,0.225]);	% 3/4
hold on;
%contourf(ht_timewd,ht_tie,ti_tie,11);
contourf(htmid_time,ht_mid,ti_tie,tnlev);
[min1,max1] = caxis;
min1,max1
caxis([200.,1850.]);
plot(xtot,z300,':k');
xlabel(' ','FontSize',fsize);
set (gca,'xlim',[x1,x2],'ylim',[alt1,alt2],'Xtick',([x1:0.25:x2]),'XTickLabel',{' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';},'FontSize',fsize,'box','on');
text(0.01*(x2-x1)+x2, altcb, 'Ti (K)','FontSize',fsize);
pos1=get(h1,'position');
set(h1,'position',pos1)
hp=colorbar;
pos=get(hp,'position');
pos(4) = 0.19;
set(hp,'position',pos);
pos1(3)=0.80;
set(gca,'position',pos1);

h1=subplot('position',[0.08,0.06,0.87,0.225]);	% 4/4
hold on;
%contourf(timewd,htav_tie,tn_tie,13);
%contourf(htmid_time,ht_mid,tn_tie,tnlev);
 contourf(htmid_time,ht_mid,wi_tie,wilev);
[min1,max1] = caxis;
min1,max1
%caxis([200.,1850.]);
caxis([-100.,400.]);
plot(xtot,z300,':k');
xlabel ('Daynumber','FontSize',fsize);
set (gca,'xlim',[x1,x2],'ylim',[alt1,alt2],'Xtick',([x1:0.25:x2]),'FontSize',fsize,'box','on');
%text(0.01*(x2-x1)+x2, altcb, 'Tn (K)','FontSize',fsize);
text(0.01*(x2-x1)+x2, altcb, 'Wi m/s','FontSize',fsize);
hold off;
pos1=get(h1,'position');
set(h1,'position',pos1)
hp=colorbar;
pos=get(hp,'position');
pos(4) = 0.19;
set(hp,'position',pos);
pos1(3)=0.80;
set(gca,'position',pos1);

hold off;

set(gcf,'PaperPosition',[0.5, 0.5, 7.5, 10.5]);
set(gcf,'PaperPositionmode', 'auto');
saveas(gcf,plottie,'epsc')


%% Plot mad, residuals, and median over npmad pts (Te/Ti/Tn, Ne300/NmF2/TEC?, hmF2)
ymadx = 300.;
yr1 = -199.;
yr2 = 199.;


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

subplot('position',[0.08,0.720,0.87,0.240])	% 1/3
hold on;
title(daytit,'FontSize',fsize);
plot (timemad,tnmad_wd,'Color','k','Linestyle','--');
plot (timemad,timad_wd,'Color','c','Linestyle','--');
plot (timemad,temad_wd,'Color','r','Linestyle','--');
 plot (timemadti_obs,timad_obs,'Color','c');
 plot (timemadte_obs,temad_obs,'Color','r');
ylabel('median abs deviation (m/s)', 'FontSize',fsize);
xlabel(' ','FontSize',fsize);
set (gca,'xlim',[x1,x2],'ylim',[0.,ymadx],'XTickLabel',{' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';},'FontSize',fsize,'box','on');

subplot('position',[0.08,0.480,0.87,0.240])	% 2/3
hold on;
plot (timemad,tnresid_wd,'Color','k','Linestyle','--');
plot (timemad,tiresid_wd,'Color','c','Linestyle','--');
plot (timemad,teresid_wd,'Color','r','Linestyle','--');
 plot (timemadti_obs,tiresid_obs,'Color','c','LineStyle','-');
 plot (timemadte_obs,teresid_obs,'Color','r','LineStyle','-');
plot(xtot,zero,'--k');
ylabel('Residual (m/s)', 'FontSize',fsize);
xlabel(' ','FontSize',fsize);
set (gca,'xlim',[x1,x2],'ylim',[yr1,yr2],'XTickLabel',{' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';},'FontSize',fsize,'box','on');

subplot('position',[0.08,0.050,0.87,0.430])	% 3/3
hold on;
text (x1+0.80*(x2-x1),y95,ht300str,'Color','k','FontSize',f9);
text (x1+0.80*(x2-x1),y90,'20min 4h med Obs','Color','k','FontSize',f9);
text (x1+0.01*(x2-x1),y95,'Tn','Color','k','FontSize',f9);
text (x1+0.09*(x2-x1),y95,'Ti','Color','c','FontSize',f9);
text (x1+0.19*(x2-x1),y95,'Te','Color','r','FontSize',f9);
 madobs = sprintf('%4.2f',mad_tiobs);
 plot (timemadti_obs,timed_obs,'Color','c','LineStyle','-');
 text (x1+0.10*(x2-x1),y90,madobs,'Color','c','FontSize',f9);
 madobs = sprintf('%4.2f',mad_teobs);
 plot (timemadte_obs,temed_obs,'Color','r','LineStyle','-');
 text (x1+0.19*(x2-x1),y90,madobs,'Color','r','FontSize',f9);
 madwd = sprintf('%4.2f',mad_tn_wd);
 plot (timemad,tnmed_wd,'Color','k','LineStyle','--');
text (x1+0.01*(x2-x1),y85,madwd,'Color','k','FontSize',f9);
 madwd = sprintf('%4.2f',mad_ti_wd);
 plot (timemad,timed_wd,'Color','c','LineStyle','--');
text (x1+0.10*(x2-x1),y85,madwd,'Color','c','FontSize',f9);
 madwd = sprintf('%4.2f',mad_te_wd);
 plot (timemad,temed_wd,'Color','r','LineStyle','--');
text (x1+0.19*(x2-x1),y85,madwd,'Color','r','FontSize',f9);
text (x1+0.80*(x2-x1),y85,'--av TIEGCM ','Color','k','FontSize',f9);
plot(xtot,zero,'--k');
ytit = sprintf ('%d hour median T (K)',numhr);
ylabel (ytit,'FontSize',fsize);
xlabel('Daynumber','FontSize',fsize);
set (gca,'xlim',[x1,x2],'ylim',[y1,y2],'FontSize',fsize,'box','on');
hold off;

set(gcf,'PaperPosition',[0.5, 0.5, 7.5, 10.5]);
set(gcf,'PaperPositionmode', 'auto');
saveas(gcf,plotmadt,'epsc')


% Plot mad of NmF2, Ne300km, (TEC?,hmax?)
figure;
% Set size to be 8.5x11 inches
set(gcf,'position',[142    99   647   822])
H = gcf;

ymadx = 0.2;
subplot('position',[0.08,0.720,0.87,0.240])	% 1/3
hold on;
title(daytit,'FontSize',fsize);
plot (timemad,nemad_wd,'Color','k','Linestyle','--');
plot (timemad,tecmad_wd,'Color','c','Linestyle','--');
plot (timemad,nmaxmad_wd,'Color','r','Linestyle','--');
 plot (timemadne_obs,nemad_obs,'Color','k');
 plot (timemadn_obs,nmaxmad_obs,'Color','r');
 plot (timemadt_obs,tecmad_obs,'Color','c');
ylabel('median abs deviation (m/s)', 'FontSize',fsize);
xlabel(' ','FontSize',fsize);
set (gca,'xlim',[x1,x2],'ylim',[0.,ymadx],'XTickLabel',{' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';},'FontSize',fsize,'box','on');

yr1 = -0.2;
yr2 = 0.2;
subplot('position',[0.08,0.480,0.87,0.240])	% 2/3
hold on;
plot (timemad,neresid_wd,'Color','k','Linestyle','--');
plot (timemad,tecresid_wd,'Color','c','Linestyle','--');
plot (timemad,nmaxresid_wd,'Color','r','Linestyle','--');
 plot (timemadne_obs,neresid_obs,'Color','k','LineStyle','-');
 plot (timemadt_obs,tecresid_obs,'Color','c','LineStyle','-');
 plot (timemadn_obs,nmaxresid_obs,'Color','r','LineStyle','-');
plot(xtot,zero,'--k');
ylabel('Residual (m/s)', 'FontSize',fsize);
xlabel(' ','FontSize',fsize);
set (gca,'xlim',[x1,x2],'ylim',[yr1,yr2],'XTickLabel',{' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';},'FontSize',fsize,'box','on');

y1 = 9.01;
y2 = 17.99;
y95 = y1 + 0.62*(y2-y1);
y90 = y1 + 0.55*(y2-y1);
y85 = y1 + 0.48*(y2-y1);

subplot('position',[0.08,0.050,0.87,0.430])	% 3/3
hold on;
text (x1+0.80*(x2-x1),y95,ht300str,'Color','k','FontSize',f9);
text (x1+0.80*(x2-x1),y90,'20min 4h med Obs','Color','k','FontSize',f9);
text (x1+0.01*(x2-x1),y95,'Ne','Color','k','FontSize',f9);
text (x1+0.09*(x2-x1),y95,'TEC','Color','c','FontSize',f9);
text (x1+0.19*(x2-x1),y95,'Nmax','Color','r','FontSize',f9);
 madobs = sprintf('%4.2f',mad_neobs);
plot (timemadne_obs,nemed_obs,'Color','k','LineStyle','-');
%plot (timemadne_obs,10.^nemed_obs,'Color','k','LineStyle','-');
 text (x1+0.01*(x2-x1),y90,madobs,'Color','k','FontSize',f9);
 madobs = sprintf('%4.2f',mad_tecobs);
 plot (timemadt_obs,tecmed_obs,'Color','c','LineStyle','-');
% plot (timemadt_obs,10.^tecmed_obs,'Color','c','LineStyle','-');
 text (x1+0.10*(x2-x1),y90,madobs,'Color','c','FontSize',f9);
 madobs = sprintf('%4.2f',mad_nmaxobs);
 plot (timemadn_obs,nmaxmed_obs,'Color','r','LineStyle','-');
% plot (timemadn_obs,10.^nmaxmed_obs,'Color','r','LineStyle','-');
 text (x1+0.19*(x2-x1),y90,madobs,'Color','r','FontSize',f9);
 madwd = sprintf('%4.2f',mad_ne_wd);
 plot (timemad,nemed_wd,'Color','k','LineStyle','--');
% plot (timemad,10.^nemed_wd,'Color','k','LineStyle','--');
text (x1+0.01*(x2-x1),y85,madwd,'Color','k','FontSize',f9);
 madwd = sprintf('%4.2f',mad_tec_wd);
 plot (timemad,tecmed_wd,'Color','c','LineStyle','--');
% plot (timemad,10.^tecmed_wd,'Color','c','LineStyle','--');
text (x1+0.10*(x2-x1),y85,madwd,'Color','c','FontSize',f9);
 madwd = sprintf('%4.2f',mad_nmax_wd);
 plot (timemad,nmaxmed_wd,'Color','r','LineStyle','--');
% plot (timemad,10.^nmaxmed_wd,'Color','r','LineStyle','--');
text (x1+0.19*(x2-x1),y85,madwd,'Color','r','FontSize',f9);
text (x1+0.80*(x2-x1),y85,'-- TIEGCM ','Color','k','FontSize',f9);
plot(xtot,zero,'--k');
ytit = sprintf ('%d hour median log10Ne,Nmax(m-3) TEC(m-2)',numhr);
ylabel (ytit,'FontSize',fsize);
xlabel('Daynumber','FontSize',fsize);
set (gca,'xlim',[x1,x2],'ylim',[y1,y2],'FontSize',fsize,'box','on');
hold off;


set(gcf,'PaperPosition',[0.5, 0.5, 7.5, 10.5]);
set(gcf,'PaperPositionmode', 'auto');
saveas(gcf,plotmadn,'epsc')
end


