close all

% Part 1 to read Omni data needs to be done only once per run
enddoy=90;                 % always reads from doy=1
noflines = enddoy*24*12;   % 12 = 60min/h / 5 min

time = zeros(noflines,5);  % year,doy,ut,min
imf  = zeros(noflines,4);  % bt,bx,by,bz

% specify file to read in & open
file_omni = sprintf('/hao/aim3/dmsp/omni/omni_5min2005.asc');
fid_omni  = fopen(file_omni,'r');

% read on line at a time
% checked values with omni website http://omniweb.gsfc.nasa.gov/cgi/nx1.cgi
for nline=1:noflines
   line_in = fgets(fid_omni);
   data    = sscanf(line_in,'%d %d %d %d %*f %*f %*f %*f %*f %*f %*f %*f %*f %f %f %*f %*f %f %f');
   %data
   time(nline,1) = data(1);  %  year
   time(nline,2) = data(2);  %  doy
   time(nline,3) = data(3);  %  ut
   time(nline,4) = data(4);  %  min
   time(nline,5) = data(3)*3600+data(4)*60;  %  ut in sec
   imf(nline,1)  = data(5);  %  bt
   imf(nline,2)  = data(6);  %  bx
   imf(nline,3)  = data(7);  %  by
   imf(nline,4)  = data(8);  %  bz
end   % end reading in data from omni file

% find time with fill value
for nline=1:noflines
   if(imf(nline,:)>1000)
     ifound = 0;
     i=0;
     while ifound<1 %find first real value after this fill value 
       if(imf(nline+i,:)<1000)
         ifound = 1;
	 inext=nline+i;
       end
       i=i+1;  % increase by 1
     end
     x1= time(nline-1,2)*86400+time(nline-1,5);  % 24*60*60 sec/day
     x2= time(inext,2)*86400+time(inext,5);  % 24*60*60 sec/day
     x3= time(nline,2)*86400+time(nline,5);
     x= [x1,x2]; 
     for nimf=1:4 
        y= [imf(nline-1,nimf),imf(inext,nimf)];   
        imf(nline,nimf) = interp1(x,y,x3);
     end	
   end
end   % end reading in data from omni file

% close omni file
fclose(fid_omni) 
% end reading in omni data / part 1

iyearN = 2005;
ndayN1=60;
ndayN2=90;

% Clear number of eq, pol, peaks to plot for 6 and 18 MLT
for n=1:6
 n6(n) = 0;
end
clear mlt06eq;
clear mlt06pk;
clear mlt06pol;
clear mlt18eq;
clear mlt18pk;
clear mlt18pol;
clear rf2oeq06;
clear rf2oeq18;
clear rf2opol06;
clear rf2opol18;
clear peak06;
clear peak18;
clear delh06eq;
clear delh06pol;
clear delh18eq;
clear delh18pol;

PlotTitle = sprintf('Peak and Fall Off Comparisons for %d %3.3d-%3.3d',iyearN,ndayN1,ndayN2);
PlotTitle1 = sprintf('E flux Peak Comparisons for %d %3.3d-%3.3d',iyearN,ndayN1,ndayN2);

for ndayN=ndayN1:ndayN2
%ndayN = ndayN1

% Open read file of calculated smoothed peaks, 1/e fall-offs, dmlat, etc (w revised flags and revised numh!)
%  NOTE:  After the flags are changed (often for no polar cap, polar cap arcs, or low peaks), then add 'r' to name of file so is not destroyed accidentally



%%filcalc = sprintf('calc_ssj4_%4.0d_%3.3dr_sh',iyearN,ndayN);  % Test - needed to clear
filcalc = sprintf('/hao/aim3/dmsp/dmsp_rich/txt_dmspssj/calc_ssj4_%4.0d_%3.3dr',iyearN,ndayN);
CALC = fopen(filcalc,'r');
%  Read first line
 linhead = fgetl(CALC);
 headdat = sscanf(linhead,'%d %d %d %d %f %d %d %f');
 jyearN = headdat(1)
 jdayN = headdat(2)
 numh(1) = headdat(3);
 numh(2) = headdat(4);
 fxmin = headdat(5);
 nsm = headdat(6);
 tossx = headdat(7);
 meanflux = headdat(8);


% Read numh(1) SH and numh(2) NH polar passes (revised downwards for passes skipped if ixmlat<55deg)

