%%
%% This program is written for plotting the TIEGCM dres TEC data at 8 lonsra
%% BA Emery, October 2011 from /aim3/emery/tiegcm_dres/plots/tec8lon.m 
%% Revised for climatology output November 2011
%% March 2012 revised for IGS 'global' TEC for 8 glons
clear all; close all

dark_green = [0 0.5 0];
green_blue = [0 0.75 0.50];
purple = [0.48 0.06 0.89];
light_orange = [1 0.69 0.39];
yellow_green = [0.75 0.75 0];
light_purple = [0.7 0.78 1];
aqua = [0.04 0.52 0.78];

% Data
 loncutsE=[25,90,140,175,200,250,285,345]; % specify climatology longitudes to plot

% Model

 year1 = 2007;
 nday1 = 325;
% nday2 = 366+26;
 nday2 = 365+20+1;

% Open Viz averages from Jicamarca magnetometer
filviz = sprintf('/hao/aim3/emery/cetic_climatology/jicmag_viz1121')
datviz = dlmread(filviz,'',0,0);
% Open Viz averages from JULIA
filjul = sprintf('/hao/aim3/emery/cetic_climatology/julia_viz1121')
datjul = dlmread(filjul,'',0,0);
% Calculate averages
julam = 0.;
julpm = 0.;
magam = 0.;
magpm = 0.;
njula = 0;
njulp = 0;
nmag = 0;
for nd=1:61
if (datjul(nd,3)>-99.)
 njula = njula + 1;
 julam = julam + datjul(nd,3);
end
if (datjul(nd,5)>-99.)
 njulp = njulp + 1;
 julpm = julpm + datjul(nd,5);
end
if (datviz(nd,3)>-99.)
 nmag = nmag + 1;
 magam = magam + datviz(nd,3);
 magpm = magpm + datviz(nd,5);
end
end
julam = julam/njula
julpm = julpm/njulp
magam = magam/nmag
magpm = magpm/nmag



%  Open 8 glon 'global' TEC daily values from IGS TEC
for j=1:8
filetec = sprintf('/hao/aim3/emery/cetic_climatology/igs_tec/global_tec_igs_07325_08020_lon%3.3d-%3.3d.dat',loncutsE(j),loncutsE(j)+5);
datnam = dlmread(filetec,'',0,0);
for nd=1:61
 datigs(nd,j) = 0.5*(datnam(nd,6)+datnam(nd,8));
end
end		% j=1:8
nd = 20
datigs(nd,1),datigs(nd,4)

% Open 8 glon 'global' TEC daily values from MIT GPS TEC
for j=1:8
filemit = sprintf('/hao/aim3/emery/cetic_climatology/mit_gps/global_tec_mit_07325_08020_lon%3.3d.dat',loncutsE(j));
datgps = dlmread(filemit,'',0,0);
for nd=1:61
 datmit(nd,j) = datgps(nd,3);
 tecmit(nd) = datgps(nd,3);
 areapc(nd) = datgps(nd,4);
end
 areamit(j) = median(areapc);
 mittecmed(j) = median(tecmit);
end		% j=1:8
areamit,mittecmed

% Open 'global' CHAMP density files from ascent/descent LTs (even wts) (08006-9 set to missing)
ifsat = 0;
if (ifsat==0)
 filcad = sprintf('/hao/aim3/emery/cetic_climatology/champ_density/global_400km_LTden_07325_08020c.dat')
else
% filcad = sprintf('/hao/aim3/emery/cetic_climatology/champ_density/global_sat_LTden_07300_08045_a.dat')
 filcad = sprintf('/hao/aim3/emery/cetic_climatology/champ_density/global_sat_LTden_07325_08020_b.dat')
end
datcad = dlmread(filcad,'',0,0);
nl=length(datcad)
for nd=1:nl
% CHAMP 'global' density from LT paths in ascent/descent
 champday(nd) = datcad(nd,2)+0.5;
 if (datcad(nd,1)==2008)
  champday(nd) = datcad(nd,2)+365.5;
 end
 if (ifsat==1)
 dencad(nd) = datcad(nd,13);
 denasc(nd) = datcad(nd,5);
 dendesc(nd) = datcad(nd,10);
 dencad2(nd) = datcad(nd,14);
 denasc2(nd) = datcad(nd,6);
 dendesc2(nd) = datcad(nd,11);
 dencadmod(nd) = datcad(nd,15);
 denascmod(nd) = datcad(nd,7);
 dendescmod(nd) = datcad(nd,12);
 rmdasc(nd) = datcad(nd,6)/datcad(nd,5);
 rmddesc(nd) = datcad(nd,11)/datcad(nd,10);
 rmoddasc(nd) = datcad(nd,7)/datcad(nd,5);
 rmodddesc(nd) = datcad(nd,12)/datcad(nd,10);
 else
% 400 km data from Sutton (no MSIS00 results at 400 km)
 dencad(nd) = datcad(nd,11);
 denasc(nd) = datcad(nd,5);
 dendesc(nd) = datcad(nd,9);
 end
end
% Read Luehr CHAMP data 61-4=57 days (missing 07356, 08006-8)
% filuhr = sprintf('/hao/aim3/emery/cetic_climatology/champ_density/global_sat_400km_LTden_luehr_07325_08020_wo4.dat')
%datluhr = dlmread(filuhr,'',1,0);
% Read Eelco Doornbos data as of 28 Sep 2012
 filuhr = sprintf('/hao/aim3/emery/cetic_climatology/champ_density/global_LTden_doornbos_07325_08020.dat2')
datluhr = dlmread(filuhr,'',0,0);
nll = length(datluhr)
for nd=1:nll
% CHAMP 'global' density from LT paths in ascent/descent
 luhrday(nd) = datluhr(nd,2)+0.5;
 if (datluhr(nd,1)==2008)
  luhrday(nd) = datluhr(nd,2)+365.5;
 end
 if ifsat==1
% 332-354 km data from Luehr (lower in value and in ht from Sutton at 337-368!)
%  because Katrin Haeusler used a code to convert sat xyz with spherical earth
  ldenasc(nd) = datluhr(nd,10);
  ldendesc(nd) = datluhr(nd,19);
  lmsisasc(nd) = datluhr(nd,11);
  lmsisdesc(nd) = datluhr(nd,20);
  ljb08asc(nd) = datluhr(nd,12);
  ljb08desc(nd) = datluhr(nd,21);
  lhasdmasc(nd) = datluhr(nd,13);
  lhasdmdesc(nd) = datluhr(nd,22);
 else
% 400 km data from Luehr
  ldenasc(nd) = datluhr(nd,14);
  ldendesc(nd) = datluhr(nd,23);
  lmsisasc(nd) = datluhr(nd,15);
  lmsisdesc(nd) = datluhr(nd,24);
  ljb08asc(nd) = datluhr(nd,16);
  ljb08desc(nd) = datluhr(nd,25);
  lhasdmasc(nd) = datluhr(nd,17);
  lhasdmdesc(nd) = datluhr(nd,26);
 end
  ldencad(nd) = 0.5*(ldenasc(nd)+ldendesc(nd));
  lmsiscad(nd) = 0.5*(lmsisasc(nd)+lmsisdesc(nd));
  ljb08cad(nd) = 0.5*(ljb08asc(nd)+ljb08desc(nd));
  lhasdmcad(nd) = 0.5*(lhasdmasc(nd)+lhasdmdesc(nd));
end

% Open Emmert's 400km file
fil400e = sprintf('/home/emery/cetic/tledens_2007_2010_400km.txt')
dat400e = dlmread(fil400e,'',300,0);
for nd=1:300
  day400(nd) = dat400e(nd,2)+0.5;
  if (dat400e(nd,1)==2008)
   day400(nd) = dat400e(nd,2)+365.5;
  end
  drag400(nd) = 10.^dat400e(nd+1,7);
  msis400(nd) = 10.^dat400e(nd,6);
end
 


% Open global density file (last 2 lines are median and mean)
% daynum(2007) 250km_drag,msis,model 400km_drag,msis,model 550km_drag,msis,model 
%filden = sprintf('/hao/aim3/emery/cetic_climatology/daily_denc')
% daynum(2007) 250km_drag,msis,model 400km_drag,msis,model 550km_drag,msis,model globleTEC
filden = sprintf('/hao/aim3/emery/cetic_climatology/daily_den_tec')
datden = dlmread(filden,'',1,0);

% 11/12 Open CCMC TIE 1.93 Heelis Kp
filden3 =sprintf('/hao/aim3/emery/cetic_climatology/daily_den_tec_tie193')
datden3 = dlmread(filden3,'',1,0);

% 11/12 Open global UAM files and TIME-AMIE
% Skip UAM5 since results are poor for smin
%filden5 = sprintf('/hao/aim3/emery/cetic_climatology/daily_den_tec_UAM5')
%datden5 = dlmread(filden5,'',2,0);
filtec4 = sprintf('/hao/aim3/emery/cetic_climatology/daily_tec500_tecgps_UAM4')
dattec4 = dlmread(filtec4,'',2,0);
filden4 = sprintf('/hao/aim3/emery/cetic_climatology/daily_den_tec_UAM4')
datden4 = dlmread(filden4,'',2,0);
fildenta = sprintf('/hao/aim3/emery/cetic_climatology/daily_den_tec_TIMA')
datdenta = dlmread(fildenta,'',2,0);

% Open daily IMF file
filimf = sprintf('/hao/aim3/emery/cetic_climatology/daily_imf_07325_08026')
datimf = dlmread(filimf,'',0,0);

% 25 daily Kp>=2, Vsw>=500: 325-330(6), 345-346(2), 351-357(7), 05-08(4), 14-19(6); 25(+1)
kp2x = 2.0;
vsw2x = 500.;
% 25 daily Kp<=1, vsw<=450: 332-343(12), 349-350(2), 358-360(3), 362-364(3), 01-04(4), 11(1); 22-23(+2)
kp1n = 1.0;
vsw1n = 450.;
nd = 0;
nd1 = 0;
nd2 = 0;
for ida=nday1:nday2-1
 nd = nd + 1;
 midday(nd) = ida+0.5;
 kp(nd) = datimf(nd,4); 
 vsw(nd) = datimf(nd,36);
 dsw(nd) = datimf(nd,32);
 bt(nd) = datimf(nd,15);
 bz(nd) = datimf(nd,24);
 by(nd) = datimf(nd,21);
 f107(nd) = datimf(nd,10);
