close all

fxmin = 0.25;

% Read SSJ4 .txt file for a day 
hemisphere = 'N';
hemisphere = 'S';
%filname = '/home/liamk/data/cenergy/cenergy_all_f15_20050514.txt';
%filname = '/home/emery/reu11/cenergy_all_f15_20050121.txt';
iyear = 2005;
iday = 21;
Bzav = 7;
Byav = 12;
nskp = 70801;
nrd = 72587-nskp;
poymax = 20.028;
poyx = [20.0 20.13];
emax = 17.;
%filname = '/home/emery/reu11/cenergy_all_f15_20041107.txt';
%poyname = '/hao/aim3/dmsp/poynting_flux/pfomni_04312.txt';
% 19UT 67763-69404 (fracday mlat MLT ? ? ? PoynF ? ? By Bz Vsw Dsw ? of 14)
iyear = 2004;
iday = 312;
Bzav = 35;
Byav = 25;
nskp = 67211;
nrd = 68965-nskp;
nskpp = 67763
nrdp = 69404-nskpp;
x1 = 18.8;
x2 = 19.4;
%poymax = 19.15;
%poyx = [19.1 19.3];
emax = 10.
poyn = -200.
filname = '/home/emery/reu11/cenergy_all_f15_20001108.txt';
poyname = '/hao/aim3/dmsp/poynting_flux/pfomni_00313.txt';
iyear = 2000;
iday = 313;
% 22104-23648 for SH 6.3UT minpoy=-0.03!!  (21549-23226ssj)  and nn=0 problems!
hemisphere = 'S';
nskp = 21549;
nrd = 14851-nskp;
nskpp = 22104
nrdp = 23648-nskpp;
x1 = 6.0;
x2 = 7.0;
emax = 5.;
poyn = -50.;
Bzav = 10;
Byav = 0;
%12914-14504 for NH 3.7UT (13156-14851ssj4) 
hemisphere = 'N';
nskp = 13156;
nrd = 14851-nskp;
nskpp = 12914
nrdp = 14504-nskpp;
x1 = 3.6;
x2 = 4.2;
emax = 1.;
poyn = -12.;
kmax = 9.;
kmax = 12.;
fxmin = 0.03;
fxmin = 0.04;


yearS = int2str(iyear);
daynumS = int2str(iday);
PlotTitle = ([yearS,' ',daynumS,' ',hemisphere,'H  DMSP-F15 SSJ4'] );

% Read Poynting day 1 sec file with OMNI data
rdall=0;
if rdall==1
 datpf=dlmread(poyname,'');
 jj = 0;
% This loop goes very slowly - Why?
 for n=nskpp+1:nskpp+nrdp
  jj = jj + 1;
  poyf(jj,1) = datpf(n,6);
  by(jj,1) = datpf(n,7);
  bz(jj,1) = datpf(n,8);
  vsw(jj,1) = datpf(n,9);
  dsw(jj,1) = datpf(n,10);
 uthp(jj,1) = (datpf(n,1)-iday)*24.;
 mlatp(jj,1) = datpf(n,2);
 mltp(jj,1) = datpf(n,3);
 end
jj
else
PFO=fopen(poyname,'r');
%nj = 1;
for n=1:nskpp
 linpf=fgetl(PFO);
end
%  Changing to array(jj,nj) makes the loop over 1590 go VERY slowly after ~950 - WHY???
for jj=1:nrdp
 linpf=fgetl(PFO);
if (iyear > 2000) 
 pfdat=sscanf(linpf,'%f %f %f %f %f %f %f %s %s %f %f %f %f %f');
 poyf(jj) = pfdat(7);
 by(jj) = pfdat(10);
 bz(jj) = pfdat(11);
 vsw(jj) = pfdat(12);
 dsw(jj) = pfdat(13);
else
 pfdat=sscanf(linpf,'%f');
 poyf(jj) = pfdat(6);
 by(jj) = pfdat(7);
 bz(jj) = pfdat(8);
 vsw(jj) = pfdat(9);
 dsw(jj) = pfdat(10);
end
 uthp(jj) = (pfdat(1)-iday)*24.;
 mlatp(jj) = pfdat(2);
 mltp(jj) = pfdat(3);
%jj
end
jj=jj-1
uthlen=length(uthp)
utsp = uthp(jj)*3600.
maxpoy = max(poyf);
minpoy = min(poyf);
minmlatp = min(mlatp);
maxmlatp = max(mlatp);
end
% Get 5 min plasma data (every 300 sec)
n5 = int8((jj+299)/300)
for n=1:n5
 m1s = (n-1)*300 + 1;
 ut5m(n) = uthp(m1s);
 mlat5m(n) = mlatp(m1s);
 mlt5m(n) = mltp(m1s);
 dsw5m(n) = dsw(m1s);
 vsw5m(n) = vsw(m1s);
 by5m(n) = by(m1s);
 bz5m(n) = bz(m1s);
end

% Read SSJ4 data
SSJ=fopen(filname,'r');

for n=1:nskp
 topssj = fgetl(SSJ);
end
nn = 0;
for n=1:nrd
 lindat = fgetl(SSJ);
 weidat = sscanf(lindat,'%f');
 if abs(weidat(8)) > 45
  nn = nn + 1;
  uts(nn) = weidat(2);
  mlat(nn) = weidat(8);
  alat(nn) = abs(weidat(8));
  mlt(nn) = weidat(10);
  eflx(nn) = weidat(11);
  ekev(nn) = weidat(13)*0.001;
  iflx(nn) = weidat(14);
  ikev(nn) = weidat(16)*0.001;
  if (eflx(nn)<fxmin)
   ekev(nn) = 0.;
  end
  if (iflx(nn)<fxmin)
   ikev(nn) = 0.;
  end
 end