for ih=1:2
 for num=1:numh(ih)
% Read 6 flags (delh06eq peak06 delh06pol delh18pol peak18 peak18eq)
  linflg = fgetl(CALC);
  flgdat = sscanf(linflg,'%d');
  for i=1:6
   iflg(i) = flgdat(i);
   if (iflg(i)==0)
    n6(i) = n6(i) + 1;
   end
  end
% Read max mlat (must be >minmlat or 55deg - if not, then have a missing polar pass!)
% num,ih,ixmlat,uts(ixmlat),mlat(ixmlat),mlt(ixmlat),aveflx(ixmlat)
  linxmlat = fgetl(CALC);
  xmlatdat = sscanf(linxmlat,'%f');
  jnum=xmlatdat(1);		% won't be the same as num if a polar pass was skipped

% Read 06 mlt eq side
% delh(i),origsum0(i),origsum(i),nsmsum0(i),nsmsum(i),gaussum0(i),gaussum(i),istart(i),uts(istart(i)),mlat(istart(i)),mlt(istart(i)),aveflx(istart(i)),fall(i),uts(fall(i)),mlat(fall(i)),mlt(fall(i)),aveflx(fall(i)),ipeak1,uts(ipeak1),mlat(ipeak1),mlt(ipeak1),aveflx(ipeak1),peak1
  lineq = fgetl(CALC);
  eqdat = sscanf(lineq,'%f');
  utsec_peak06 = eqdat(19);   % since peak and fall of UT are close together and we use 5 min omni data
                              % we use the time of the peak for the imf conitions
                              % in case only one has good quailty we get utsec from both
  if (iflg(1)==0)
   delh06eq(n6(1)) = eqdat(1);
   rf2oeq06(n6(1)) = eqdat(7)/eqdat(3);
   rf2o0eq06(n6(1)) = eqdat(6)/eqdat(2);
   mlt06eq(n6(1)) = eqdat(16);
  end
  if (iflg(2)==0)
   peak06(n6(2)) = eqdat(23);
   mlt06pk(n6(2)) = eqdat(21);
  end
% Read 06 MLT pol side
% delh(i),origsum0(i),origsum(i),nsmsum0(i),nsmsum(i),gaussum0(i),gaussum(i),istart(i),uts(istart(i)),mlat(istart(i)),mlt(istart(i)),aveflx(istart(i)),fall(i),uts(fall(i)),mlat(fall(i)),mlt(fall(i)),aveflx(fall(i))
  linpol = fgetl(CALC);
  poldat = sscanf(linpol,'%f');
  if (iflg(3)==0)
   delh06pol(n6(3)) = poldat(1);
   rf2opol06(n6(3)) = poldat(7)/poldat(3);
   rf2o0pol06(n6(3)) = poldat(6)/poldat(2);
   mlt06pol(n6(3)) = poldat(16);
   utsec_peak06 = eqdat(19); 
  end
% Read 18 MLT pol side
  linpol = fgetl(CALC);
  poldat = sscanf(linpol,'%f');
  if (iflg(4)==0)
   delh18pol(n6(4)) = poldat(1);
   rf2opol18(n6(4)) = poldat(7)/poldat(3);
   rf2o0pol18(n6(4)) = poldat(6)/poldat(2);
   mlt18pol(n6(4)) = poldat(16);
  end
