%%
%% This program is written for comparing the TIEGCM and MLH 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

 year1 = 2005;
 year1 = 2006;
if (year1 < 2005.5)
 day1 = 243;
 da1 = 243 + 10./24.;
 day2 = 245;
 da2 = 244 + 12./24.;
 n15m = 67;
else
 day1 = 348;
 da1 = 348 + 12./24.;
 day2 = 350;
 da2 = 350 + 0./24.;
 n15m = 192;
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;


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

% Look at altitude 300 km or range from 280-320 km 
ht1 = 280.e+5;
ht2 = 320.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


% Millstone Hill coordinates (42.62N, 288.51E, j,i=26,23)  36 lats from -87.5 to +87.5, 72 lons from -180 to +175(?)
jj = 26;
ii = 23;


% Have MLH obs start at 0020 UT for the 20 min values to be the same as the TIEGCM
nalt = 34;
n300 = 22;
for n=1:nobs
 n300_tn_20min(n) = 0;
 n300_ti_20min(n) = 0;
 n300_te_20min(n) = 0;
 n300_ne_20min(n) = 0;
 tn_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;
% 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)
 obs1fil = ['/hao/aim2/emery/cetic11/mlh050901_nmax'];
 obs2fil = ['/hao/aim2/emery/cetic11/mlh050901_profiles'];
else
 obs1fil = ['/hao/aim2/emery/cetic11/mlh061214_15_nmax'];
 obs2fil = ['/hao/aim2/emery/cetic11/mlh061214_15_profiles'];
end

%%  Read MLH 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);
% Skip over data description in mlh061214_15_nmax
if (year1 > 2005.5)
 for n=1:60
  linobs1 = fgetl(OBS1);
 end
end

for n=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) 
 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


% Find nearest 20 min time 
 ntfrac = fix( ((timeobs(n)-time_obs_20min(1))*1440. + 9.99)/20. ) + 1;
if (ntfrac >= 1)
 hmax_obs_nan_20min(ntfrac) = hmaxobs(n);
 nmax_obs_nan_20min(ntfrac) = nmaxobs(n);
 tec_obs_nan_20min(ntfrac) = tecobs(n);
end

% Read 2D profiles
 tnobs300(n) = NaN;
 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 %f');
   if (n == 1)
     altobs(n9) = obs2dat(4);
   end
   tnobs(n9,n) = NaN;
   tiobs(n9,n) = NaN;
   teobs(n9,n) = NaN;
   neobs(n9,n) = NaN;
 if (obs2dat(5) > 0.)
  tnobs(n9,n) = obs2dat(5);
 end
 if (obs2dat(6) > 0.)
  tiobs(n9,n) = obs2dat(6);
 end
 if (obs2dat(8) > 0.)
  teobs(n9,n) = obs2dat(8);
 end
 if (obs2dat(10) > 0.)
  neobs(n9,n) = log10(obs2dat(10));
 end
if (n9 == n300)
   if (obs2dat(5) > 0.)
  tnobs300(n) = obs2dat(5);
 end
 if (obs2dat(6) > 0.)
  tiobs300(n) = obs2dat(6);
 end
 if (obs2dat(8) > 0.)
  teobs300(n) = obs2dat(8);
 end
 if (obs2dat(10) > 0.)
  neobs300(n) = log10(obs2dat(10));
 end

% Find nearest 20 min time every 15 min
 ntfrac = fix( ((timeobs(n)-time_obs_20min(1))*1440. + 9.99)/20. ) + 1;
if (ntfrac >= 1)
 tn_obs_nan_20min(ntfrac) = tnobs300(n);
 ti_obs_nan_20min(ntfrac) = tiobs300(n);
 te_obs_nan_20min(ntfrac) = teobs300(n);
 ne_obs_nan_20min(ntfrac) = neobs300(n);
end
end	% end of n9=n300
end	% end of if daynumobs between day1 and day2-1
end	% end of for n=1:n15m

end	% end n9=1:nalt

