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_20050123.txt';
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 = 300;
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';
filname = '/hao/aim3/dmsp/dmsp_rich/txt_dmspssj/2005/13_20050123.txt';
% 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;
poyn = -50.;
Bzav = 10;
Byav = 0;
%12914-14504 for NH 3.7UT (13156-14851ssj4) 
hemisphere = 'S';
nskp = 13156;
nrd = 14851-nskp;
nskpp = 12914;
nrdp = 14504-nskpp;
x1 = 3.6;
x2 = 4.2;
emax = 5;
poyn = -12.;
kmax = 9.;
kmax = 12.;
fxmin = 0.03;
fxmin = 0.04;
cos_85 = 0.087155743;
Sa = 128.30;

%Getting the date from the filename

inst   = filname(43:44);
iyear  = filname(46:49);  
iday   = filname(52:53);
imonth = filname(50:51);

PlotTitle = ([imonth, ' ', iyear,' ',iday,' ',hemisphere,'H  DMSP-F', inst,'SSJ4'] );

iyearN = str2num(iyear);
idayN  = str2num(iday);

% Read Poynting day 1 sec file with OMNI data

% 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);
  glat(nn) = weidat(6); 
  glon(nn) = weidat(7);
  mlat(nn) = weidat(8);
  mlon(nn) = weidat(9);
  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
eflx
uts1 = uts(1);
% 05134 nn=1675 but uts is filled to 21264 values for some reason.  Therefore limit plot to nn
utslen  = length(uts);
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));


% Particle Induced Conductivity Calculations

for n = 1:nn
	condp(n) = 40 * (ekev(n)) * (eflx(n))^(1/2) / (16 + (ekev(n))^2);

	condh(n) = 0.45 * (ekev(n))^0.85 * condp(n);

end 

% Calculating conductivities from UV/EUV
%   csza = solar zenith angle
%   slat = co-latitude of the sub solar point
%   slon = longitude of the sub solar point
%   Sa   = solar radio flux in W/m^2 Hz
%   condUVp = conductivity from solar UV, EUV, and cosmic rays
%   condp_Tot = total Pederson conductivity
	
 for n = 1:nn
	
	[slat(n) slong(n)] = subsol(iyearN, idayN, uts(n)/3600, uts(n)/60, uts(n));

	csza(n) = cossza(glat(n), glon(n), slat(n), slong(n));	

	if csza(n) > cos_85

  		condUVp(n) = 0.88 * (Sa * csza(n))^(1/2) + 0.1; 
		condUVh(n) = 1.5 * (Sa * csza(n))^(1/2) +0.2;

	else 
		f(n)  = 0.06 + exp(1.803 * tanh(3.833 * csza(n)) + (0.5 * csza(n)) - 2.332);  
		
		Bp(n) = (1 - (0.99524 * (sin(slat(n)))^2)) * (1 + (0.3 * cos(slat(n))^2)); 

		condUVp(n) = 12.5 * (Sa/180)^(1/2) * f(n) / Bp(n);
	end	
	
	condp_Tot(n) = (condp(n)^2 + condUVp(n)^2)^(1/2);		
end	


% Averaging eflx using tsmovavg function
% smooth(y) does a moving smoothing over default of 5 values (20 sec for 4 sec data); 
% smooth(y,span) where span MUST BE ODD - hence smooth(y,10) is converted to the default 5

aveflx = smooth(eflx, 19);

% Setting maximum value for eflx on plot
emax = max(aveflx)+1;

%Averaging condh data

avecondh = smooth(condh,5);

%Filtering bad data on Ekev plot
% Do not plot E-kev if E-flux less than 0.20

for  n = 1:nn 

	if aveflx(n) <= 0.20

		ekev(n) = 0;
	end

end 

% will label start of pass in hours=uts/3600.
% data1 = [uts/3600. alat mlt eflx];  
for n=1:nn
data2(n,1) = uts(n)/3600.;
data2(n,2) = mlat(n);
data2(n,3) = mlt(n);
data2(n,4) = eflx(n);
end
%data2 = [uthp(:,1) mlatp(:,1) mltp(:,1) abs(poyf(:,1))];	% works!
   [out,track,lines] = satelliteTrackPlot(data2,[0 10], 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])
	hold on;
	plot(uts(1:nn)/3600.,aveflx(1:nn),'g');
	
% Plotting peaks on top of eflx plot


eflxmin = .5; 


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])
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',{' ';});

    
% Plotting conductivities
subplot ('Position',[0.1,0.30,0.83,0.10])
hold on;
	
	plot(uts(1:nn)/3600., condp(1:nn),'r'); 
%	plot(uts(1:nn)/3600., avecondh(1:nn), 'y');
	plot(uts(1:nn)/3600., condp_Tot(1:nn), 'k');
	plot(uts(1:nn)/3600., condUVp(1:nn), 'm');
	text (uts(1)/3600.,10,'\Sigma_P','Color','r','fontsize',fsize, 'fontweight','bold');
	text (uts(100)/3600.,10,'\Sigma_PTot','color','k','fontsize',fsize, 'fontweight','bold');
	text (uts(300)/3600.,10,'\Sigma_PUV/EUV','color','m','fontsize',fsize, 'fontweight','bold');
        
grid;
        ylabel('\Sigma_P, \Sigma_PTot', 'fontsize',fsiz,'fontweight','bold')
	set (gca,'Xlim',[x1 x2],'Ylim',[0 20],'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');

grid;
        ylabel('MLT','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');

grid;
        ylabel('MLAT','fontsize',fsiz,'fontweight','bold')
        xlabel('UT (h)','fontsize',fsize,'fontweight','bold')
	set (gca,'Xlim',[x1 x2]);


hold off;