% Read 18 mlt eq side
  lineq = fgetl(CALC);
  eqdat = sscanf(lineq,'%f');
  utsec_peak18 = eqdat(19);   % since peak and fall of UT are close together and we use 5 min omni data
                              % we use the time of the peak for the imf conitions
                              % in case only one has good quailty we get utsec from both
  if (iflg(6)==0)
   delh18eq(n6(6)) = eqdat(1);
   rf2oeq18(n6(6)) = eqdat(7)/eqdat(3);
   rf2o0eq18(n6(6)) = eqdat(6)/eqdat(2);
   mlt18eq(n6(6)) = eqdat(16);
  end
  if (iflg(5)==0)
   peak18(n6(5)) = eqdat(23);
   mlt18pk(n6(5)) = eqdat(21);
  end
  % calculates IMF condition for DMSP data based on Omni data
  % make  idoy_dmsp, utsec_dmsp variable depending on your data

  %idoy_dmsp = 9             % day number of DMSP data
  utsec_dmsp = 45541         % UT in sec of DMSP data
  idoy_dmsp = jdayN
  utsec_dmsp = utsec_peak06

  % here the time of the dmsp data is used
  % needs to propagate the solar wind measurement from 1 AU to Earth
  % calculate index to use for OMNI data based on DMSP data OMNI is 5 min data
  % 288: # of index for one day based on 5 min data
  % 300: to convert utsec in 5min slices-> 1min/60s/5min=1/300s
  
  i5m = (idoy_dmsp-1)*288 + cast((utsec_dmsp/300.),'int32') + 1;
  if (i5m < 5) 
      
  %  Fix this later for first 25 min of day
    bx15 = imf(i5m,2);
    by15 = imf(i5m,3);
    bz15 = imf(i5m,4);
    bt15 = imf(i5m,1);
  else
    bx15 = (imf(i5m-2,2) + imf(i5m-3,2) + imf(i5m-4,2))/3.;
    by15 = (imf(i5m-2,3) + imf(i5m-3,3) + imf(i5m-4,3))/3.;
    bz15 = (imf(i5m-2,4) + imf(i5m-3,4) + imf(i5m-4,4))/3.;
    bt15 = (imf(i5m-2,1) + imf(i5m-3,1) + imf(i5m-4,1))/3.;
  end; % get IMF conditions

 end	% of num=1:numh(ih)
end	% of ih=1:2

end	% of for ndayN=ndayN1:ndayN2


% Bins the peaks,delh,and ratio according to 1 hour MLT

ende06 = length(mlt06eq);
endp06 = length(mlt06pol);
ende18 = length(mlt18eq);
endp18 = length(mlt18pol);
endpk06 = length(peak06);
endpk18 = length(peak18);
je06 = 0;
je18 = 0;
jp06 = 0;
jp18 = 0;
jpk06 = 0;
jpk18 = 0;


for j=1:12
	for i = 1:ende06
		if mlt06eq(i) >= j && mlt06eq(i) < (j+1)
			je06 = je06+1;
			bmlt06eq(je06,j) = mlt06eq(i);
			brfeq06(je06,j) = rf2oeq06(i);
			bdelh06eq(je06,j) = delh06eq(i);
			
		end
	end
	je06 = 0;
end

for j=1:12	
	for i = 1:endp06
		if mlt06pol(i) >= j && mlt06pol(i) < (j+1)
			jp06 = jp06+1;
			bmlt06pol(jp06,j) = mlt06pol(i);
			brfpol06(jp06,j) = rf2opol06(i);
			bdelh06pol(jp06,j) = delh06pol(i);
		end
	end
	jp06=0;
end

j2=1;

for j=12:23;
	for i = 1:ende18
		if mlt18eq(i) >= j && mlt18eq(i) < (j+1)
			je18 = je18+1;
			bmlt18eq(je18,j2) = mlt18eq(i);
			brfeq18(je18,j2) = rf2oeq18(i);
			bdelh18eq(je18,j2) = delh18eq(i);
			
		end
	end
	je18 = 0;
	j2 = j2+1;
end

j2=1;
for j=12:23	
	for i = 1:endp18
		if mlt18pol(i) >= j && mlt18pol(i) < (j+1)
			jp18 = jp18+1;
			bmlt18pol(jp18,j2) = mlt18pol(i);
			brfpol18(jp18,j2) = rf2opol18(i);
			bdelh18pol(jp18,j2) = delh18pol(i);
		end
	end
	jp18=0;
	j2 = j2+1;
end

j2=1;
for j=1:12
	for i = 1:endpk06
		if mlt06pk(i) >= j && mlt06pk(i) < (j+1)
			jpk06 = jpk06+1;
			bpeak06(jpk06,j2) = peak06(i);
		end
	end
	jpk06=0;
	j2=j2+1;
end

j2=1;
for j=12:24
	for i = 1:endpk18
		if mlt18pk(i) >= j && mlt18pk(i) < (j+1)
			jpk18 = jpk18+1;
			bpeak18(jpk18,j2) = peak18(i);
		end
	end
	jpk18=0;
	j2=j2+1;
end

% Setting index lengths for plots

if length(rf2opol06) > length(rf2oeq06)
	nend = length(rf2oeq06);
else
	nend = length(rf2opol06);
end

if length(rf2opol18) > length(rf2oeq18)
	nend2 = length(rf2oeq18);
else
	nend2 = length(rf2opol18);
