close all

iyearN = 2005;
ndayN1=1;
ndayN2=90;

% Clear number of eq, pol, peaks to plot for 6 and 18 MLT
for n=1:4
 n4(n) = 0;
end
clear kev06eq;
clear kev06pol;
clear kev18eq;
clear kev18pol;
clear efx06eq;
clear efx06pol;
clear efx18eq;
clear efx18pol;

PlotTitle = sprintf('Smoothed E-Flux and Mean Energy Relation for %d %3.3d-%3.3d',iyearN,ndayN1,ndayN2);

for ndayN=ndayN1:ndayN2

if ndayN ~= 31 && ndayN ~= 59

k=1;
clear data;
fil1 = sprintf('/hao/aim3/dmsp/dmsp_rich/txt_dmspssj/%4.0d%3.3deq06.txt',iyearN,ndayN);
OUT1 = fopen(fil1,'r');
% Read in entire file skipping first 2 lines with dlmread
data = dlmread(fil1,'',2,0);
nrd = length(data);
for i=1:nrd
 n4(k) = n4(k) + 1;
 kev06eq(n4(k)) = data(i,6);
 efx06eq(n4(k)) = data(i,3);
end
k=2;
clear data;
fil2 = sprintf('/hao/aim3/dmsp/dmsp_rich/txt_dmspssj/%4.0d%3.3deq18.txt',iyearN,ndayN);
OUT2 = fopen(fil2,'r');
% Read in entire file skipping first 2 lines with dlmread
data = dlmread(fil2,'',2,0);
nrd = length(data);
for i=1:nrd
 n4(k) = n4(k) + 1;
 kev18eq(n4(k)) = data(i,6);
 efx18eq(n4(k)) = data(i,3);
end
k=3;
clear data;
fil3 = sprintf('/hao/aim3/dmsp/dmsp_rich/txt_dmspssj/%4.0d%3.3dpo06.txt',iyearN,ndayN);
OUT3 = fopen(fil3,'r');
% Read in entire file skipping first 2 lines with dlmread
data = dlmread(fil3,'',2,0);
nrd = length(data);
for i=1:nrd
 n4(k) = n4(k) + 1;
 kev06pol(n4(k)) = data(i,6);
 efx06pol(n4(k)) = data(i,3);
end
k=4;
clear data;
fil4 = sprintf('/hao/aim3/dmsp/dmsp_rich/txt_dmspssj/%4.0d%3.3dpo18.txt',iyearN,ndayN);
OUT4 = fopen(fil4,'r');
% Read in entire file skipping first 2 lines with dlmread
data = dlmread(fil4,'',2,0);
nrd = length(data);
for i=1:nrd
 n4(k) = n4(k) + 1;
 kev18pol(n4(k)) = data(i,6);
 efx18pol(n4(k)) = data(i,3);
end

end     % of if statement
end	% of for ndayN=ndayN1:ndayN2


cp18eflx =length(efx18pol)
cp18kev  =length(kev18pol)
cp06eflx =length(efx06pol)
ce18eflx =length(efx18eq)
ce06eflx =length(efx06eq)

% Make 4 plots of Mean smoothed (nsm) ekeV vs eflux (mW/m3)

fsize = 12;
fsiz = 11;

figure
emax = fix(max(efx06eq(efx06eq<100))+0.99);
kmax = fix(max(kev06eq)+0.99);
	plot(kev06eq,efx06eq, 'r+');
hold on;
	text (0.5,0.9*emax,'06MLT Equatorward','Color','r','fontsize',fsize, 'fontweight','bold');
	title(PlotTitle,'fontsize',fsize,'fontweight','bold');
	ylabel('E flux (mW/m^2)','fontsize',fsiz,'fontweight','bold');
	xlabel('<E> (keV)','fontsize',fsiz,'fontweight','bold');
	set (gca,'Xlim',[0 kmax],'Ylim',[0 emax]);
hold off;
emax4 = emax;
kmax4 = kmax;

figure
emax = fix(max(efx06pol(efx06pol<100))+0.99);
kmax = fix(max(kev06pol)+0.99);
	plot(kev06pol,efx06pol, 'b+');
hold on;
	text (0.5,0.9*emax,'06MLT Poleward','Color','b','fontsize',fsize, 'fontweight','bold');
	title(PlotTitle,'fontsize',fsize,'fontweight','bold');
	ylabel('E flux (mW/m^2)','fontsize',fsiz,'fontweight','bold');
	xlabel('<E> (keV)','fontsize',fsiz,'fontweight','bold');
	set (gca,'Xlim',[0 kmax],'Ylim',[0 emax]);
hold off;
emax4 = max(emax4,emax);
kmax4 = max(kmax4,kmax);


figure
emax = fix(max(efx18eq(efx18eq<100))+0.99);
kmax = fix(max(kev18eq)+0.99);
	plot(kev18eq,efx18eq, 'r+');
hold on;
	text (0.5,0.9*emax,'18MLT Equatorward','Color','r','fontsize',fsize, 'fontweight','bold');
	title(PlotTitle,'fontsize',fsize,'fontweight','bold');
	ylabel('E flux (mW/m^2)','fontsize',fsiz,'fontweight','bold');
	xlabel('<E> (keV)','fontsize',fsiz,'fontweight','bold');
	set (gca,'Xlim',[0 kmax],'Ylim',[0 emax]);
hold off;
emax4 = max(emax4,emax);
kmax4 = max(kmax4,kmax);


figure
emax = fix(max(efx18pol(efx18pol<100))+0.99);
kmax = fix(max(kev18pol)+0.99);
	plot(kev18pol,efx18pol, 'b+');
hold on;
	text (0.5,0.9*emax,'18MLT Poleward','Color','b','fontsize',fsize, 'fontweight','bold');
	title(PlotTitle,'fontsize',fsize,'fontweight','bold');
	ylabel('E flux (mW/m^2)','fontsize',fsiz,'fontweight','bold');
	xlabel('<E> (keV)','fontsize',fsiz,'fontweight','bold');
	set (gca,'Xlim',[0 kmax],'Ylim',[0 emax]);
hold off;
emax4 = max(emax4,emax);
kmax4 = max(kmax4,kmax);

% Plot all together
figure
	plot(kev06eq,efx06eq, 'r+');
hold on;
	plot(kev06pol,efx06pol, 'b+');
	plot(kev18eq,efx18eq, 'm+');
	plot(kev18pol,efx18pol, 'c+');
	text (0.70*kmax4,0.95*emax4,'Poleward MLT','Color','k','fontsize',fsize, 'fontweight','bold');
	text (0.75*kmax4,0.90*emax4,'06','Color','r','fontsize',fsize, 'fontweight','bold');
	text (0.85*kmax4,0.90*emax4,'18','Color','m','fontsize',fsize, 'fontweight','bold');
	text (0.70*kmax4,0.85*emax4,'Equatorward MLT','Color','k','fontsize',fsize, 'fontweight','bold');
	text (0.75*kmax4,0.80*emax4,'06','Color','b','fontsize',fsize, 'fontweight','bold');
	text (0.85*kmax4,0.80*emax4,'18','Color','c','fontsize',fsize, 'fontweight','bold');
	title(PlotTitle,'fontsize',fsize,'fontweight','bold');
	ylabel('E flux (mW/m^2)','fontsize',fsiz,'fontweight','bold');
	xlabel('<E> (keV)','fontsize',fsiz,'fontweight','bold');
	set (gca,'Xlim',[0 kmax4],'Ylim',[0 emax4]);
hold off;