% Set daynumber for 3-day satellite drag 1 day earlier (den from nd+1)
 lg10drag(nd) = datden(nd+1,5);
 lg10msis(nd) = datden(nd,6);
 lg10denmod(nd) = datden(nd,7);
 drag(nd) = 10^datden(nd+1,5);
 dragviz(nd) = drag(nd)*1.e+13;
 msis(nd) = 10^datden(nd,6);
 denmod(nd) = 10^datden(nd,7);
% For files with global TEC from model
 gtec(nd) = datden(nd,11);
% For UAM5 and TIME-AMIE and CCMC TIE1.93 Heelis Kp
 gtec3(nd) = datden3(nd,11);
 gden3(nd) = 10^datden3(nd,7)
 gtecta(nd) = datdenta(nd,8);
 gdenta(nd) = 10^datdenta(nd,2);
%gtec5(nd) = datden5(nd,8);
%gden5(nd) = 10^datden5(nd,2);
% Use gps ht for TEC UAM4
 gtec4(nd) = dattec4(nd,8);
 gtec4shp(nd) = dattec4(nd,9);
 gtec4shm(nd) = dattec4(nd,10);
 gtec4eqb(nd) = dattec4(nd,11);
 gtec4nhm(nd) = dattec4(nd,12);
 gtec4nhp(nd) = dattec4(nd,13);
% cf w 500km for TECl for low
 gtec4l(nd) = dattec4(nd,2);
 gtec4shpl(nd) = dattec4(nd,3);
 gtec4shml(nd) = dattec4(nd,4);
 gtec4eqbl(nd) = dattec4(nd,5);
 gtec4nhml(nd) = dattec4(nd,6);
 gtec4nhpl(nd) = dattec4(nd,7);
% cf w 200km for TECm for medium
 gtec4m(nd) = datden4(nd,8);
 gtec4shpm(nd) = datden4(nd,9);
 gtec4shmm(nd) = datden4(nd,10);
 gtec4eqbm(nd) = datden4(nd,11);
 gtec4nhmm(nd) = datden4(nd,12);
 gtec4nhpm(nd) = datden4(nd,13);
 gden4(nd) = 10^datden4(nd,2);

% Divide into days Kp>=kp2(26) and Kp<=kp1(28)
%  Add vsw2(25) and vsw1(25) conditions
if (kp(nd)>=kp2x & vsw(nd)>=vsw2x)
  nd2 = nd2 + 1;
  midday2(nd2) = midday(nd);
  f1072(nd2) = f107(nd);
  vsw2(nd2) = vsw(nd);
  dsw2(nd2) = dsw(nd);
  bt2(nd2) = bt(nd);
  bz2(nd2) = bz(nd);
  by2(nd2) = by(nd);
  kp2(nd2) = kp(nd);
  drag2(nd2) = drag(nd);
  msis2(nd2) = msis(nd);
  denmod2(nd2) = denmod(nd);
  gtec2(nd2) = gtec(nd);
end
if (kp(nd)<=kp1n & vsw(nd)<=vsw1n)
  nd1 = nd1 + 1;
  midday1(nd1) = midday(nd);
  f1071(nd1) = f107(nd);
  vsw1(nd1) = vsw(nd);
  dsw1(nd1) = dsw(nd);
  bt1(nd1) = bt(nd);
  bz1(nd1) = bz(nd);
  by1(nd1) = by(nd);
  kp1(nd1) = kp(nd);
  drag1(nd1) = drag(nd);
  msis1(nd1) = msis(nd);
  denmod1(nd1) = denmod(nd);
  gtec1(nd1) = gtec(nd);
end

end	% for ida=nday1:nday2
mean(gden4),mean(gtec4)
nd2,mean(vsw2),mean(dsw2),mean(bt2),mean(bz2),mean(by2),mean(kp2),mean(drag2),mean(msis2),mean(denmod2),mean(f1072),mean(gtec2)
midday2
nd1,mean(vsw1),mean(dsw1),mean(bt1),mean(bz1),mean(by1),mean(kp1),mean(drag1),mean(msis1),mean(denmod1),mean(f1071),mean(gtec1)
midday1

% Plot global densities
% 09027 msis=[-10.49505 -12.15793 -13.44664] drag=[-10.62822 -12.38078 -13.71003]
% tiesaber g/cm3 = [3.88e-14 7.85e-16 NA], lg10(kg/m3) = [-10.4106 -12.1053 NA]  (like MSIS)
dark_green = [0 0.5 0];

figure
zero = [0.,0.];
xtot = [nday1,nday2];

fsize = 12;
% 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.54,0.87,0.40])	% 1/2
subplot('position',[0.10,0.82,0.87,0.12])	% 1/4 of 1/2
hold on;
%title('Median Daily Solar Wind Velocity','FontSize',15,'FontWeight','bold');
title('Median Daily Geophysical Indices','FontSize',15,'FontWeight','bold');
plot (midday,kp,'x-k');
ylabel('Kp','FontSize',fsize,'FontWeight','bold');
xlabel(' ','FontSize',fsize,'FontWeight','bold');
grid;
set (gca,'xlim',[nday1,nday2],'ylim',[0,4.5],'FontSize',fsize,'box','on');
subplot('position',[0.10,0.70,0.87,0.12])	% 2/4 of 1/2
hold on;
plot (midday,bz,'x-k');
plot (xtot,zero,'--k');
ylabel('Bz nT','FontSize',fsize,'FontWeight','bold');
xlabel(' ','FontSize',fsize,'FontWeight','bold');
grid;
set (gca,'xlim',[nday1,nday2],'ylim',[-1.9,1.9],'FontSize',fsize,'box','on');
subplot('position',[0.10,0.58,0.87,0.12])	% 3/4 of 1/2
hold on;
plot (midday,bt,'x-k');
ylabel('Bt nT','FontSize',fsize,'FontWeight','bold');
xlabel(' ','FontSize',fsize,'FontWeight','bold');
grid;
set (gca,'xlim',[nday1,nday2],'ylim',[0,11],'FontSize',fsize,'box','on');
subplot('position',[0.10,0.46,0.87,0.12])	% 4/4 of 1/2
hold on;
plot (midday,vsw,'x-k');
ylabel('Vsw km/s','FontSize',fsize,'FontWeight','bold');
xlabel('2007 Daynumber','FontSize',fsize,'FontWeight','bold');
grid;
set (gca,'xlim',[nday1,nday2],'ylim',[250,750],'FontSize',fsize,'box','on');

tecx = 7.0;
tecx = 15.0;
tecn = 4.5;
yt = tecn + 0.90*(tecx-tecn);
ymit = 5.5;
%subplot('position',[0.10,0.30,0.87,0.10])	% 2/3
subplot('position',[0.10,0.26,0.87,0.14])	% 2/3
hold on;
%plot (midday,gtec,'x-k');
plot (midday,datigs(:,1),'r');
plot (midday,datigs(:,2),'--r');
plot (midday,datigs(:,3),'m');
plot (midday,datigs(:,4),'--m');
plot (midday,datigs(:,5),'g');
plot (midday,datigs(:,6),'--g');
plot (midday,datigs(:,7),'b');
plot (midday,datigs(:,8),'--b');
plot (midday,datmit(:,1),'-.r');	% 72.6% area 8.5TECU
%plot (midday,datmit(:,2),':r');	% 67.5% area bad hi 13.1TECU
plot (midday,datmit(:,3),'-.m');	% 88.4% area 7.8TECU
plot (midday,datmit(:,4),':m');		% 74.1% area 8.1TECU
%plot (midday,datmit(:,5),'-.g');	% 55.4% area 8.4TECU
%plot (midday,datmit(:,6),':g');	% 55.6% area bad lo 6.3TECU
plot (midday,datmit(:,7),'-.b');	% 96.5% area 8.2TECU
%plot (midday,datmit(:,9),':b');	% 42.8% area 7.9TECU lo
text(nday1+0.03*(nday2-nday1),ymit,'MIT __.,..','Color','k','FontSize',fsize,'FontWeight','bold');
text(nday1+0.17*(nday2-nday1),ymit,'73%','Color','r','FontSize',fsize,'FontWeight','bold');
text(nday1+0.24*(nday2-nday1),ymit,'88,74%','Color','m','FontSize',fsize,'FontWeight','bold');
text(nday1+0.36*(nday2-nday1),ymit,'96%','Color','b','FontSize',fsize,'FontWeight','bold');
text(nday1+0.03*(nday2-nday1),yt,'IGS __,--','Color','k','FontSize',fsize,'FontWeight','bold');
text(nday1+0.18*(nday2-nday1),yt,'25,90E','Color','r','FontSize',fsize,'FontWeight','bold');
text(nday1+0.47*(nday2-nday1),yt,'140,175E','Color','m','FontSize',fsize,'FontWeight','bold');
text(nday1+0.62*(nday2-nday1),yt,'200,250E','Color','g','FontSize',fsize,'FontWeight','bold');
text(nday1+0.82*(nday2-nday1),yt,'285,345E','Color','b','FontSize',fsize,'FontWeight','bold');
ylabel('TECU','FontSize',fsize,'FontWeight','bold');
xlabel(' ','FontSize',fsize,'FontWeight','bold');
grid;
set (gca,'xlim',[nday1,nday2],'ylim',[tecn,tecx],'FontSize',fsize,'box','on');
hold off;

y1 = 0.35e-12;
y2 = 1.35e-12;
yt1 = y2 - 0.1*(y2-y1);
yt2 = y2 - 0.23*(y2-y1);
%subplot('position',[0.10,0.050,0.87,0.30])	% 2/2
subplot('position',[0.10,0.050,0.87,0.15])	% 2/3
hold on;
title('Global and CHAMP Neutral Density at 400 km','FontSize',15,'FontWeight','bold');
%plot (midday,drag,'d-k');
plot (day400,drag400,'d-k');
plot (champday,dencad,'*-','Color',dark_green);
plot (luhrday,ldencad,'*-b');
plot(nday1+0.03*(nday2-nday1),yt1,'dk');
text(nday1+0.05*(nday2-nday1),yt1,'3-d global sat drag','Color','k','FontSize',fsize,'FontWeight','bold');
text(nday1+0.41*(nday2-nday1),yt1,'CHAMP av asc(8-3LT)+desc(20-15LT)','Color','k','FontSize',fsize,'FontWeight','bold');
text(nday1+0.41*(nday2-nday1),yt2,'* Doornbos','Color','b','FontSize',fsize,'FontWeight','bold');
text(nday1+0.65*(nday2-nday1),yt2,'* AFRL','Color',dark_green,'FontSize',fsize,'FontWeight','bold');
ylabel('kg/m^3','FontSize',fsize,'FontWeight','bold');
xlabel('2007 Daynumber','FontSize',fsize,'FontWeight','bold');
grid;
set (gca,'xlim',[nday1,nday2],'ylim',[y1,y2],'FontSize',fsize,'box','on');
hold off;