end

% Creates new array of peaks that hold certain conditions:
% The ratio of the fit over data must be less than 1.25 and greater than 0.75 for both polward and equatorward halfwidths 
% These new peak values will be plotted along with the old 

nu=0;
for i = 1:nend
	if rf2oeq06(i) <= 1.25 && rf2opol06(i) <= 1.25 && rf2oeq06(i) >= 0.75 && rf2opol06(i) >= 0.75
		nu=nu+1;	
		nmlt06pk(nu) = mlt06pk(i);
		npeak06(nu) = peak06(i);
		nrfeq06(nu) = rf2oeq06(i);
		nrfpol06(nu) = rf2opol06(i);
	end
end

nu2=0;

for i = 1:nend2
	if rf2oeq18(i) <= 1.25 && rf2opol18(i) <= 1.25 && rf2oeq18(i) >= 0.75 && rf2opol18(i) >= 0.75
		nu2 = nu2+1;
		nmlt18pk(nu2) = mlt18pk(i);
		npeak18(nu2) = peak18(i);
		nrfeq18(nu2) = rf2oeq18(i);
		nrfpol18(nu2) = rf2opol18(i);
	end
end

% This creates the BINNED arrays for the new peak values (for 06 and 18MLT)

nendpk06 = length(npeak06);
nendpk18 = length(npeak18);

j2=1;
jpk06=0;
for j=1:12
	for i = 1:nendpk06
		if nmlt06pk(i) >= j && nmlt06pk(i) < (j+1)
			jpk06 = jpk06+1;
			nbpeak06(jpk06,j2) = npeak06(i);
		end
	end
	jpk06=0;
	j2=j2+1;
end

j2=1;
jpk18=0;
for j=12:24
	for i = 1:nendpk18
		if nmlt18pk(i) >= j && nmlt18pk(i) < (j+1)
			jpk18 = jpk18+1;
			nbpeak18(jpk18,j2) = npeak18(i);
		end
	end
	jpk18=0;
	j2=j2+1;
end


% From binned data, take the mean, median, and standard deviation
%Examples:
%mlt18 = sum(bmlt18eq,1)./sum(bmlt18eq~=0,1);
%mlt06 = sum(bmlt06eq,1)./sum(bmlt06eq~=0,1)

%Calculates the mean for the hourly MLT

mpk06 = sum(bpeak06,1)./sum(bpeak06~=0,1);
mpk18 = sum(bpeak18,1)./sum(bpeak18~=0,1);
nmpk06 = sum(nbpeak06,1)./sum(nbpeak06~=0,1);
nmpk18 = sum(nbpeak18,1)./sum(nbpeak18~=0,1);
mrfeq06 = sum(brfeq06,1)./sum(brfeq06~=0,1);
mrfeq18 = sum(brfeq18,1)./sum(brfeq18~=0,1);
mrfpol18 = sum(brfpol18,1)./sum(brfpol18~=0,1);
mrfpol06 = sum(brfpol06,1)./sum(brfpol06~=0,1);
mdh06eq = sum(bdelh06eq,1)./sum(bdelh06eq~=0,1);
mdh18eq = sum(bdelh18eq,1)./sum(bdelh18eq~=0,1);
mdh06pol = sum(bdelh06pol,1)./sum(bdelh06pol~=0,1);
mdh18pol = sum(bdelh18pol,1)./sum(bdelh18pol~=0,1);
mrfeq06 = sum(brfeq06,1)./sum(brfeq06~=0,1);
mrfeq18 = sum(brfeq18,1)./sum(brfeq18~=0,1);
mrfpol06 = sum(brfpol06,1)./sum(brfpol06~=0,1);
mrfpol18 = sum(brfpol18,1)./sum(brfpol18~=0,1);

% Standard Deviation
% For peak raw data: stpk06 & stpk18 (uses binned values bpeak06 & bpeak18)
% For peak modified data: mstpk06 & mstpk18 (uses binned values nbpeak06 & nbpeak18)

nu=0;
for j=1:11
	for i=1:length(bpeak06)
		if bpeak06(i,j)~=0
			nu=nu+1;
			pkst06(nu) = bpeak06(i,j);
		end
	end
	stpk06(j) = std(pkst06);
end

nu=0;
for j=1:11
	for i=1:length(bpeak18)
		if bpeak18(i,j)~=0
			nu=nu+1;
			pkst18(nu) = bpeak18(i,j);
		end
	end
	stpk18(j) = std(pkst18);