% Calculate mad from tn_obs_nan_20min
nmad1 = fix((timeobs(1)-time_obs_20min(1))*1440./20.) + 1;
nmad2 = min(nobs-npmad+1,fix((timeobs(n)-time_obs_20min(1))*1440./20.) + 1);
minmad1 = max(1,nmad1);
maxmad2 = min(nobs-midm1-1,nmad2);
iob = 0;
nmad = 0;
kob = 0;
nob = 0;
job = 0;
jmad = 0;
for i=minmad1:maxmad2
 icen = i+midm1;
 [obsmed,obsresid,obsmad] = get_mad(tn_obs_nan_20min,npmad,i,midm1);
 if ~isnan(obsmed)
   iob = iob + 1;
   nmad = nmad + 1;
   tntimemad_obs(iob) = time_obs_20min(icen);
   tnmed_obs(iob) = obsmed;
   tnresid_obs(iob) = obsresid;
   tnmad_obs(iob) = obsmad;
 end
 [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(obstimed)
   job = job + 1;
   jmad = jmad + 1;
   timemad_obs(job) = time_obs_20min(icen);
   timed_obs(job) = obstimed;
   tiresid_obs(job) = obstiresid;
   timad_obs(job) = obstimad;
   temed_obs(job) = obstemed;
   teresid_obs(job) = obsteresid;
   temad_obs(job) = obstemad;
   nemed_obs(job) = obsnemed;
   neresid_obs(job) = obsneresid;
   nemad_obs(job) = 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
iob
% Calculate total mads:
mad_tnobs = 0.;
if (iob > 0)
 for i=1:iob
  if ~isnan(tnmad_obs(i))
   mad_tnobs = mad_tnobs + tnmad_obs(i);
  end
 end
 mad_tnobs = mad_tnobs/nmad;
end
nmad_tnobs = nmad;
mad_tiobs = 0.;
mad_teobs = 0.;
mad_neobs = 0.;
if (job > 0)
 for i=1:job
  if ~isnan(timad_obs(i))
   mad_tiobs = mad_tiobs + timad_obs(i);
   mad_teobs = mad_teobs + temad_obs(i);
   mad_neobs = mad_neobs + nemad_obs(i);
  end
 end
 mad_tiobs = mad_tiobs/jmad;
 mad_teobs = mad_teobs/jmad;
 mad_neobs = mad_neobs/jmad;
end
nmad_tiobs = 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.s05243_244_001.nc'];
else
%    Fname=[fpath '/wd194/EMERY.tiegcm_dyncrit.s06348_350_001.nc'];
    Fname=[fpath '/wd194/wd.s06348_349_001.nc'];
end
    [nt,dayfrac,tnav,tiav,teav,neav,nmax,hmax,tec,nz,ht,htne,htte,htti,htmid,httn] = get_netein_nc(Fname,ht1,ht2,ii,jj);
 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);
 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)
   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.s05243_244_002.nc'];
else
%    Fname=[fpath '/wd194/EMERY.tiegcm_dyncrit.s06348_350_002.nc'];
    Fname=[fpath '/wd194/wd.s06348_349_002.nc'];
end
    [nt,dayfrac,tnav,tiav,teav,neav,nmax,hmax,tec,nz,ht,htne,htte,htti,htmid,httn] = get_netein_nc(Fname,ht1,ht2,ii,jj);
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);
 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)
   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);


%% calculate rms and print from obs every 20 min 

 filerms = sprintf('rms_dec06_mlh_tie');
RMS = fopen(filerms,'w');
fprintf(RMS, 'RMS of 20-min TIEGCM values of hmF2, h=NmF2, TEC from avs every 20 min of Millstone Hill ISR obs from 280-320 km \n');
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);
 rmstn300wd = 0.;
 [rmstn300wd,ntn300rms] = get_rms(tn_obs_nan_20min,tn_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,rmstn300wd,ntn300rms);