% TEMP
%stop

% Plot <=kp1 and >=kp2 periods together
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.82,0.87,0.12])	% 1/4 of 1/2
hold on;
title('Median Daily Geophysical Indices','FontSize',15,'FontWeight','bold');
plot (midday2,kp2,'*-r');
plot (midday1,kp1,'s--b');
ylabel('Kp','FontSize',fsize,'FontWeight','bold');
xlabel(' ','FontSize',fsize,'FontWeight','bold');
text (nday1+0.9,3.9,'__ Kp>=2 (av=2.79)','Color','r','FontSize',fsize,'FontWeight','bold'); 
text (nday1+0.9,1.5,'---- Kp<=1 (av=0.46)','Color','b','FontSize',fsize,'FontWeight','bold'); 
grid;
set (gca,'xlim',[nday1,nday2],'ylim',[0,4.5],'FontSize',fsize,'box','on');
subplot('position',[0.10,0.70,0.87,0.12])	% 2/4 of 1/2
hold on;
plot (midday2,bz2,'*-r');
plot (midday1,bz1,'s--b');
plot (xtot,zero,'--k');
ylabel('Bz nT','FontSize',fsize,'FontWeight','bold');
xlabel(' ','FontSize',fsize,'FontWeight','bold');
text (nday1+0.9,1.4,'av=-0.08','Color','r','FontSize',fsize,'FontWeight','bold'); 
text (nday1+0.9,-1.6,'av=+0.09','Color','b','FontSize',fsize,'FontWeight','bold'); 
grid;
set (gca,'xlim',[nday1,nday2],'ylim',[-1.9,1.9],'FontSize',fsize,'box','on');
subplot('position',[0.10,0.58,0.87,0.12])	% 3/4 of 1/2
hold on;
plot (midday2,bt2,'*-r');
plot (midday1,bt1,'s--b');
ylabel('Bt nT','FontSize',fsize,'FontWeight','bold');
xlabel(' ','FontSize',fsize,'FontWeight','bold');
text (nday1+0.9,9.8,'av=4.83','Color','r','FontSize',fsize,'FontWeight','bold'); 
text (nday1+0.9,0.9,'av=3.25','Color','b','FontSize',fsize,'FontWeight','bold'); 
grid;
set (gca,'xlim',[nday1,nday2],'ylim',[0,11],'FontSize',fsize,'box','on');
subplot('position',[0.10,0.46,0.87,0.12])	% 4/4 of 1/2
hold on;
plot (midday2,vsw2,'*-r');
plot (midday1,vsw1,'s--b');
ylabel('Vsw km/s','FontSize',fsize,'FontWeight','bold');
xlabel('2007 Daynumber','FontSize',fsize,'FontWeight','bold');
text (nday1+0.9,680,'av=606','Color','r','FontSize',fsize,'FontWeight','bold'); 
text (nday1+0.9,290,'av=359','Color','b','FontSize',fsize,'FontWeight','bold'); 
grid;
set (gca,'xlim',[nday1,nday2],'ylim',[250,750],'FontSize',fsize,'box','on');


y2 = -11.7;
y1 = -12.5;
yt = y2-0.1;
y1 = 0.4e-12;
y2 = 1.6e-12;
yt = 1.53e-12;
subplot('position',[0.10,0.050,0.87,0.30])	% 2/2
hold on;
title('Global Neutral Density at 400 km','FontSize',15,'FontWeight','bold');
plot (midday2,denmod2,'*-k');
plot (midday1,denmod1,'s--k');
plot (midday2,drag2,'*-r');
plot (midday1,drag1,'s--r');
plot (midday2,msis2,'*-c');
plot (midday1,msis1,'s--c');
text(nday1+0.01*(nday2-nday1),yt,'satellite drag (0.759,0.524)','Color','r','FontSize',fsize,'FontWeight','bold');
text(nday1+0.35*(nday2-nday1),yt,'MSIS (1.072,0.862)','Color','c','FontSize',fsize,'FontWeight','bold');
text(nday1+0.60*(nday2-nday1),yt,'TIEGCM-WT (1.090,0.881e-12)','Color','k','FontSize',fsize,'FontWeight','bold');
text(nday1+0.01*(nday2-nday1),1.45e-12,'10.7 cm flux: __ 72.8, ---- 74.3','Color','k','FontSize',fsize,'FontWeight','bold');
ylabel('kg/m3','FontSize',fsize,'FontWeight','bold');
xlabel('2007 Daynumber','FontSize',fsize,'FontWeight','bold');
grid;
set (gca,'xlim',[nday1,nday2],'ylim',[y1,y2],'FontSize',fsize,'box','on');
hold off;

% Plot JULIA and mag Viz from Jicamarca with new and full moon dates
fmoon = 15;
julp = [julpm,julpm];
jula = [julam,julam];
magp = [magpm,magpm];
maga = [magam,magam];

figure
vmin = -10.;
vmax = 30.;
vt1 = 0.95*(vmax-vmin)+vmin;
vt2 = 0.90*(vmax-vmin)+vmin;
vt3 = 0.85*(vmax-vmin)+vmin;
title('Jicamarca Viz','FontSize',15,'FontWeight','bold');
hold on;
plot(midday,datjul(:,3),'dr');
plot(midday,datjul(:,5),'db');
plot(midday,datviz(:,3),'pr');
plot(midday,datviz(:,5),'pb');
%plot(midday,datviz(:,7),'p-g');
plot (midday,dragviz,'x-k');
plot(xtot,julp,'-.b');
plot(xtot,jula,'-.r');
plot(xtot,magp,'--b');
plot(xtot,maga,'--r');
text(328.7,vt3,'f','Color','k','FontSize',fmoon,'FontWeight','bold');
text(357.7,vt3,'f','Color','k','FontSize',fmoon,'FontWeight','bold');
text(342.7,vt3,'n','Color','k','FontSize',fmoon,'FontWeight','bold');
text(372.7,vt3,'n','Color','k','FontSize',fmoon,'FontWeight','bold');
text(nday1+0.01*(nday2-nday1),vt2,'<noon LT','Color','r','FontSize',fsize,'FontWeight','bold');
text(nday1+0.20*(nday2-nday1),vt2,'>noon LT','Color','b','FontSize',fsize,'FontWeight','bold');
text(nday1+0.01*(nday2-nday1),vt1,'x satellite drag','Color','k','FontSize',fsize,'FontWeight','bold');
plot(nday1+0.45*(nday2-nday1),vt1,'d','Color','k');
text(nday1+0.47*(nday2-nday1),vt1,'JULIA','Color','k','FontSize',fsize,'FontWeight','bold');
plot(nday1+0.65*(nday2-nday1),vt1,'p','Color','k');
text(nday1+0.67*(nday2-nday1),vt1,'magnetometer','Color','k','FontSize',fsize,'FontWeight','bold');
ylabel('Viz m/s','FontSize',fsize,'FontWeight','bold');
xlabel('2007 Daynumber','FontSize',fsize,'FontWeight','bold');
grid;
set (gca,'xlim',[nday1,nday2],'ylim',[vmin,vmax],'FontSize',fsize,'box','on');
hold off;

% Plot w/o JULIA data - only mag Viz
figure;
vmin = -5.;
vmax = 20.;
vt1 = 0.95*(vmax-vmin)+vmin;
vt3 = 0.90*(vmax-vmin)+vmin;
title('Jicamarca Viz','FontSize',15,'FontWeight','bold');
hold on;
plot(midday,datviz(:,3),'pr');
plot(midday,datviz(:,5),'pb');
%plot(midday,datviz(:,7),'p-g');
plot (midday,dragviz,'x-k');
plot(xtot,magp,'--b');
plot(xtot,maga,'--r');
text(328.7,vt3,'f','Color','k','FontSize',fmoon,'FontWeight','bold');
text(357.7,vt3,'f','Color','k','FontSize',fmoon,'FontWeight','bold');
text(342.7,vt3,'n','Color','k','FontSize',fmoon,'FontWeight','bold');
text(372.7,vt3,'n','Color','k','FontSize',fmoon,'FontWeight','bold');
text(nday1+0.28*(nday2-nday1),vt1,'<noon LT','Color','r','FontSize',fsize,'FontWeight','bold');
text(nday1+0.45*(nday2-nday1),vt1,'>noon LT','Color','b','FontSize',fsize,'FontWeight','bold');
text(nday1+0.01*(nday2-nday1),vt1,'x satellite drag','Color','k','FontSize',fsize,'FontWeight','bold');
plot(nday1+0.65*(nday2-nday1),vt1,'p','Color','k');
text(nday1+0.67*(nday2-nday1),vt1,'magnetometer','Color','k','FontSize',fsize,'FontWeight','bold');
ylabel('Viz m/s','FontSize',fsize,'FontWeight','bold');
xlabel('2007 Daynumber','FontSize',fsize,'FontWeight','bold');
grid;
set (gca,'xlim',[nday1,nday2],'ylim',[vmin,vmax],'FontSize',fsize,'box','on');
hold off;

if (ifsat==1)
% Revise nday1,2
%nday1 = nday1 - 25;
%nday2 = nday2 + 25;
end

