%%
%% This program is written for plotting the Dec06 'global' TEC for 8 glons 
%% May 2012, BA Emery, HAO/NCAR

clear all; close all

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

% Model

 year1 = 2006;
 nday1 = 347;
 nday2 = 350;

%  Open 'global' TEC daily values from 5 IGS TEC

filecode = sprintf('/hao/aim3/emery/dec06/igs/cddis/global_tec_CODE_06347-9.dat');
fileesoc = sprintf('/hao/aim3/emery/dec06/igs/cddis/global_tec_ESOC_06347-9.dat');
filecJPL = sprintf('/hao/aim3/emery/dec06/igs/cddis/global_tec_cJPL_06347-9.dat');
filecUPC = sprintf('/hao/aim3/emery/dec06/igs/cddis/global_tec_cUPC_06347-9.dat');
filecIGS = sprintf('/hao/aim3/emery/dec06/igs/cddis/global_tec_cIGS_06347-9.dat');
dcode = dlmread(filecode,'',0,0);
desoc = dlmread(fileesoc,'',0,0);
dcjpl = dlmread(filecJPL,'',0,0);
dcupc = dlmread(filecUPC,'',0,0);
dcigs = dlmread(filecIGS,'',0,0);
%2006 347 12 13 27.50 12.71 92.50 12.42 142.50 12.48 177.50 13.05 202.50 13.40 252.50 14.25 287.50 14.20 347.50 13.58
for nd=1:3
 igsday(nd) = dcigs(nd,2)+0.5;
 for j=1:8
  datcode(nd,j) = dcode(nd,4+2*j);
  datesoc(nd,j) = desoc(nd,4+2*j);
  datcjpl(nd,j) = dcjpl(nd,4+2*j);
  datcupc(nd,j) = dcupc(nd,4+2*j);
  datcigs(nd,j) = dcigs(nd,4+2*j);
 end
end		% nd=1:3




% Open Emmert's 400km file
fil400e = sprintf('/hao/aim3/emery/dec06/tledens_dec_2006_400km.txt')
dat400e = dlmread(fil400e,'',12,0);
for nd=1:5
  day400(nd) = dat400e(nd,2)+0.5;
  drag400(nd) = 10.^dat400e(nd+1,7);
  msis400(nd) = 10.^dat400e(nd,6);
end
 
% Open hourly IMF file
filimf = sprintf('/hao/aim3/emery/dec06/hourly_imf_06346_350')
datimf = dlmread(filimf,'',0,0);

for nd=1:24*5
 uth = datimf(nd,3)/100.;
 midday(nd) = datimf(nd,2)+(uth+0.5)/24.;
 kp(nd) = datimf(nd,4); 
 vsw(nd) = datimf(nd,19);
 dsw(nd) = datimf(nd,18);
 bt(nd) = datimf(nd,11);
 bz(nd) = datimf(nd,15);
 by(nd) = datimf(nd,14);
 f107(nd) = datimf(nd,6);
end

% 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('Hourly 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],'xticklabel',{' '},'ylim',[0,9],'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],'xticklabel',{' '},'ylim',[-19 +16],'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],'xticklabel',{' '},'ylim',[0,21],'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('2006 Daynumber','FontSize',fsize,'FontWeight','bold');
grid;
set (gca,'xlim',[nday1,nday2],'ylim',[500,1000],'FontSize',fsize,'box','on');

tecx = 20;
tecn = 8;
yt = tecn + 0.90*(tecx-tecn);
yb = tecn + 0.10*(tecx-tecn);
ymit = 6.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 (igsday,datcigs(:,1),'r');
plot (igsday,datcigs(:,2),'--r');
plot (igsday,datcigs(:,3),'m');
plot (igsday,datcigs(:,4),'--m');
plot (igsday,datcigs(:,5),'g');
plot (igsday,datcigs(:,6),'--g');
plot (igsday,datcigs(:,7),'b');
plot (igsday,datcigs(:,8),'--b');
plot (igsday,datcode(:,3),'*m');
plot (igsday,datesoc(:,3),'om');
plot (igsday,datcjpl(:,3),'xm');
plot (igsday,datcupc(:,3),'sm');
text(nday1+0.01*(nday2-nday1),yb,'IGS 140E (*CODE oESOC xJPL sqUPC)','Color','k','FontSize',fsize,'FontWeight','bold');
%text(nday1+0.01*(nday2-nday1),yt,'IGS','Color','k','FontSize',fsize,'FontWeight','bold');
text(nday1+0.20*(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.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.1e-11;
y2 = 0.250e-11;
yt = y1 + 0.92*(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 Neutral Density at 400 km','FontSize',15,'FontWeight','bold');
plot (day400,drag400,'d-r');
plot (day400,msis400,'s-g');
text(nday1+0.01*(nday2-nday1),yt,'satellite drag','Color','r','FontSize',fsize,'FontWeight','bold');
text(nday1+0.25*(nday2-nday1),yt,'MSIS','Color','g','FontSize',fsize,'FontWeight','bold');
ylabel('kg/m3','FontSize',fsize,'FontWeight','bold');
xlabel('2006 Daynumber','FontSize',fsize,'FontWeight','bold');
grid;
set (gca,'xlim',[nday1,nday2],'ylim',[y1,y2],'FontSize',fsize,'box','on');
hold off;








figure
hold on;
tecx = 20.;
tecn = 10.;
yt = tecn + 0.95*(tecx-tecn);
yb = tecn + 0.05*(tecx-tecn);

plot (igsday,datcigs(:,1),'r');
plot (igsday,datcigs(:,2),'--r');
plot (igsday,datcigs(:,3),'m');
plot (igsday,datcigs(:,4),'--m');
plot (igsday,datcigs(:,5),'g');
plot (igsday,datcigs(:,6),'--g');
plot (igsday,datcigs(:,7),'b');
plot (igsday,datcigs(:,8),'--b');
plot (igsday,datcode(:,3),'*m');
plot (igsday,datesoc(:,3),'om');
plot (igsday,datcjpl(:,3),'xm');
plot (igsday,datcupc(:,3),'sm');

%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');
%text(nday1+0.05*(nday2-nday1),ytie,'x TIE-WT','Color','k','FontSize',fsize,'FontWeight','bold');
text(nday1+0.01*(nday2-nday1),yb,'IGS 140E (*CODE oESOC xJPL sqUPC)','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(' ','FontSize',fsize,'FontWeight','bold');
grid;
set (gca,'xlim',[nday1,nday2],'ylim',[tecn,tecx],'FontSize',fsize,'box','on');
hold off;