end
uts1 = uts(1)
% 05134 nn=1675 but uts is filled to 21264 values for some reason.  Therefore limit plot to nn
nn
utslen=length(uts)
utsnn = uts(nn)
minuts = min(uts(1:nn))
maxuts = max(uts(1:nn))
maxeflx = max(eflx(1:nn))
maxiflx = max(iflx(1:nn))
maxekev = max(ekev(1:nn))
maxikev = max(ikev(1:nn))
minmlat = min(mlat(1:nn))
maxmlat = max(mlat(1:nn))

uthp1=uthp(1)
mlatp1=mlatp(1)
mltp1=mltp(1)
poyf1=poyf(1)

length(poyf)
% will label start of pass in hours=uts/3600.
%data1 = [uts/3600. alat mlt eflx];  
for n=1:jj
data2(n,1) = uthp(n);
data2(n,2) = mlatp(n);
data2(n,3) = mltp(n);
data2(n,4) = abs(poyf(n));
end
%data2 = [uthp(:,1) mlatp(:,1) mltp(:,1) abs(poyf(:,1))];	% works!
   [out,track,lines] = satelliteTrackPlot(data2,[0 10], 15, PlotTitle, 1,0,{});
for n=1:jj
 data2(n,4) = dsw(n);
end
%for n=1:n5
% data5(n,1) = ut5m(n);
% data5(n,2) = mlat5m(n);
% data5(n,3) = mlt5m(n);
% data5(n,4) = dsw5m(n);
%end
   [out,track,lines] = satelliteTrackPlot(data2,[0 20], 15, PlotTitle, 1,0,{});
%stop

% top plot
dark_green = [0 0.5 0];
fsize = 12;
fsiz = 11;
	H1=figure;
		subplot ('Position',[0.1,0.70,0.83,0.20])
	plot(uts(1:nn)/3600.,iflx(1:nn),'c');
	hold on;
	plot(uts(1:nn)/3600.,eflx(1:nn),'g');
	text (uts(1)/3600.,0.65*emax,'e-','Color','g','fontsize',fsize, 'fontweight','bold');
	text (uts(100)/3600.,0.65*emax,'p+','color','c','fontsize',fsize, 'fontweight','bold');
grid;
                        ylabel('E flux (mW/m^2)','fontsize',fsize,'fontweight','bold')
	                title (PlotTitle,'fontsize',fsize,'fontweight','bold');
	set (gca,'Xlim',[x1 x2],'Ylim',[0 emax],'XTickLabel',{' ';});
		subplot ('Position',[0.1,0.50,0.83,0.20])
	plot(uts(1:nn)/3600.,ikev(1:nn),'c');
hold on;
	plot(uts(1:nn)/3600.,ekev(1:nn),'g');
% MLT and Dsw plots to save space
%	plot(uts(1:nn)/3600.,mlt(1:nn),'b');
%       plot(uthp,dsw,'k');
grid;
         ylabel('<E> (keV)','fontsize',fsize,'fontweight','bold')
	set (gca,'Xlim',[x1 x2],'Ylim',[0 kmax],'XTickLabel',{' ';});
% Plot Poynting flux
		subplot ('Position',[0.1,0.40,0.83,0.10])
hold on;
	plot(uthp,poyf,'r');
grid;
         ylabel('Poy','fontsize',fsize,'fontweight','bold')
	set (gca,'Xlim',[x1 x2],'Ylim',[poyn 0],'XTickLabel',{' ';});
% Plot hourly By and Bz (nT)
		subplot ('Position',[0.1,0.30,0.83,0.10])
hold on;
	plot(uthp,bz,'r');
        plot(uthp,by,'color',dark_green);
	text (uts(1)/3600.,10.,'y','Color',dark_green,'fontsize',fsize, 'fontweight','bold');
	text (uts(1)/3600.,-10.,'z','Color','r','fontsize',fsize, 'fontweight','bold');
grid;
        ylabel('By,z','fontsize',fsiz,'fontweight','bold')
	set (gca,'Xlim',[x1 x2],'Ylim',[-15 15],'XTickLabel',{' ';});
% bottom plots MLAT and MLT
% Plot MLT and Dsw 
		subplot ('Position',[0.1,0.20,0.83,0.10])
	hold on;
	plot(uts(1:nn)/3600.,mlt(1:nn),'b');
%       plot(uthp,mltp,'k');
        plot(uthp,dsw,'k');
grid;
        ylabel('MLT,D','fontsize',fsiz,'fontweight','bold')
	set (gca,'Xlim',[x1 x2],'Ylim',[0 24],'XTickLabel',{' ';});
% Plot MLAT axis
		subplot ('Position',[0.1,0.10,0.83,0.10])
	hold on;
	plot(uts(1:nn)/3600.,mlat(1:nn),'b');
        plot(uthp,mlatp,'k');
% Plot position of max Poynting flux
%poyy = [-80.,-80.];
%	plot(poyx,poyy,'r');
%        plot(poymax,-80.,'xr','MarkerSize',10);
%	text(uts(1)/3600.+0.05,-50.,'Poynting flux <-20mW/m^2 (X=-180)','fontsize',fsize, 'fontweight','bold','color','r')
grid;
        ylabel('MLAT','fontsize',fsiz,'fontweight','bold')
        xlabel('UT (h)','fontsize',fsize,'fontweight','bold')
	set (gca,'Xlim',[x1 x2]);


hold off;