% Plot CHAMP at 400 km
avlden = mean(ldencad)
avljb08 = mean(ljb08cad)
avlhasdm = mean(lhasdmcad)
avlmsis = mean(lmsiscad)
rxnld = max(ldencad)/min(ldencad)
rxnljb = max(ljb08cad)/min(ljb08cad)
rxnlhd = max(lhasdmcad)/min(lhasdmcad)
rxnlm = max(lmsiscad)/min(lmsiscad)
max(ldencad),min(ldencad)
max(ljb08cad),min(ljb08cad)
max(lhasdmcad),min(lhasdmcad)
max(lmsiscad),min(lmsiscad)
dxnld = max(ldencad)-min(ldencad)
dxnljb = max(ljb08cad)-min(ljb08cad)
dxnlhd = max(lhasdmcad)-min(lhasdmcad)
dxnlm = max(lmsiscad)-min(lmsiscad)
nxnld = (max(ldencad)-min(ldencad))/avlden
nxnljb = (max(ljb08cad)-min(ljb08cad))/avljb08
nxnlhd = (max(lhasdmcad)-min(lhasdmcad))/avlhasdm
nxnlm = (max(lmsiscad)-min(lmsiscad))/avlmsis
avden = mean(dencad)
rxnd = max(dencad)/min(dencad)
nxnd = (max(dencad)-min(dencad))/avden
dxnd = max(dencad)-min(dencad)
max(dencad),min(dencad)
ndm = length(midday)
avdrag = 0.;
avmsis = 0.;
nr = 0;
xdr = 0.;
ndr = 100.e-13;
xms = 0.;
nms = 100.e-13;
for nd=1:ndm
 if (midday(nd)>325 & midday(nd)<386)
  nr = nr + 1;
  avdrag = avdrag + drag400(nd);
  avmsis = avmsis + msis400(nd);
  xdr = max(xdr,drag400(nd));
  ndr = min(ndr,drag400(nd));
  xms = max(xms,msis400(nd));
  nms = min(nms,msis400(nd));
 end
end
nr
avdrag = avdrag/nr
avmsis = avmsis/nr
rxndr = xdr/ndr
rxnms = xms/nms
dxndr = xdr-ndr
dxnms = xms-nms
nxndr = (xdr-ndr)/avdrag
nxnms = (xms-nms)/avmsis
xdr,ndr,xms,nms
% Larger than Luehr densities are: drag 25%, U CO 33%, MSIS(L) 84%


y1 = 0.0e-12;
y2 = 1.55e-12;
y2 = 1.75e-12;
if (ifsat==1)
 y2 = 5.5e-12;
end
yt1 = y2 - 0.06*(y2-y1);
yt2 = y2 - 0.12*(y2-y1);
yt3 = y2 - 0.18*(y2-y1);
yt4 = y2 - 0.24*(y2-y1);
yt5 = y2 - 0.30*(y2-y1);
yb1 = y1 + 0.03*(y2-y1);
yb2 = y1 + 0.08*(y2-y1);
yb3 = y1 + 0.13*(y2-y1);

% Plot everything
figure
hold on;
if (ifsat==0)
 title('Global and CHAMP Neutral Density at 400 km','FontSize',15,'FontWeight','bold');
else
 title('Global Neutral Density at CHAMP (337-368km)','FontSize',15,'FontWeight','bold');
end
plot (midday,gdenta,'x-','Color',purple);	% TIME-AMIE
%plot (midday,gden5,'s-c');	% UAM5
plot (midday,gden4,'x-','Color',aqua);	% UAM4 w MSIS
plot (midday,denmod,'x-','Color',light_orange);	% TIEGCM-Weimer-TIMED
plot (midday,gden3,'x-','Color',yellow_green);	% TIEGCM-1.93 Heelis Kp (CCMC)
%plot (midday,drag,'d-r');
 plot (day400,drag400,'d-k');
plot (champday,dencad,'*-','Color',dark_green);
plot (luhrday,ldencad,'*-b');
plot (luhrday,lmsiscad,'s-m');
plot (luhrday,ljb08cad,'s-g');
plot (luhrday,lhasdmcad,'s-c');
%plot (champday,denasc,'p-b');
%plot (champday,dendesc,'p-m');
%plot (luhrday,ldencad,'o-g');
%plot (luhrday,ldenasc,'o-b');
%plot (luhrday,ldendesc,'o-m');
%plot (luhrday,lmsiscad,'s-g');
%plot (luhrday,lmsisasc,'s-b');
%plot (luhrday,lmsisdesc,'s-m');
if (ifsat==1)
 plot (champday,dencad2,'o-g');
 plot (champday,denasc2,'o-b');
 plot (champday,dendesc2,'o-m');
text(nday1+0.03*(nday2-nday1),yt2,'o MSIS00','Color','k','FontSize',fsize,'FontWeight','bold');
text(nday1+0.21*(nday2-nday1),yt2,'x TIE','Color','k','FontSize',fsize,'FontWeight','bold');
 plot (champday,dencadmod,'x-g');
 plot (champday,denascmod,'x-b');
 plot (champday,dendescmod,'x-m');
text(nday1+0.03*(nday2-nday1),yt3,'+ MSIS/CHAMP, x TIE/CHAMP','Color','k','FontSize',fsize,'FontWeight','bold');
 plot (champday,rmdasc*1e-12,'+-k');
 plot (champday,rmddesc*1e-12,'+-.k');
plot (champday,rmoddasc*1e-12,'x-k');
plot (champday,rmodddesc*1e-12,'x-.k');
else
plot (day400,msis400,'o-r');
text(nday1+0.03*(nday2-nday1),yt2,'o global MSIS00','Color','r','FontSize',fsize,'FontWeight','bold');
text(nday1+0.41*(nday2-nday1),yt2,'* Doornbos','Color','b','FontSize',fsize,'FontWeight','bold');
text(nday1+0.65*(nday2-nday1),yt2,'* AFRL','Color',dark_green,'FontSize',fsize,'FontWeight','bold');
plot (nday1+0.41*(nday2-nday1),yt3,'sm');
text(nday1+0.45*(nday2-nday1),yt3,'MSIS00','Color','m','FontSize',fsize,'FontWeight','bold');
plot (nday1+0.61*(nday2-nday1),yt3,'sg');
text(nday1+0.65*(nday2-nday1),yt3,'JB2008','Color','g','FontSize',fsize,'FontWeight','bold');
plot (nday1+0.81*(nday2-nday1),yt3,'sc');
text(nday1+0.85*(nday2-nday1),yt3,'HASDM','Color','c','FontSize',fsize,'FontWeight','bold');
text (nday1+0.03*(nday2-nday1),yt3,'x UAM4','Color',aqua,'FontSize',fsize,'FontWeight','bold');
text (nday1+0.17*(nday2-nday1),yt4,'x TIEWT','Color',light_orange,'FontSize',fsize,'FontWeight','bold');
text (nday1+0.17*(nday2-nday1),yt3,'x TIME-AMIE','Color',purple,'FontSize',fsize,'FontWeight','bold');
text (nday1+0.03*(nday2-nday1),yt4,'x TIEKp','Color',yellow_green,'FontSize',fsize,'FontWeight','bold');
text(nday1+0.01*(nday2-nday1),yb3,'ave*e-13: 5.4D 6.4drag 6.8AFRL 7.0H 7.9JB 9.4DM 9.8MSIS','Color','k','FontSize',fsize,'FontWeight','bold');
text(nday1+0.01*(nday2-nday1),yb1,'(x-n)/av: .96D .81drag .84AFRL 1.01H .87JB .56DM .47MSIS','Color','k','FontSize',fsize,'FontWeight','bold');
text(nday1+0.01*(nday2-nday1),yb2,'x-n*e-13: 5.2D 5.2drag 5.7AFRL 7.1H 7.0JB 5.3DM 4.6MSIS','Color','k','FontSize',fsize,'FontWeight','bold');
%text(nday1+0.01*(nday2-nday1),yb1,'max/min: 2.4D 2.2drag 2.2AFRL 2.5H 2.3JB 1.7DM 1.6MSIS','Color','k','FontSize',fsize,'FontWeight','bold');
end
plot(nday1+0.03*(nday2-nday1),yt1,'dk');
text(nday1+0.05*(nday2-nday1),yt1,'global sat drag','Color','k','FontSize',fsize,'FontWeight','bold');
text(nday1+0.41*(nday2-nday1),yt1,'CHAMP av asc(8-3LT)+desc(20-15LT)','Color','k','FontSize',fsize,'FontWeight','bold');
%text(nday1+0.48*(nday2-nday1),yt1,'descending','Color','m','FontSize',fsize,'FontWeight','bold');
%text(nday1+0.48*(nday2-nday1),yt2,'20-15LT','Color','m','FontSize',fsize,'FontWeight','bold');
%text(nday1+0.71*(nday2-nday1),yt1,'ascending','Color','b','FontSize',fsize,'FontWeight','bold');
%text(nday1+0.71*(nday2-nday1),yt2,'8-3LT','Color','b','FontSize',fsize,'FontWeight','bold');
%text(nday1+0.90*(nday2-nday1),yt1,'both','Color','g','FontSize',fsize,'FontWeight','bold');
ylabel('kg/m3','FontSize',fsize,'FontWeight','bold');
xlabel('2007 Daynumber','FontSize',fsize,'FontWeight','bold');
grid;
set (gca,'xlim',[nday1,nday2],'ylim',[y1,y2],'FontSize',fsize,'box','on');
hold off;

% Plot global and CHAMP obs and empirical models only
figure
hold on;
 title('Global and CHAMP Neutral Density at 400 km','FontSize',15,'FontWeight','bold');
plot (day400,drag400,'d-k');
plot (champday,dencad,'*-','Color',dark_green);
plot (luhrday,ldencad,'*-b');
plot (luhrday,lmsiscad,'s-m');
plot (luhrday,ljb08cad,'s-g');
plot (luhrday,lhasdmcad,'s-c');
plot (day400,msis400,'o-r');
text(nday1+0.03*(nday2-nday1),yt2,'o global MSIS00','Color','r','FontSize',fsize,'FontWeight','bold');
text(nday1+0.41*(nday2-nday1),yt2,'* Doornbos','Color','b','FontSize',fsize,'FontWeight','bold');
text(nday1+0.65*(nday2-nday1),yt2,'* AFRL','Color',dark_green,'FontSize',fsize,'FontWeight','bold');
plot (nday1+0.41*(nday2-nday1),yt3,'sm');
text(nday1+0.45*(nday2-nday1),yt3,'MSIS00','Color','m','FontSize',fsize,'FontWeight','bold');
plot (nday1+0.61*(nday2-nday1),yt3,'sg');
text(nday1+0.65*(nday2-nday1),yt3,'JB2008','Color','g','FontSize',fsize,'FontWeight','bold');
plot (nday1+0.81*(nday2-nday1),yt3,'sc');
text(nday1+0.85*(nday2-nday1),yt3,'HASDM','Color','c','FontSize',fsize,'FontWeight','bold');
plot(nday1+0.03*(nday2-nday1),yt1,'dk');
text(nday1+0.05*(nday2-nday1),yt1,'global sat drag','Color','k','FontSize',fsize,'FontWeight','bold');
text(nday1+0.41*(nday2-nday1),yt1,'CHAMP av asc(8-3LT)+desc(20-15LT)','Color','k','FontSize',fsize,'FontWeight','bold');
ylabel('kg/m3','FontSize',fsize,'FontWeight','bold');
xlabel('2007 Daynumber','FontSize',fsize,'FontWeight','bold');
grid;
set (gca,'xlim',[nday1,nday2],'ylim',[y1,y2],'FontSize',fsize,'box','on');
hold off;

