%%
%% This program is written for plotting CP and HP fort.51 outputs from amie.F or advance.F
%%  for the study on different drivers for the TIEGCM presented at pre-SWW CCMC Apr 7, 2014
%% 4/14 BA Emery
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];

% Models or data/models (SuperDARN)

 year1 = 2006;
 nday1 = 347;
 nday2 = 351;

nt = 0;
nt5 = 0;
nt10 = 0;
for nd=nday1:nday2-1
  for nmin=1:1440
   nt = nt + 1;
   minutes(nt) = nt-1;
  end
  for m5=1:288
   nt5 = nt5 + 1;
   min5(nt5) = (m5-1)*5+(nd-nday1)*1440;
  end
  for m10=1:144
   nt10 = nt10 + 1;
   min10(nt10) = 5+(m10-1)*10+(nd-nday1)*1440;
  end
end
minutes(nt+1) = nt;
min5(nt5+1) = nt;


% # daynum hr min sec SH_HP NH_HP power av_power SH_CP NH_CP CP  
fheelis = sprintf('/hao/aim3/emery/tiegcm_kameleon/tiegcm-linux/heelis_2p5d_06345-50_hp_cp')
dheelis = dlmread(fheelis,'',2879,0);
fweimer = sprintf('/hao/aim3/emery/tiegcm_kameleon/tiegcm-linux/weimer05_dec06_hp_cp')
dweimer = dlmread(fweimer,'',2879,0);
fastra = sprintf('/hao/aim3/emery/tiegcm_kameleon/tiegcm-linux/astra_dec06_cp_hp')
dastra = dlmread(fastra,'',2,0);
fsdd = sprintf('/hao/aim3/emery/tiegcm_kameleon/tiegcm-linux/superdarn_dec06_cp')
dsdd = dlmread(fsdd,'',2,0);
fsdvt = sprintf('/hao/aim3/emery/tiegcm_kameleon/tiegcm-linux/sd_vt_dec06_cp')
dsdvt = dlmread(fsdvt,'',2,0);
nll=size(dsdvt)
for nm=1:nll(1)
 minsdvt(nm) = dsdvt(nm,4) + (dsdvt(nm,3)-nday1)*1440;
end

fswmf2 = sprintf('/hao/aim3/emery/tiegcm_kameleon/tiegcm-linux/swmf_2p5d_06347-9_cat_fort51')
infswmf2=fopen(fswmf2,'r');
lineskp = fgetl(infswmf2);
for nm=1:4320
% Skip 2 lines, read 1
lineskp = fgetl(infswmf2);
lineskp = fgetl(infswmf2);
linein2 = fgetl(infswmf2);
inswmf2=sscanf(linein2,'%f');
minswmf2(nm) = (inswmf2(2)-nday1)*1440 + inswmf2(3)*60 + inswmf2(4);
cpswmf2(nm) = inswmf2(12);
hpswmf2(nm) = inswmf2(8);
end



figure
zero = [0.,0.];
xtot = [nday1,nday2];
yx = 250.;
ytot = [0,yx];
xpk = [349+1.75/24.,349+1.75/24.];

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.50,0.87,0.40])	% 1/2
hold on;
title('High-Latitude Drivers for December 13-16, 2006','FontSize',15,'FontWeight','bold');
plot (nday1+minutes/1440,dheelis(:,12),'k');
plot (nday1+minutes/1440,dweimer(:,11),'b');
plot (nday1+min5/1440,dastra(:,8),'r');
plot (nday1+min10/1440,dsdd(:,8),'g');
plot (nday1+minsdvt/1440,dsdvt(:,7),'Color',dark_green);
plot (nday1+minswmf2/1440,cpswmf2,'m');
plot (xpk,ytot,'--m');
ylabel('CP (kV)','FontSize',fsize,'FontWeight','bold');
xlabel(' ','FontSize',fsize,'FontWeight','bold');
grid;
set (gca,'xlim',[nday1,nday2],'Xtick',([347:0.5:351]),'XTickLabel',{' ',' ',' ',' ',' ',' ',' ',' ',' '},'ylim',[0,yx],'FontSize',fsize,'box','on');
subplot('position',[0.10,0.10,0.87,0.40])	% 2/2
hold on;
plot (nday1+minutes/1440,dheelis(:,8),'k');
plot (nday1+minutes/1440,dweimer(:,7),'b');
plot (nday1+min5/1440,dastra(:,11),'r');
plot (nday1+minswmf2/1440,hpswmf2,'m');
plot (xpk,ytot,'--m');
text (nday1+0.03*(nday2-nday1),225,'Heelis Kp','Color','k','FontSize',fsize,'FontWeight','bold');
text (nday1+0.03*(nday2-nday1),200,'Weimer05 solar wind','Color','b','FontSize',fsize,'FontWeight','bold');
text (nday1+0.03*(nday2-nday1),175,'SuperDARN VT','Color',dark_green,'FontSize',fsize,'FontWeight','bold');
text (nday1+0.03*(nday2-nday1),150,'SuperDARN 50mlat','Color','g','FontSize',fsize,'FontWeight','bold');
text (nday1+0.03*(nday2-nday1),125,'ASTRA AMIE (SD+mags)','Color','r','FontSize',fsize,'FontWeight','bold');
text (nday1+0.03*(nday2-nday1),100,'SWMF','Color','m','FontSize',fsize,'FontWeight','bold');
ylabel('HP (GW)','FontSize',fsize,'FontWeight','bold');
xlabel(' ','FontSize',fsize,'FontWeight','bold');
grid;
set (gca,'xlim',[nday1,nday2],'ylim',[0,yx],'FontSize',fsize,'box','on');

hold off;