end


nu=0;
for j=1:11
	for i=1:length(nbpeak06)
		if nbpeak06(i,j)~=0
			nu=nu+1;
			npkst06(nu) = nbpeak06(i,j);
		end
	end
	nstpk06(j) = std(npkst06);
end

nu=0;
for j=1:11
	for i=1:length(nbpeak18)
		if nbpeak18(i,j)~=0
			nu=nu+1;
			npkst18(nu) = nbpeak18(i,j);
		end
	end
	if exist('npkst18')==1
		nstpk18(j) = std(npkst18);
	else
		nstpk18(j) = 0;
	end
end



%MLT array
mltarr1 = 1:11;
mltarr2 = 13:23;




% Third Plot: plots the MLAT delta h for each 1/e fall off

fsize = 12;
fsiz = 11;

x3 = 0;
x4 = 24;
mmax = 24;
emax06 = max(peak06);
emax18 = max(peak18);
emax = fix(max(emax06,emax18)+0.99)
delmax06eq = max(delh06eq);
delmax06pol = max(delh06pol);
delmax06 = max(delmax06eq,delmax06pol)
delmax18eq = max(delh18eq);
delmax18pol = max(delh18pol);
delmax18 = max(delmax18eq,delmax18pol)
delmax = fix(max(delmax06,delmax18)+0.99)
delmax = min(15.,delmax)
rmax = 4.

figure(1);

% Plot of ratio of fit to original energy flux
	subplot ('Position',[0.1,0.60,0.83,0.30])
hold on;
	plot(mlt06eq,rf2oeq06,'r*')
	plot(mlt18eq,rf2oeq18,'r*')
	plot(mlt06pol,rf2opol06,'c*')
	plot(mlt18pol,rf2opol18,'c*')
	text (10.5,0.9*rmax,'Poleward fall off','Color','c','fontsize',fsize, 'fontweight','bold');
	text (10.5,0.75*rmax,'Equatorward fall off','Color','r','fontsize',fsize, 'fontweight','bold');
grid;
hold on;
	plot((0:mmax), ones(1,25),'color','k');
	title (PlotTitle,'fontsize',fsize,'fontweight','bold');
	ylabel('Ratio Fit/Data','fontsize',fsiz,'fontweight','bold');
	set (gca,'Xlim',[x3 x4], 'Ylim',[0 rmax]);	
	yL = get(gca,'YLim');
	line([6 6],yL,'color','k','LineStyle','--');
	line([18 18],yL,'color','k','LineStyle','--');
	text (6.5,0.9*rmax,'6 MLT','fontsize',fsize, 'fontweight','bold');
	text (18.5,0.9*rmax,'18 MLT','fontsize',fsize, 'fontweight','bold');

% Plot of peaks
	subplot ('Position',[0.1,0.40,0.83,0.20])
grid;
hold on;
	set (gca,'Xlim',[x3 x4], 'Ylim',[0 30]);
	plot(mlt06pk, peak06,'md','MarkerFaceColor','m');
	plot(mlt18pk, peak18,'md','MarkerFaceColor','m');
hold on;

	plot(nmlt06pk, npeak06,'kd','MarkerFaceColor','k');
	plot(nmlt18pk, npeak18,'kd','MarkerFaceColor','k');
	yL = get(gca,'YLim');
	line([6 6],yL,'color','k','LineStyle','--');
	line([18 18],yL,'color','k','LineStyle','--');	
	text (2,.8*emax,'Peaks','Color','m','fontsize',fsize, 'fontweight','bold');
	ylabel('E Flux (mW/m^2)','fontsize',fsiz,'fontweight','bold');

% Plot of delta(h) in mlat
	subplot ('Position',[0.1,0.10,0.83,0.30])
hold on;
	set (gca,'Xlim',[x3 x4], 'Ylim',[0,delmax]);
	plot(mlt06eq,delh06eq,'r*')
hold on;
	plot(mlt18eq,delh18eq,'r*')
	plot(mlt06pol,delh06pol,'c*')
	plot(mlt18pol,delh18pol,'c*')	
	yL = get(gca,'YLim');
	line([6 6],yL,'color','k','LineStyle','--');
	line([18 18],yL,'color','k','LineStyle','--');
	text (10.5,0.90*delmax,'Poleward \Deltah','Color','c','fontsize',fsize, 'fontweight','bold');
	text (10.5,0.75*delmax,'Equatorward \Deltah','Color','r','fontsize',fsize, 'fontweight','bold');
	ylabel('MLAT \Deltah','fontsize',fsiz,'fontweight','bold');
	xlabel('MLT','fontsize',fsiz,'fontweight','bold');
		