% Plot global model comparisons
figure
hold on;
 title('Global Neutral Density at 400 km','FontSize',15,'FontWeight','bold');
plot (midday,gdenta,'x-','Color',purple);	% TIME-AMIE
%plot (midday,gden5,'s-c');	% UAM5
plot (midday,gden4,'x-','Color',aqua);	% UAM4 w MSIS
plot (midday,denmod,'x-m');	% TIEGCM-Weimer-TIMED
plot (midday,gden3,'x-','Color',yellow_green);	% TIEGCM-1.93 Heelis Kp (CCMC)
% 61 midday 325.5-385.5
% 300 day400 300.5-599.50 (so want 61 days from #26-86)
 plot (day400,drag400,'d-k');
% 58 champday = 325.5-385.5 missing 371-373
%plot (champday,dencad,'*-','Color',dark_green);
% 60 luhrday = 325.5-385.5 missing 372
%plot (luhrday,ldencad,'*-b');
plot (day400,msis400,'o-r');
plot(nday1+0.03*(nday2-nday1),yt1,'dk');
text(nday1+0.05*(nday2-nday1),yt1,'global sat drag','Color','k','FontSize',fsize,'FontWeight','bold');
text(nday1+0.03*(nday2-nday1),yt2,'o global MSIS00','Color','r','FontSize',fsize,'FontWeight','bold');
%text(nday1+0.41*(nday2-nday1),yt1,'CHAMP av asc(8-3LT)+desc(20-15LT)','Color','k','FontSize',fsize,'FontWeight','bold');
%text(nday1+0.41*(nday2-nday1),yt2,'* Doornbos','Color','b','FontSize',fsize,'FontWeight','bold');
%text(nday1+0.70*(nday2-nday1),yt2,'* AFRL','Color',dark_green,'FontSize',fsize,'FontWeight','bold');
text (nday1+0.03*(nday2-nday1),yb1,'x UAM4','Color',aqua,'FontSize',fsize,'FontWeight','bold');
text (nday1+0.40*(nday2-nday1),yb1,'x TIEWT','Color','m','FontSize',fsize,'FontWeight','bold');
text (nday1+0.60*(nday2-nday1),yb1,'x TIME-AMIE','Color',purple,'FontSize',fsize,'FontWeight','bold');
text (nday1+0.22*(nday2-nday1),yb1,'x TIEKp','Color',yellow_green,'FontSize',fsize,'FontWeight','bold');
ylabel('kg/m3','FontSize',fsize,'FontWeight','bold');
xlabel('2007 Daynumber','FontSize',fsize,'FontWeight','bold');
grid;
set (gca,'xlim',[nday1,nday2],'ylim',[y1,y2],'FontSize',fsize,'box','on');
hold off;

% Plot CHAMP model comparisons
av58afrl=mean(dencad)
av60doorn=mean(ldencad)
av58doorn=(mean(ldencad)*60-ldencad(47)-ldencad(48))/58.
av61uam4=mean(gden4)
av61tiekp=mean(gden3)
av61tiewt=mean(denmod)
av61tima=mean(gdenta)
av58uam4=(mean(gden4)*61-gden4(47)-gden4(48)-gden4(49))/58.
av58tiekp=(mean(gden3)*61-gden3(47)-gden3(48)-gden3(49))/58.
av58tiewt=(mean(denmod)*61-denmod(47)-denmod(48)-denmod(49))/58.
av58tima=(mean(gdenta)*61-gdenta(47)-gdenta(48)-gdenta(49))/58.
% Calculate 58-d and 61-d means,max,min of 61 days 
% Will subtract means from the global (and CHAMP obs) 
av61drag = 0.;
av61msis = 0.;
av58drag = 0.;
av58msis = 0.;
max61drag = 0.;
max58drag = 0.;
mindrag = 1.;
maxmsis = 0.;
minmsis = 1.;
max58doorn = 0.;
for it=26:86
 nd = it - 25;
 av61drag = av61drag + drag400(it);
 av61msis = av61msis + msis400(it);
 if (max61drag<drag400(it))
   max61drag = drag400(it);
 end
 if (mindrag>drag400(it))
   mindrag = drag400(it);
 end
 if (maxmsis<msis400(it))
   maxmsis = msis400(it);
 end
 if (minmsis>msis400(it))
   minmsis = msis400(it);
 end
 if (day400(it)<371. | day400(it)>374.)
  av58drag = av58drag + drag400(it);
  av58msis = av58msis + msis400(it);
  if (max58drag<drag400(it))
   max58drag = drag400(it);
  end
  if (nd<61)
   if (max58doorn<ldencad(nd))
    max58doorn = ldencad(nd);
   end
  end 
 end
end
av61drag=av61drag/61.
av61msis=av61msis/61.
av58drag=av58drag/58.
av58msis=av58msis/58.
max61drag
max58drag
mindrag
maxmsis
minmsis
xn61drag = max61drag-mindrag
xn58drag = max58drag-mindrag
xnmsis = maxmsis-minmsis
xnafrl = max(dencad) - min(dencad)
xn60doorn = max(ldencad) - min(ldencad)
xn58doorn = max58doorn - min(ldencad)
xntiekp = max(gden3) - min(gden3)
xntiewt = max(denmod) - min(denmod)
xnuam4 = max(gden4) - min(gden4)
xntima = max(gdenta) - min(gdenta)

dragav58=[av58drag av58drag];
msisav58=[av58msis av58msis];
doornav58=[av58doorn av58doorn];
afrlav58=[av58afrl av58afrl];
timaav58=[av58tima av58tima];
kpav58=[av58tiekp av58tiekp];
wtav58=[av58tiewt av58tiewt];
uam4av58=[av58uam4 av58uam4];
dragxn=[xn58drag xn58drag];
msisxn=[xnmsis xnmsis];
doornxn=[xn58doorn xn58doorn];
afrlxn=[xnafrl xnafrl];
timaxn=[xntima xntima];
kpxn=[xntiekp xntiekp];
wtxn=[xntiewt xntiewt];
uam4xn=[xnuam4 xnuam4];
% Plot averages and ranges (max-min)
y2 = 1.3e-12;
y1 = 0.;
yt1 = y2 - 0.06*(y2-y1);
yt2 = y2 - 0.12*(y2-y1);
yt3 = y2 - 0.18*(y2-y1);
yt4 = y2 - 0.24*(y2-y1);
yb1 = y1 + 0.03*(y2-y1);
yb2 = y1 + 0.08*(y2-y1);
yb3 = y1 + 0.13*(y2-y1);
figure
hold on;
 title('Av Global Neutral Density and --(max-min)','FontSize',14,'FontWeight','bold');
plot (xtot,dragav58,'d-k');
plot (xtot,dragxn,'d--k');
plot (xtot,msisav58,'o-r');
plot (xtot,msisxn,'o--r');
plot (xtot,doornav58,'*-b');
plot (xtot,doornxn,'*--b');
plot (xtot,afrlav58,'*-','Color',dark_green);
plot (xtot,afrlxn,'*--','Color',dark_green);
plot (xtot,uam4av58,'-','Color',aqua);
plot (xtot,uam4xn,'--','Color',aqua);
plot (xtot,timaav58,'-','Color',purple);
plot (xtot,timaxn,'--','Color',purple);
plot (xtot,kpav58,'-','Color',yellow_green);
plot (xtot,kpxn,'--','Color',yellow_green);
plot (xtot,wtav58,'-m');
plot (xtot,wtxn,'--m');
plot(nday1+0.06*(nday2-nday1),yt1,'dk');
text(nday1+0.08*(nday2-nday1),yt1,'global sat drag 0-0','Color','k','FontSize',fsize,'FontWeight','bold');
text(nday1+0.06*(nday2-nday1),yt2,'o global MSIS00 6-2','Color','r','FontSize',fsize,'FontWeight','bold');
text(nday1+0.41*(nday2-nday1),yt1,'CHAMP av asc(8-3LT)+desc(20-15LT)','Color','k','FontSize',fsize,'FontWeight','bold');
text(nday1+0.41*(nday2-nday1),yt2,'* Doornbos 3-4','Color','b','FontSize',fsize,'FontWeight','bold');
text(nday1+0.70*(nday2-nday1),yt2,'* AFRL 1-3','Color',dark_green,'FontSize',fsize,'FontWeight','bold');
text (nday1+0.03*(nday2-nday1),yb1,'x UAM4 2-7','Color',aqua,'FontSize',fsize,'FontWeight','bold');
text (nday1+0.40*(nday2-nday1),yb1,'x TIEWT 7-5','Color','m','FontSize',fsize,'FontWeight','bold');
text (nday1+0.60*(nday2-nday1),yb1,'x TIME-AMIE 4-6','Color',purple,'FontSize',fsize,'FontWeight','bold');
text (nday1+0.22*(nday2-nday1),yb1,'x TIEKp 5-1','Color',yellow_green,'FontSize',fsize,'FontWeight','bold');
ylabel('kg/m3','FontSize',fsize,'FontWeight','bold');
xlabel('2007 Daynumber','FontSize',fsize,'FontWeight','bold');
grid;
set (gca,'xlim',[nday1,nday2],'ylim',[y1,y2],'FontSize',fsize,'box','on');
hold off;

% Find 3-d av to cf w drag (and MSIS)
for nd=1:61
 na = nd;
 if (nd>46)
  na = nd - 3;
 end
 nl = nd;
 if (nd>47)
  nl = nd - 1;
 end
 if (nd==1)
  den3afrl(na) = 0.5*(dencad(na)+dencad(na+1));
  den3doorn(nl) = 0.5*(ldencad(nl)+ldencad(nl+1));
  den3msisl(nl) = 0.5*(lmsiscad(nl)+lmsiscad(nl+1));
  den3kp(nd) = 0.5*(gden3(nd)+gden3(nd+1));
  den3uam(nd) = 0.5*(gden4(nd)+gden4(nd+1));
  den3wt(nd) = 0.5*(denmod(nd)+denmod(nd+1));
  den3ta(nd) = 0.5*(gdenta(nd)+gdenta(nd+1));
 end
 if (nd==61)
  den3afrl(na) = 0.5*(dencad(na)+dencad(na-1));
  den3doorn(nl) = 0.5*(ldencad(nl)+ldencad(nl-1));
  den3msisl(nl) = 0.5*(lmsiscad(nl)+lmsiscad(nl-1));
  den3kp(nd) = 0.5*(gden3(nd)+gden3(nd-1));
  den3uam(nd) = 0.5*(gden4(nd)+gden4(nd-1));
  den3wt(nd) = 0.5*(denmod(nd)+denmod(nd-1));
  den3ta(nd) = 0.5*(gdenta(nd)+gdenta(nd-1));
 end
 if (nd>1 & nd<61)
  if (nd<46 | nd>50)
   den3afrl(na) = (dencad(na-1)+dencad(na)+dencad(na+1))/3.;
  else
  if (nd==46)
   den3afrl(na) = 0.5*(dencad(na)+dencad(na-1));
  end
  if (nd==50)
   den3afrl(na) = 0.5*(dencad(na)+dencad(na+1));
  end
  end
  if (nd<47 | nd>49)
   den3doorn(nl) = (ldencad(nl-1)+ldencad(nl)+ldencad(nl+1))/3.;
   den3msisl(nl) = (lmsiscad(nl-1)+lmsiscad(nl)+lmsiscad(nl+1))/3.;
  else
   if (nd==47)
    den3doorn(nl) = 0.5*(ldencad(nl)+ldencad(nl-1));
    den3msisl(nl) = 0.5*(lmsiscad(nl)+lmsiscad(nl-1));
   end
   if (nd==49)
    den3doorn(nl) = 0.5*(ldencad(nl)+ldencad(nl+1));
    den3msisl(nl) = 0.5*(lmsiscad(nl)+lmsiscad(nl+1));
   end
  end
  den3kp(nd) = (gden3(nd-1)+gden3(nd)+gden3(nd+1))/3.;
  den3uam(nd) = (gden4(nd-1)+gden4(nd)+gden4(nd+1))/3.;
  den3wt(nd) = (denmod(nd-1)+denmod(nd)+denmod(nd+1))/3.;
  den3ta(nd) = (gdenta(nd-1)+gdenta(nd)+gdenta(nd+1))/3.;
 end
end	% for nd=1:61

% Plot global model comparisons
y2 = 4.5e-13;
y1 = -3.5e-13;
yt1 = y2 - 0.06*(y2-y1);
yt2 = y2 - 0.12*(y2-y1);
yt3 = y2 - 0.18*(y2-y1);
yt4 = y2 - 0.24*(y2-y1);
yb1 = y1 + 0.03*(y2-y1);
yb2 = y1 + 0.08*(y2-y1);
yb3 = y1 + 0.13*(y2-y1);
figure
hold on;
 title('3-d Av Global Neutral Density Minus 58-d Av','FontSize',14,'FontWeight','bold');
plot (xtot,zero,'--k');
plot (midday,den3ta-av58tima,'x-','Color',purple);	% TIME-AMIE
%plot (midday,gden5,'s-c');	% UAM5
plot (midday,den3uam-av58uam4,'x-','Color',aqua);	% UAM4 w MSIS
plot (midday,den3wt-av58tiewt,'x-m');	% TIEGCM-Weimer-TIMED
plot (midday,den3kp-av58tiekp,'x-','Color',yellow_green);	% TIEGCM-1.93 Heelis Kp (CCMC)
% 61 midday 325.5-385.5
% 300 day400 300.5-599.50 (so want 61 days from #26-86)
 plot (day400,drag400-av58drag,'d-k');
% 58 champday = 325.5-385.5 missing 371-373
plot (champday,den3afrl-av58afrl,'*-','Color',dark_green);
% 60 luhrday = 325.5-385.5 missing 372
plot (luhrday,den3doorn-av58doorn,'*-b');
plot (day400,msis400-av58msis,'o-r');
plot(nday1+0.06*(nday2-nday1),yt1,'dk');
text(nday1+0.08*(nday2-nday1),yt1,'global sat drag 0','Color','k','FontSize',fsize,'FontWeight','bold');
text(nday1+0.05*(nday2-nday1),yt2,'o global MSIS00 1','Color','r','FontSize',fsize,'FontWeight','bold');
text(nday1+0.41*(nday2-nday1),yt1,'CHAMP av asc(8-3LT)+desc(20-15LT)','Color','k','FontSize',fsize,'FontWeight','bold');
text(nday1+0.41*(nday2-nday1),yt2,'* Doornbos 3','Color','b','FontSize',fsize,'FontWeight','bold');
text(nday1+0.70*(nday2-nday1),yt2,'* AFRL 5','Color',dark_green,'FontSize',fsize,'FontWeight','bold');
text (nday1+0.03*(nday2-nday1),yb1,'x UAM4 4','Color',aqua,'FontSize',fsize,'FontWeight','bold');
text (nday1+0.40*(nday2-nday1),yb1,'x TIEWT 7','Color','m','FontSize',fsize,'FontWeight','bold');
text (nday1+0.60*(nday2-nday1),yb1,'x TIME-AMIE 2','Color',purple,'FontSize',fsize,'FontWeight','bold');
text (nday1+0.22*(nday2-nday1),yb1,'x TIEKp 6','Color',yellow_green,'FontSize',fsize,'FontWeight','bold');
ylabel('kg/m3','FontSize',fsize,'FontWeight','bold');
xlabel('2007 Daynumber','FontSize',fsize,'FontWeight','bold');
grid;
set (gca,'xlim',[nday1,nday2],'ylim',[y1,y2],'FontSize',fsize,'box','on');
hold off;

% Find the Root-Mean-Square RMS=sqrt[sum(xobs-xmod)^2/N] 
rmsmsis=0.;
rmsdoorn=0.;
rmsafrl=0.;
rmstima=0.;
rmsuam4=0.;
rmskp=0.;
rmswt=0.;
for nd=1:61
 it = nd+25;
 rmsmsis = rmsmsis + (drag400(it)-av58drag-msis400(it)+av58msis)^2;
 rmstima = rmstima + (drag400(it)-av58drag-den3ta(nd))^2;
 rmsuam4 = rmsuam4 + (drag400(it)-av58drag-den3uam(nd))^2;
 rmskp = rmskp + (drag400(it)-av58drag-den3kp(nd))^2;
 rmswt = rmswt + (drag400(it)-av58drag-den3wt(nd))^2;
 na = nd;
 if (nd>46)
  na = nd - 3;
 end
 nl = nd;
 if (nd>47)
  nl = nd - 1;
 end
 if (nd<47 | nd>49)
 rmsdoorn = rmsdoorn + (drag400(it)-av58drag-den3doorn(nl))^2;
 rmsafrl = rmsafrl + (drag400(it)-av58drag-den3afrl(na))^2;
 else
  if (nd~=48)
   rmsdoorn = rmsdoorn + (drag400(it)-av58drag-den3doorn(nl))^2;
  end
 end
end
 rmsmsis = sqrt(rmsmsis/61.)
 rmstima = sqrt(rmstima/61.)
 rmsuam4 = sqrt(rmsuam4/61.)
 rmskp = sqrt(rmskp/61.)
 rmswt = sqrt(rmswt/61.)
 rmsdoorn = sqrt(rmsdoorn/60.)
 rmsafrl = sqrt(rmsafrl/58.)


% Read TEC from SAMI3 (global given by Sarah McDonald, and calculated from 8 glons)
filsam0 = sprintf('/hao/aim3/emery/cetic_climatology/sami3_output/sami_globaltec_07325_08020.txt');
datsam0 = dlmread(filsam0,'',1,0);
for nd=1:61
 gtecsam0(nd) = datsam0(nd,4);
end
% Open 8 glon 'global' TEC daily values from SAMI3 TEC
for j=1:8
filesam = sprintf('/hao/aim3/emery/cetic_climatology/sami3_output/global_tec_sami3_07325_08020_lon%3.3d.dat',loncutsE(j));
datsamtnh = dlmread(filesam,'',0,0);
for nd=1:61
 datsam(nd,j) = datsamtnh(nd,5);
end
end		% j=1:8

% Average IGS and SAMI3 over 8 glons and cf w true global TEC
for nd=1:61
 sumsam(nd) = 0.;
 sumigs(nd) = 0.;
 for j=1:8
  sumsam(nd) = sumsam(nd)+datsam(nd,j);
  sumigs(nd) = sumigs(nd)+datigs(nd,j);
 end
 sumsam(nd) = sumsam(nd)/8.;
 sumigs(nd) = sumigs(nd)/8.;
end
% Average MIT over 4 glons only (1,3,4,7) , with 140(3) and 175(4) combined as 1
for nd=1:61
 summit(nd) = (datmit(nd,1)+0.5*(datmit(nd,3)+datmit(nd,4))+datmit(nd,7))/3.;
end

figure
hold on;
tecx = 20.;
tecn = 0.;
yt = tecn + 0.95*(tecx-tecn);
ytie = 1.5;
yigs = 12.5;
 title('Global TEC 07325-08020','FontSize',15,'FontWeight','bold');
plot (midday,gtecta,'x-','Color',purple);
plot (midday,gtec3,'x-','Color',yellow_green);
plot (midday,gtec4,'*-','Color',aqua);
plot (midday,gtec4l,'x-','Color',aqua);
plot (midday,gtec,'x-m');
plot (midday,gtecsam0,'*-','Color',light_orange);
plot (midday,sumsam,'o-','Color',light_orange);
plot (midday,sumigs,'d-k','LineWidth',2);
plot (midday,summit,'d-r','LineWidth',2);
text(nday1+0.03*(nday2-nday1),0.70*(tecx-tecn)+tecn,'IGS 8 lons','Color','k','FontSize',fsize,'FontWeight','bold');
text(nday1+0.25*(nday2-nday1),0.70*(tecx-tecn)+tecn,'MIT 4 lons','Color','r','FontSize',fsize,'FontWeight','bold');
text(nday1+0.03*(nday2-nday1),0.95*(tecx-tecn)+tecn,'* UAM4 gps hts','Color',aqua,'FontSize',fsize,'FontWeight','bold');
text(nday1+0.30*(nday2-nday1),0.95*(tecx-tecn)+tecn,'* SAMI3 gps hts and 8 glons (o)','Color',light_orange,'FontSize',fsize,'FontWeight','bold');
text(nday1+0.03*(nday2-nday1),0.03*(tecx-tecn)+tecn,'x UAM4 500km','Color',aqua,'FontSize',fsize,'FontWeight','bold');
text (nday1+0.50*(nday2-nday1),0.03*(tecx-tecn)+tecn,'x TIEWT','Color','m','FontSize',fsize,'FontWeight','bold');
text (nday1+0.70*(nday2-nday1),0.03*(tecx-tecn)+tecn,'x TIME-AMIE','Color',purple,'FontSize',fsize,'FontWeight','bold');
text (nday1+0.30*(nday2-nday1),0.03*(tecx-tecn)+tecn,'x TIEKp','Color',yellow_green,'FontSize',fsize,'FontWeight','bold');
ylabel('TECU','FontSize',fsize,'FontWeight','bold');
xlabel('2007 Daynumber','FontSize',fsize,'FontWeight','bold');
grid;
set (gca,'xlim',[nday1,nday2],'ylim',[tecn,tecx],'FontSize',fsize,'box','on');
hold off;

% Delete first day of UAM4
for nd=2:61
middaym1(nd-1) = midday(nd);
gtec4m1(nd-1)= gtec4(nd);
gtec4lm1(nd-1) = gtec4l(nd);
end


% Calc metrics of baseline, max-min, RMS
avigs=mean(sumigs)
avmit=mean(summit)
avsami=mean(gtecsam0)
avuam4=mean(gtec4)
avuam4l=mean(gtec4l)
avuam4m1=mean(gtec4m1)
avuam4lm1=mean(gtec4lm1)
avwt=mean(gtec)
avkp=mean(gtec3)
avtima=mean(gtecta)
xnigs=max(sumigs)-min(sumigs)
xnmit=max(summit)-min(summit)
xnsami=max(gtecsam0)-min(gtecsam0)
xnuam4=max(gtec4)-min(gtec4)
xnuam4l=max(gtec4l)-min(gtec4l)
xnuam4m1=max(gtec4m1)-min(gtec4m1)
xnuam4lm1=max(gtec4lm1)-min(gtec4lm1)
xntima=max(gtecta)-min(gtecta)
xnwt=max(gtec)-min(gtec)
xnkp=max(gtec3)-min(gtec3)
% Find the Root-Mean-Square RMS=sqrt[sum(xobs-xmod)^2/N] 
rmsigs=0.;
rmsmit=0.;
rmsigssami=0.;
rmsmitsami=0.;
rmsigsuam4=0.;
rmsmituam4=0.;
rmsigsuam4l=0.;
rmsmituam4l=0.;
rmsigsuam4m1=0.;
rmsmituam4m1=0.;
rmsigsuam4lm1=0.;
rmsmituam4lm1=0.;
rmsigswt=0.;
rmsmitwt=0.;
rmsigskp=0.;
rmsmitkp=0.;
rmsigstima=0.;
rmsmittima=0.;
for it=1:61
 rmsigs = rmsigs + (summit(it)-avmit-sumigs(it)+avigs)^2;
 rmsmit = rmsmit + (sumigs(it)-avigs-summit(it)+avmit)^2;
 rmsmitsami = rmsmitsami + (summit(it)-avmit-gtecsam0(it)+avsami)^2;
 rmsigssami = rmsigssami + (sumigs(it)-avigs-gtecsam0(it)+avsami)^2;
 rmsmituam4 = rmsmituam4 + (summit(it)-avmit-gtec4(it)+avuam4)^2;
 rmsigsuam4 = rmsigsuam4 + (sumigs(it)-avigs-gtec4(it)+avuam4)^2;
 rmsmituam4l = rmsmituam4l + (summit(it)-avmit-gtec4l(it)+avuam4l)^2;
 rmsigsuam4l = rmsigsuam4l + (sumigs(it)-avigs-gtec4l(it)+avuam4l)^2;
 if (it>1)
 rmsmituam4m1 = rmsmituam4m1 + (summit(it)-avmit-gtec4m1(it-1)+avuam4m1)^2;
 rmsigsuam4m1 = rmsigsuam4m1 + (sumigs(it)-avigs-gtec4m1(it-1)+avuam4m1)^2;
 rmsmituam4lm1 = rmsmituam4lm1 + (summit(it)-avmit-gtec4lm1(it-1)+avuam4lm1)^2;
 rmsigsuam4lm1 = rmsigsuam4lm1 + (sumigs(it)-avigs-gtec4lm1(it-1)+avuam4lm1)^2;
 end
 rmsmittima = rmsmittima + (summit(it)-avmit-gtecta(it)+avtima)^2;
 rmsigstima = rmsigstima + (sumigs(it)-avigs-gtecta(it)+avtima)^2;
 rmsmitwt = rmsmitwt + (summit(it)-avmit-gtec(it)+avwt)^2;
 rmsigswt = rmsigswt + (sumigs(it)-avigs-gtec(it)+avwt)^2;
 rmsmitkp = rmsmitkp + (summit(it)-avmit-gtec3(it)+avkp)^2;
 rmsigskp = rmsigskp + (sumigs(it)-avigs-gtec3(it)+avkp)^2;
end
rmsigs=rmsigs/61.
rmsmit=rmsmit/61.
rmsmitsami=rmsmitsami/61.
rmsigssami=rmsigssami/61.
rmsmituam4=rmsmituam4/61.
rmsigsuam4=rmsigsuam4/61.
rmsmituam4l=rmsmituam4l/61.
rmsigsuam4l=rmsigsuam4l/61.
rmsmituam4m1=rmsmituam4m1/60.
rmsigsuam4m1=rmsigsuam4m1/60.
rmsmituam4lm1=rmsmituam4lm1/60.
rmsigsuam4lm1=rmsigsuam4lm1/60.
rmsmittima=rmsmittima/61.
rmsigstima=rmsigstima/61.
rmsmitkp=rmsmitkp/61.
rmsigskp=rmsigskp/61.
rmsmitwt=rmsmitwt/61.
rmsigswt=rmsigswt/61.

tecx=4.;
tecn=-4.;
figure
hold on;
title('Global TEC-Av 07325-08020','FontSize',15,'FontWeight','bold');
plot (xtot,zero,'--k');
plot (midday,gtecta-avtima,'x-','Color',purple);
plot (midday,gtec3-avkp,'x-','Color',yellow_green);
plot (middaym1,gtec4m1-avuam4m1,'*-','Color',aqua);
plot (middaym1,gtec4lm1-avuam4lm1,'x-','Color',aqua);
plot (midday,gtec-avwt,'x-m');
plot (midday,gtecsam0-avsami,'*-','Color',light_orange);
plot (midday,sumigs-avigs,'d-k','LineWidth',2);
plot (midday,summit-avmit,'d-r','LineWidth',2);
text(nday1+0.03*(nday2-nday1),0.95*(tecx-tecn)+tecn,'IGS 8 lons 3,3,1(MIT)','Color','k','FontSize',fsize,'FontWeight','bold');
text(nday1+0.45*(nday2-nday1),0.95*(tecx-tecn)+tecn,'MIT 4 lons 1,3,1(IGS)','Color','r','FontSize',fsize,'FontWeight','bold');
text(nday1+0.03*(nday2-nday1),0.28*(tecx-tecn)+tecn,'* UAM4 gps hts 3,1,6(IGS) 1,2,6(MIT)','Color',aqua,'FontSize',fsize,'FontWeight','bold');
text(nday1+0.03*(nday2-nday1),0.22*(tecx-tecn)+tecn,'* SAMI3 gps hts 2,2,7(IGS) 7,6,7(MIT)','Color',light_orange,'FontSize',fsize,'FontWeight','bold');
text(nday1+0.03*(nday2-nday1),0.03*(tecx-tecn)+tecn,'x UAM4 500km','Color',aqua,'FontSize',fsize,'FontWeight','bold');
text (nday1+0.25*(nday2-nday1),0.15*(tecx-tecn)+tecn,'MIT: 6,4,4     2,1,5      4,5,3     5,7,2','Color','r','FontSize',fsize,'FontWeight','bold');
text (nday1+0.25*(nday2-nday1),0.09*(tecx-tecn)+tecn,'IGS: 7,5,4     4,4,2      5,6,3     6,7,5','Color','k','FontSize',fsize,'FontWeight','bold');
text (nday1+0.50*(nday2-nday1),0.03*(tecx-tecn)+tecn,'x TIEWT','Color','m','FontSize',fsize,'FontWeight','bold');
text (nday1+0.70*(nday2-nday1),0.03*(tecx-tecn)+tecn,'x TIME-AMIE','Color',purple,'FontSize',fsize,'FontWeight','bold');
text (nday1+0.30*(nday2-nday1),0.03*(tecx-tecn)+tecn,'x TIEKp','Color',yellow_green,'FontSize',fsize,'FontWeight','bold');
ylabel('TECU-Av','FontSize',fsize,'FontWeight','bold');
xlabel('2007 Daynumber','FontSize',fsize,'FontWeight','bold');
grid;
set (gca,'xlim',[nday1,nday2],'ylim',[tecn,tecx],'FontSize',fsize,'box','on');
hold off;

stop

tecn=1;
tecx=20;
figure
hold on;
tecx = 20.;
yt = tecn + 0.95*(tecx-tecn);
ytie = 1.5;
yigs = 12.5;
 title('Global and 8 glon TEC 07325-08020','FontSize',15,'FontWeight','bold');
plot (midday,gtecta,'x-b');
plot (midday,gtec3,'x-g');
plot (midday,gtec4,'sc');
plot (midday,gtec4l,'s-c');
%plot (midday,gtec4shp,'s-g');
%plot (midday,gtec4eqb,'s-m');
%plot (midday,gtec4nhp,'s-b');
plot (midday,gtec,'x-k');
plot (midday,gtecsam0,'*-k');
plot (midday,datsam(:,1),'or');
plot (midday,datsam(:,2),'o--r');
plot (midday,datsam(:,3),'om');
plot (midday,datsam(:,4),'o--m');
plot (midday,datsam(:,5),'og');
plot (midday,datsam(:,6),'o--g');
plot (midday,datsam(:,7),'ob');
plot (midday,datsam(:,8),'o--b');
plot (midday,datigs(:,1),'r');
plot (midday,datigs(:,2),'--r');
plot (midday,datigs(:,3),'m');
plot (midday,datigs(:,4),'--m');
plot (midday,datigs(:,5),'g');
plot (midday,datigs(:,6),'--g');
plot (midday,datigs(:,7),'b');
plot (midday,datigs(:,8),'--b');
plot (midday,datmit(:,1),'-.r');	% 72.6% area 8.5TECU
%plot (midday,datmit(:,2),':r');	% 67.5% area bad hi 13.1TECU
plot (midday,datmit(:,3),'-.m');	% 88.4% area 7.8TECU
plot (midday,datmit(:,4),':m');		% 74.1% area 8.1TECU
%plot (midday,datmit(:,5),'-.g');	% 55.4% area 8.4TECU
%plot (midday,datmit(:,6),':g');	% 55.6% area bad lo 6.3TECU
plot (midday,datmit(:,7),'-.b');	% 96.5% area 8.2TECU
%plot (midday,datmit(:,9),':b');	% 42.8% area 7.9TECU lo
text(nday1+0.01*(nday2-nday1),ymit,'MIT','Color','k','FontSize',fsize,'FontWeight','bold');
text(nday1+0.08*(nday2-nday1),ymit,'73%','Color','r','FontSize',fsize,'FontWeight','bold');
text(nday1+0.16*(nday2-nday1),ymit,'88,74%','Color','m','FontSize',fsize,'FontWeight','bold');
text(nday1+0.29*(nday2-nday1),ymit,'96%','Color','b','FontSize',fsize,'FontWeight','bold');
plot(nday1+0.035*(nday2-nday1),ytie,'sc');
text(nday1+0.06*(nday2-nday1),ytie,'UAM4gps,__500km','Color','c','FontSize',fsize,'FontWeight','bold');
text(nday1+0.65*(nday2-nday1),ytie,'x TIE-WT','Color','k','FontSize',fsize,'FontWeight','bold');
text(nday1+0.80*(nday2-nday1),ytie,'x TIME-AMIE','Color','b','FontSize',fsize,'FontWeight','bold');
text(nday1+0.50*(nday2-nday1),ytie,'x TIE-Kp','Color','g','FontSize',fsize,'FontWeight','bold');
text(nday1+0.01*(nday2-nday1),yigs,'IGS','Color','k','FontSize',fsize,'FontWeight','bold');
text(nday1+0.01*(nday2-nday1),yt,'*global SAMI3 o glon','Color','k','FontSize',fsize,'FontWeight','bold');
text(nday1+0.39*(nday2-nday1),yt,'25,90E','Color','r','FontSize',fsize,'FontWeight','bold');
text(nday1+0.50*(nday2-nday1),yt,'140,175E','Color','m','FontSize',fsize,'FontWeight','bold');
text(nday1+0.66*(nday2-nday1),yt,'200,250E','Color','g','FontSize',fsize,'FontWeight','bold');
text(nday1+0.82*(nday2-nday1),yt,'285,345E','Color','b','FontSize',fsize,'FontWeight','bold');
ylabel('TECU','FontSize',fsize,'FontWeight','bold');
xlabel('2007 Daynumber','FontSize',fsize,'FontWeight','bold');
grid;
set (gca,'xlim',[nday1,nday2],'ylim',[tecn,tecx],'FontSize',fsize,'box','on');
hold off;


for nd=1:61
 rgtec4(nd) = gtec4(nd)/gtec4l(nd);
 rgtec4shp(nd) = gtec4shp(nd)/gtec4shpl(nd);
 rgtec4shm(nd) = gtec4shm(nd)/gtec4shml(nd);
 rgtec4eqb(nd) = gtec4eqb(nd)/gtec4eqbl(nd);
 rgtec4nhm(nd) = gtec4nhm(nd)/gtec4nhml(nd);
 rgtec4nhp(nd) = gtec4nhp(nd)/gtec4nhpl(nd);
 rgtec4m(nd) = gtec4m(nd)/gtec4l(nd);
 rgtec4shpm(nd) = gtec4shpm(nd)/gtec4shpl(nd);
 rgtec4shmm(nd) = gtec4shmm(nd)/gtec4shml(nd);
 rgtec4eqbm(nd) = gtec4eqbm(nd)/gtec4eqbl(nd);
 rgtec4nhmm(nd) = gtec4nhmm(nd)/gtec4nhml(nd);
 rgtec4nhpm(nd) = gtec4nhpm(nd)/gtec4nhpl(nd);
 dgtec4(nd) = gtec4(nd)-gtec4l(nd);
 dgtec4shp(nd) = gtec4shp(nd)-gtec4shpl(nd);
 dgtec4shm(nd) = gtec4shm(nd)-gtec4shml(nd);
 dgtec4eqb(nd) = gtec4eqb(nd)-gtec4eqbl(nd);
 dgtec4nhm(nd) = gtec4nhm(nd)-gtec4nhml(nd);
 dgtec4nhp(nd) = gtec4nhp(nd)-gtec4nhpl(nd);
 dgtec4m(nd) = gtec4m(nd)-gtec4l(nd);
 dgtec4shpm(nd) = gtec4shpm(nd)-gtec4shpl(nd);
 dgtec4shmm(nd) = gtec4shmm(nd)-gtec4shml(nd);
 dgtec4eqbm(nd) = gtec4eqbm(nd)-gtec4eqbl(nd);
 dgtec4nhmm(nd) = gtec4nhmm(nd)-gtec4nhml(nd);
 dgtec4nhpm(nd) = gtec4nhpm(nd)-gtec4nhpl(nd);
end

figure
hold on;
rn = 0.6;
rx = 3.1;
rb1 = 0.08*(rx-rn)+rn;
rb2 = 0.16*(rx-rn)+rn;
title('TEC ratios from UAM4 07325-08020','FontSize',15,'FontWeight','bold');
plot (midday,rgtec4,'x-k');
plot (midday,rgtec4shp,'x-g');
plot (midday,rgtec4shm,'x-r');
plot (midday,rgtec4eqb,'x-m');
plot (midday,rgtec4nhp,'x-b');
plot (midday,rgtec4nhm,'x-c');
plot (midday,rgtec4m,'o-k');
plot (midday,rgtec4shpm,'o-g');
plot (midday,rgtec4shmm,'o-r');
plot (midday,rgtec4eqbm,'o-m');
plot (midday,rgtec4nhpm,'o-b');
plot (midday,rgtec4nhmm,'o-c');
text (nday1+0.01*(nday2-nday1),rb2,'x=gps/500km; o=2000/500km for various mlat regions','Color','k','FontSize',fsize,'FontWeight','bold');
text (nday1+0.01*(nday2-nday1),rb1,'50-90S summer','Color','g','FontSize',fsize,'FontWeight','bold');
text (nday1+0.28*(nday2-nday1),rb1,'25-50S','Color','r','FontSize',fsize,'FontWeight','bold');
text (nday1+0.41*(nday2-nday1),rb1,'equator +/-25','Color','m','FontSize',fsize,'FontWeight','bold');
text (nday1+0.65*(nday2-nday1),rb1,'25-50N','Color','c','FontSize',fsize,'FontWeight','bold');
text (nday1+0.79*(nday2-nday1),rb1,'50-90N winter','Color','b','FontSize',fsize,'FontWeight','bold');
ylabel('Ratios','FontSize',fsize,'FontWeight','bold');
xlabel(' ','FontSize',fsize,'FontWeight','bold');
grid;
set (gca,'xlim',[nday1,nday2],'ylim',[rn rx],'FontSize',fsize,'box','on');
hold off;


figure
hold on;
rn = 0.;
rx = 7.;
rb1 = 0.08*(rx-rn)+rn;
rb2 = 0.16*(rx-rn)+rn;
rb2=4.2
rb1=1.0
title('TEC differences from UAM4 07325-08020','FontSize',15,'FontWeight','bold');
plot (midday,dgtec4,'x-k');
plot (midday,dgtec4shp,'x-g');
plot (midday,dgtec4shm,'x-r');
plot (midday,dgtec4eqb,'x-m');
plot (midday,dgtec4nhp,'x-b');
plot (midday,dgtec4nhm,'x-c');
plot (midday,dgtec4m,'o-k');
plot (midday,dgtec4shpm,'o-g');
plot (midday,dgtec4shmm,'o-r');
plot (midday,dgtec4eqbm,'o-m');
plot (midday,dgtec4nhpm,'o-b');
plot (midday,dgtec4nhmm,'o-c');
text (nday1+0.01*(nday2-nday1),rb2,'x=gps/500km; o=2000/500km for various mlat regions','Color','k','FontSize',fsize,'FontWeight','bold');
text (nday1+0.01*(nday2-nday1),rb1,'50-90S summer','Color','g','FontSize',fsize,'FontWeight','bold');
text (nday1+0.28*(nday2-nday1),rb1,'25-50S','Color','r','FontSize',fsize,'FontWeight','bold');
text (nday1+0.41*(nday2-nday1),rb1,'equator +/-25','Color','m','FontSize',fsize,'FontWeight','bold');
text (nday1+0.65*(nday2-nday1),rb1,'25-50N','Color','c','FontSize',fsize,'FontWeight','bold');
text (nday1+0.79*(nday2-nday1),rb1,'50-90N winter','Color','b','FontSize',fsize,'FontWeight','bold');
ylabel('TECU','FontSize',fsize,'FontWeight','bold');
xlabel(' ','FontSize',fsize,'FontWeight','bold');
grid;
set (gca,'xlim',[nday1,nday2],'ylim',[rn rx],'FontSize',fsize,'box','on');
hold off;