% 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_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);
 [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_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 Millston Hill ISR obs and TIEGCM values \n');
fprintf(RMS, '%     obs   wd1.94 hmax nmax tec  300km: ne te ti tn \n',year1);
fprintf(RMS, '%d-%d %4.0d %13.4e %13.4e %13.4e %13.4e %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);

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 = -20;
yi2 = 20;
yi95 = yi1 + 0.92*(yi2-yi1);
yitot = [yi1,yi2];

 y1=500;
 y2 = 2899;
y95 = y1 + 0.92*(y2-y1);
y90 = y1 + 0.85*(y2-y1);
y85 = y1 + 0.78*(y2-y1);

yh1 = 150;
yh2 = 499;

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

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

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


 x1 = day1;
 x2 = day2;
if (year1 < 2005.5)
 daytit = sprintf('Millstone Hill (42.62N, 288.51E) 01 Sep 2005');
 plottit = sprintf('mlh_05244');
 plotht = sprintf('mlhht_05244');
 plottie = sprintf('mlhtie_05244');
 plotmadt = sprintf('mlhmadt_05244');
 plotmadn = sprintf('mlhmadn_05244');
 if (notie==1)
  plottit = sprintf('mlh_data_05244');
 end
else
 daytit = sprintf('Millstone Hill (42.62N, 288.51E) 14-15 Dec 2006');
 plottit = sprintf('mlh_06348_9');
 plotht = sprintf('mlhht_06348_9');
 plottie = sprintf('mlhtie_06348_9');
 plotmadt = sprintf('mlhmadt_06348_9');
 plotmadn = sprintf('mlhmadn_06348_9');
 if (notie==1)
  plottit = sprintf('mlh_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;
%plot (time_obs_20min,tn_obs_nan_20min,'k','LineWidth',2.0);
%text (x1+0.82*(x2-x1),y95,'All obs (err<100K)','Color','k','FontSize',f9,'LineWidth',2.0);
 plot (timeobs,tnobs300,'Color','k');
 plot (timeobs,tiobs300,'Color','c');
 plot (timeobs,teobs300,'Color','r');
%plot (time_obs_20min,tn_obs_nan_20min,'Color','k','LineStyle',':');
%plot (time_obs_20min,ti_obs_nan_20min,'Color','c','LineStyle',':');
%plot (time_obs_20min,te_obs_nan_20min,'Color','r','LineStyle',':');
 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);
rms = sprintf('%4.2f',rmstn300wd);
text (x1+0.01*(x2-x1),y90,rms,'Color','k','FontSize',f9);
text (x1+0.28*(x2-x1),y95,'--av TIEGCM ','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('300 km 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/4-5
hold on;
 plot (timeobs,hmaxobs,'Color','r');
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.240,0.87,0.190])	% 4/4-5
hold on;
plot (timeobs,nmaxobs,'Color','r');
plot (timeobs,neobs300,'Color','k');
text (x1+0.01*(x2-x1),ye95,'Nmax','Color','r','FontSize',f9);
text (x1+0.10*(x2-x1),ye95,'300 km','Color','k','FontSize',f9);
text (x1+0.55*(x2-x1),ye95,'Contact John Holt for ISR data usage','Color','k','FontSize',f9);
if (notie ~= 1)
 plot (timewd,ne_wd,'Color','k','LineStyle','--');
 plot (timewd,nmax_wd,'r','LineStyle','--');
 text (x1+0.19*(x2-x1),ye95,'-- TIEGCM','Color','k','FontSize',f9);
end
ylabel('log10(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.190])	% 5/4-5
hold on;
plot (timeobs,tecobs,'Color','r');
if (notie ~= 1)
text (x1+0.19*(x2-x1),yt95,'-- TIEGCM','Color','k','FontSize',f9);
 plot (timewd,tec_wd,'r','LineStyle','--');
end
ylabel('Log10(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 TIEGCM
% Define contour levels for Ne, Te, Ti-Tn
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
contourf(timeobs,altobs,neobs,nelev);
[min1,max1] = caxis;
min1,max1
caxis([9.,12.]);
 plot (timeobs,hmaxobs,'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);
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(timeobs,altobs,teobs,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(timeobs,altobs,tiobs,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(timeobs,altobs,tnobs,tnlev);
[min1,max1] = caxis;
min1,max1
caxis([200.,1850.]);
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);
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,plotht,'epsc')

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);
plot(xtot,z300,':k');
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);
[min1,max1] = caxis;
min1,max1
caxis([200.,1850.]);
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);
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 = 200.;
yr1 = -129.;
yr2 = 129.;


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 (tntimemad_obs,tnmad_obs,'Color','k');
 plot (timemad_obs,timad_obs,'Color','c');
 plot (timemad_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 (tntimemad_obs,tnresid_obs,'Color','k','LineStyle','-');
 plot (timemad_obs,tiresid_obs,'Color','c','LineStyle','-');
 plot (timemad_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),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_tnobs);
 plot (tntimemad_obs,tnmed_obs,'Color','k','LineStyle','-');
 text (x1+0.01*(x2-x1),y90,madobs,'Color','k','FontSize',f9);
 madobs = sprintf('%4.2f',mad_tiobs);
 plot (timemad_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 (timemad_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 Tn (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 = 1.;
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 (timemad_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 (timemad_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 = ye1;
y2 = yt2;
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),y90,'20min 4h med Obs','Color','k','FontSize',f9);
text (x1+0.01*(x2-x1),y95,'Ne300','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 (timemad_obs,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','-');
 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','-');
 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','--');
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','--');
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','--');
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 log10 Ne300,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