grid;
hold off;

figure(2)

	subplot ('Position',[0.1,0.60,0.83,0.30])
hold on;
	plot(mltarr1, mrfeq06(1:11), 'r*');
	plot(mltarr2, mrfeq18(1:11), 'r*');
	plot(mltarr1, mrfpol06(1:11), 'c*');
	plot(mltarr2, mrfpol18(1:11), 'c*');
	plot((0:mmax), ones(1,25),'color','k');
	set (gca,'Xlim',[x3 x4], 'Ylim',[0 rmax]);	
	yL = get(gca,'YLim');
	line([6 6],yL,'color','k','LineStyle','--');
	line([18 18],yL,'color','k','LineStyle','--');
	ylabel('Ratio Fit/Data','fontsize',fsiz,'fontweight','bold');

hold on;
grid;

	subplot ('Position',[0.1,0.40,0.83,0.20])
grid;
hold on;
	set (gca,'Xlim',[x3 x4],'Ylim',[-5 20]);
	plot(mltarr1,mpk06(1:11),'md','MarkerFaceColor','m');
	plot(mltarr2,mpk18(1:11),'md','MarkerFaceColor','m');
	errorbar(mltarr1,mpk06(1:11),stpk06(1:11));
	errorbar(mltarr2,mpk18(1:11),stpk18(1:11));
	yL = get(gca,'YLim');
	line([6 6],yL,'color','k','LineStyle','--');
	line([18 18],yL,'color','k','LineStyle','--');	
	text (2,.8*emax,'Peaks','Color','m','fontsize',fsize, 'fontweight','bold');
	ylabel('E Flux (mW/m^2)','fontsize',fsiz,'fontweight','bold');
	xlabel('MLT','fontsize',fsiz,'fontweight','bold');

	subplot ('Position',[0.1,0.10,0.83,0.30])
hold on;

	plot(1:11, mdh06eq, 'r*');
	plot(mltarr2, mdh18eq(1:11), 'r*');
	plot(mltarr1, mdh06pol(1:11), 'c*');
	plot(mltarr2, mdh18pol(1:11), 'c*');

for i = 1:11
	if mrfeq06(i) <= 1.30 && mrfpol06(i) <= 1.30
		plot(mltarr1(i),mdh06eq(i), 'k*');
		plot(mltarr1(i),mdh06pol(i), 'k*');
	end
end

for i = 1:11
	if mrfeq18(i) <= 1.30 && mrfpol18(i) <= 1.30
		plot(mltarr2(i),mdh18eq(i), 'k*');
		plot(mltarr2(i),mdh18pol(i), 'k*');
	end
end

	set (gca,'Xlim',[x3 x4], 'Ylim',[0,delmax]);
	yL = get(gca,'YLim');
	line([6 6],yL,'color','k','LineStyle','--');
	line([18 18],yL,'color','k','LineStyle','--');
	ylabel('MLAT \Deltah','fontsize',fsiz,'fontweight','bold');
grid;
hold off;

figure(3)

subplot ('Position',[0.1,0.60,0.83,0.30])
hold on;
	plot(nmlt06pk,nrfeq06,'r*')
	plot(nmlt18pk,nrfeq18,'r*')
	plot(nmlt06pk,nrfpol06,'c*')
	plot(nmlt18pk,nrfpol18,'c*')
hold on;
	plot((0:mmax), ones(1,25),'color','k');
	set (gca,'Xlim',[x3 x4], 'Ylim',[0,4]);
	yL = get(gca,'YLim');
	line([6 6],yL,'color','k','LineStyle','--');
	line([18 18],yL,'color','k','LineStyle','--');
	title (PlotTitle1,'fontsize',fsize,'fontweight','bold');
	ylabel('Ratio Fit/Data','fontsize',fsiz,'fontweight','bold');
	text (6.5,3.5,'6 MLT','fontsize',fsize, 'fontweight','bold');
	text (18.5,3.5,'18 MLT','fontsize',fsize, 'fontweight','bold');	
grid;

subplot ('Position',[0.1,0.40,0.83,0.20])
	plot(nmlt06pk, npeak06,'kd','MarkerFaceColor','k');
hold on;
	plot(nmlt18pk, npeak18,'kd','MarkerFaceColor','k');
	set (gca,'Xlim',[x3 x4], 'Ylim',[0,35]);
	yL = get(gca,'YLim');
	line([6 6],yL,'color','k','LineStyle','--');
	line([18 18],yL,'color','k','LineStyle','--');	
	text (2,.8*emax,'Peaks','Color','m','fontsize',fsize, 'fontweight','bold');
	ylabel('E Flux (mW/m^2)','fontsize',fsiz,'fontweight','bold');
grid;

subplot ('Position',[0.1,0.10,0.83,0.30])
hold on;
	plot(mltarr1,nmpk06(1:11),'kd','MarkerFaceColor','k');
	plot(mltarr2,nmpk18(1:11),'kd','MarkerFaceColor','k');
	errorbar(mltarr1,nmpk06(1:11),nstpk06(1:11));
	errorbar(mltarr2,nmpk18(1:11),nstpk18(1:11));
	set (gca,'Xlim',[x3 x4], 'Ylim',[-5,20]);
	yL = get(gca,'YLim');
	line([6 6],yL,'color','k','LineStyle','--');
	line([18 18],yL,'color','k','LineStyle','--');	
	ylabel('E Flux Mean','fontsize',fsiz,'fontweight','bold');
	xlabel('MLT','fontsize',fsiz,'fontweight','bold');
grid;
hold off;

figure(4)

subplot ('Position',[0.1,0.60,0.83,0.30])
hold on;
	plot(mlt06eq,rf2oeq06,'r*')
	plot(mlt18eq,rf2oeq18,'r*')
	plot(mlt06pol,rf2opol06,'c*')
	plot(mlt18pol,rf2opol18,'c*')
hold on;
	plot((0:mmax), ones(1,25),'color','k');
	text (10.5,0.9*rmax,'Poleward fall off','Color','c','fontsize',fsize, 'fontweight','bold');
	text (10.5,0.75*rmax,'Equatorward fall off','Color','r','fontsize',fsize, 'fontweight','bold');
	text (6.5,3.5,'6 MLT','fontsize',fsize, 'fontweight','bold');
	text (18.5,3.5,'18 MLT','fontsize',fsize, 'fontweight','bold');	
	yL = get(gca,'YLim');
	line([6 6],yL,'color','k','LineStyle','--');
	line([18 18],yL,'color','k','LineStyle','--');
	title (PlotTitle1,'fontsize',fsize,'fontweight','bold');
	ylabel('Ratio Fit/Data','fontsize',fsiz,'fontweight','bold');
	set (gca,'Xlim',[x3 x4], 'Ylim',[0 rmax]);	
grid;
	
subplot ('Position',[0.1,0.40,0.83,0.20])
	plot(mlt06pk, peak06,'md','MarkerFaceColor','m');
hold on;
	plot(mlt18pk, peak18,'md','MarkerFaceColor','m');
	yL = get(gca,'YLim');
	line([6 6],yL,'color','k','LineStyle','--');
	line([18 18],yL,'color','k','LineStyle','--');
	set (gca,'Xlim',[x3 x4], 'Ylim',[0 35]);
	ylabel('E Flux (mW/m^2)','fontsize',fsiz,'fontweight','bold');		
grid;

subplot ('Position',[0.1,0.10,0.83,0.30])

	set (gca,'Xlim',[x3 x4],'Ylim',[-5 20]);
	plot(mltarr1,mpk06(1:11),'md','MarkerFaceColor','m');
hold on;
	plot(mltarr2,mpk18(1:11),'md','MarkerFaceColor','m');
	errorbar(mltarr1,mpk06(1:11),stpk06(1:11));
	errorbar(mltarr2,mpk18(1:11),stpk18(1:11));
	set (gca,'Xlim',[x3 x4], 'Ylim',[-5,20]);
	yL = get(gca,'YLim');
	line([6 6],yL,'color','k','LineStyle','--');
	line([18 18],yL,'color','k','LineStyle','--');
	ylabel('E Flux Mean','fontsize',fsiz,'fontweight','bold');
	xlabel('MLT','fontsize',fsiz,'fontweight','bold');
grid;
hold off;

