close all

% fxmin is the flux minimum below which we ignore the electron flux (fxmini for ion flux, and fxminkev is the electron flux where ekev is set to 0)
fxmin = 0.04
fxmini = fxmin - 0.01
fxminkev = fxmin
% eflxmin = minimum peak eflux where we plot/print the Gaussian fits (always calculated)
eflxmin = 0.50;
% Ensure that 1/e > fxmin!!
eflxmin = 2.718*fxmin;

% nsm = number of ODD seconds over which we smooth  (IDM smooths over 5*4s=20s)
nsm = 25
% ifsm=-1,0,1 if (-1) skip plotting of SH and NH orbits on polar plot; (0) plot unsmoothed eflx, (1) plot smoothed eflx 
ifsm=1

% tossx = multiplier for mean(eflx - including 0's) above which we toss spikes
tossx = 100.

% minmlat = minimum abs(mlat) above which peaks are found (after tossing tossx*mean)
minmlat = 55.

cos_85 = 0.087155743;
Sa = 128.30;

mmax = 25; 

cbarmax = 10.;

% kmax = max value of mean keV plotted
kmax = 12.;

% Read SSJ4 .txt file for a day 
filname = '/hao/aim3/dmsp/dmsp_rich/txt_dmspssj/2005/13_20050110.txt'


%Getting the date from the filename

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

iyearN = str2num(iyear)
imonN = str2num(imonth)
idaymoN = str2num(iday)
% Need to calculate daynumber of year for subsol
date3(1) = iyearN;
date3(2) = imonN;
date3(3) = idaymoN;
% datenum converts to 'serial' daynumber (so subtract 732312.5 for 2005!)
ndayN = fix(datenum(date3)-732312)

% Read Poynting day 1 sec file with OMNI data

% Read SSJ4 data 1 sec (86400 sec per day)
SSJ=fopen(filname,'r');
% Read in entire file skipping first 2 lines with dlmread
data = dlmread(filname,'',2,0);
nrd = length(data(:,1))

% DO NOT USE data(nrd,18) because large arrays are not fast - just read after skipping 2 header lines
for n=1:2
 topssj = fgetl(SSJ);
end
nn = 0;
numh(1) = 0;
numh(2) = 0;
for n=1:nrd
 lindat = fgetl(SSJ);
 weidat = sscanf(lindat,'%f');
 if abs(weidat(8)) <= 45
   nn = 0;
 end
 if abs(weidat(8)) > 45
  if (nn == 0 & weidat(8)<0) 
   ih = 1;
   numh(ih) = numh(ih) + 1
  end
  if (nn == 0 & weidat(8)>0) 
   ih = 2;
   numh(ih) = numh(ih) + 1
  end
  nn = nn + 1;
  num = numh(ih);
  nnsav(num,ih) = nn;
  uts3(nn,num,ih)  = weidat(2);
  glat(nn,num,ih) = weidat(6); 
  glon(nn,num,ih) = weidat(7);
  mlat3(nn,num,ih) = weidat(8);
  mlon(nn,num,ih) = weidat(9);
  alat(nn,num,ih) = abs(weidat(8));
  mlt3(nn,num,ih)  = weidat(10);
  eflx(nn,num,ih) = weidat(11);
  ekev(nn,num,ih) = weidat(13)*0.001;
  iflx(nn,num,ih) = weidat(14);
  ikev(nn,num,ih) = weidat(16)*0.001;
% Calculating Particle Induced Conductivity 
	condp(nn,num,ih) = 40 * (ekev(nn,num,ih)) * (eflx(nn,num,ih))^(1/2) / (16 + (ekev(nn,num,ih))^2);
	condh(nn,num,ih) = 0.45 * (ekev(nn,num,ih))^0.85 * condp(nn,num,ih);
% Calculating conductivities from UV/EUV
% Calls in variables from subsol.m and cossza.m
%   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
	[slat slong] = subsol(iyearN, ndayN, uts3(nn,num,ih)/3600, uts3(nn,num,ih)/60, uts3(nn,num,ih));
	csza = cossza(glat(nn,num,ih), glon(nn,num,ih), slat, slong);	
	if csza > cos_85
		condUVp(nn,num,ih) = 0.88 * (Sa * csza)^(1/2) + 0.1; 
		condUVh(nn,num,ih) = 1.5 * (Sa * csza)^(1/2) +0.2;
	else 
		f  = 0.06 + exp(1.803 * tanh(3.833 * csza) + (0.5 * csza) - 2.332);
		Bp = (1 - (0.99524 * (sin(slat))^2)) * (1 + (0.3 * cos(slat)^2)); 
		condUVp(nn,num,ih) = 12.5 * (Sa/180)^(1/2) * f / Bp;
	end		
	condp_Tot(nn,num,ih) = (condp(nn,num,ih)^2 + condUVp(nn,num,ih)^2)^(1/2);	

  if (eflx(nn,num,ih)<fxmin)
% DO NOT SET ekev=0 because this affects the smoothed values! - fix this later for
%  ekev(nn,num,ih) = 0.;
   eflx(nn,num,ih) = 0.;
%   ekev(nn,num,ih) = NaN;
%   eflx(nn,num,ih) = NaN;
   condp(nn,num,ih) = 0.;
   condh(nn,num,ih) = 0.;
   condp_Tot(nn,num,ih) = condUVp(nn,num,ih);
  end
  if (iflx(nn,num,ih)<fxmini)
%   ikev(nn,num,ih) = 0.;
   ikev(nn,num,ih) = NaN;
   iflx(nn,num,ih) = NaN;
  end

 end	% if abs(mlat)>45
end	% do n=1:nrd

% Save array sizes (biggest nn, biggest numh, 2)
[nnx,numx,ih2]=size(eflx)



%----------------------Tosses sharp peaks in Eflux data-------------------------
% Tossing peaks greater than 8x unsmoothed mean
 
numbad = 0;
numgood = 0;
for ih=1:2 
 for num=1:numh(ih)
  nn = nnsav(num,ih); 
 meanflx = mean(eflx(1:nn,num,ih))
 maxeflx = max(eflx(1:nn,num,ih))
 if (maxeflx>tossx*meanflx)
'max>tossx*mean!'
 end
   for n=2:nn

	if eflx(n,num,ih) < tossx*meanflx
		mflx(n,num,ih) = eflx(n,num,ih);
	        numgood = numgood + 1;
	else 
num,ih
% Fix with interpolation
%		mflx(n,num,ih) = NaN;
	if (eflx(n+1,num,ih) < tossx*meanflx)
	        mflx(n,num,ih) = 0.5*(eflx(n-1,num,ih)+eflx(n+1,num,ih));
                numbad = numbad + 1
n,eflx(n,num,ih),mflx(n,num,ih)
        else
	   for m=nn:-1:n+2
	     if (eflx(m,num,ih)<tossx*meanflx)
	      igood = m;
	     end
 	   end
	   for m=n:igood-1
	     mflx(m,num,ih) = ( (igood-m)*eflx(n-1,num,ih) + (m-n+1)*eflx(igood,num,ih) )/(igood-n+1);
             numbad = numbad + 1
m,eflx(m,num,ih),mflx(m,num,ih)
	   end
        end
	end	% end of good or interpolate

   end 
 end	% for num=1:numh(ih)
end	% for ih=1:2
numbad,numgood
% ratiobad ~0.0003 or ~0.03% or 14 in 45,685 for 05002 with tossx=100 * meanflx
ratiobad=numbad/numgood

% Linearly interpolates gaps in time series 
% R. Pawlowicz 6/Nov/99

mflx2=mflx;

%bd=isnan(mflx);
%gd=find(~bd);

% This section takes ~3 min on puffin so interpolate earlier!!

% This sets bads from NaN to 0 to then interpolate??
%bd([1:(min(gd)-1) (max(gd)+1):end])=0;

%if isnan(mflx) == 0
%mflx2(bd)=interp1(gd,mflx(gd),find(bd));
%end

% 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
% Smooths the new eflx set of values with the tossed peaks > 8x mean 
% Smooths the mean energy <E> KeV 


%  This turns a 3-parameter array of (nnx,numx,2) into a single parameter array nnx*numx*2
aveflx = smooth( mflx2,nsm); 
emax = fix(max(aveflx) + 0.99)
mekev  = smooth( ekev, nsm);
length(aveflx)

% Recalc avecondp (w EUV) from smoothed emax,mekev

% Filtering bad data on Ekev plot
% Do not plot E-kev if E-flux less than fxminkev
for ih=1:2 
 for num=1:numh(ih)
  nn = nnsav(num,ih); 
   for n=1:nn
     m = n+(num-1)*nnx + (ih-1)*numx*nnx;
	avcondp = 40 * (mekev(m)) * (aveflx(m))^(1/2) / (16 + (mekev(m))^2);
        avecondp(n,num,ih) = (avcondp^2 + condUVp(n,num,ih)^2)^(1/2);
	if aveflx(m) <= fxminkev
		ekev(n,num,ih)  = 0;
		mekev(m) = 0;
	end
   end 
 end	% for num=1:numh(ih)
end	% for ih=1:2

% Calculate median mean of whole day (meanflx=mean(eflx) gives array of meanflx(numx,2) including zeros)
mnflxnumx2 = mean(eflx);
% median(mnflxnumx2) gives 2 medians for SH and NH
medav = median(mnflxnumx2)
meanflux = 0.5*(medav(1)+medav(2))

pmaxnumx2 = max(avecondp)
pmax2 = max(pmaxnumx2)
pmax = fix(max(pmax2)+0.99)

% Open write file of calculated smoothed peaks, 1/e fall-offs, dmlat, etc
%  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.3d',iyearN,ndayN);
CALC = fopen(filcalc,'w');
fprintf (CALC,'%4.4d %3.3d %2.2d %2.2d %4.2f %2.2d %3.3d %6.2f %s year daynum numSH numNH fxmin, nsec smoothing, tossx, meanflux (delete spikes if >meanflux*tossx) filcalc\n',iyearN,ndayN,numh(1),numh(2),fxmin,nsm,tossx,meanflux,filcalc);

% Open files of smoothed ekev and eflux between 1/e fall-offs and peaks
fout1 = sprintf('/hao/aim3/dmsp/dmsp_rich/txt_dmspssj/%4.4d%3.3deq06.txt',iyearN,ndayN);
HPOUT1 = fopen( fout1, 'w' );
fprintf( HPOUT1, '%4.0d %3.3d %4.2f %2.2d %3.3d %6.2f %s year daynum fxmin, nsec smoothing, tossx, meanflux (delete spikes if >meanflux*tossx) file06 \n',iyearN,ndayN,fxmin,nsm,tossx,meanflux,fout1);
fprintf (HPOUT1,'nda UTsec  mW/m3    MLT   mlat    keV\n');
fout2 = sprintf('/hao/aim3/dmsp/dmsp_rich/txt_dmspssj/%4.4d%3.3deq18.txt',iyearN,ndayN);
HPOUT2 = fopen( fout2, 'w' );
fprintf( HPOUT2, '%4.0d %3.3d %4.2f %2.2d %3.3d %6.2f %s year daynum fxmin, nsec smoothing, tossx, meanflux (delete spikes if >meanflux*tossx) file18 \n',iyearN,ndayN,fxmin,nsm,tossx,meanflux,fout2);
fprintf (HPOUT2,'nda UTsec  mW/m3    MLT   mlat    keV\n');
fout3 = sprintf('/hao/aim3/dmsp/dmsp_rich/txt_dmspssj/%4.4d%3.3dpo06.txt',iyearN,ndayN);
HPOUT3 = fopen( fout3, 'w' );
fprintf( HPOUT3, '%4.0d %3.3d %4.2f %2.2d %3.3d %6.2f %s year daynum fxmin, nsec smoothing, tossx, meanflux (delete spikes if >meanflux*tossx) filp06 \n',iyearN,ndayN,fxmin,nsm,tossx,meanflux,fout3);
fprintf (HPOUT3,'nda UTsec  mW/m3    MLT   mlat    keV\n');
fout4 = sprintf('/hao/aim3/dmsp/dmsp_rich/txt_dmspssj/%4.4d%3.3dpo18.txt',iyearN,ndayN);
HPOUT4 = fopen( fout4, 'w' );
fprintf( HPOUT4, '%4.0d %3.3d %4.2f %2.2d %3.3d %6.2f %s year daynum fxmin, nsec smoothing, tossx, meanflux (delete spikes if >meanflux*tossx) filp18 \n',iyearN,ndayN,fxmin,nsm,tossx,meanflux,fout4);
fprintf (HPOUT4,'nda UTsec  mW/m3    MLT   mlat    keV\n');

%-------------------------Plot all tracks in SH and NH----------------------
% 

%  Plot real data or smooth data

if (ifsm~=-1)
for ih=1:2 
ih,numh(ih)
  hemisphere = 'N';
  if (ih==1)
    hemisphere = 'S';
  end
  if ifsm==1
PlotTitle = ([iyear, ' ', imonth,' ',iday,' ',hemisphere,'H DMSP-F', inst,' smoothed SSJ4'] );
  else
PlotTitle = ([iyear, ' ', imonth,' ',iday,' ',hemisphere,'H DMSP-F', inst,' original SSJ4'] );
  end

for num=1:numh(ih)
% Will label start of pass in hours=uts/3600.
% data1 = [uts/3600. alat mlt eflx]; 
 nn = nnsav(num,ih); 
for n=1:nn
     m = n+(num-1)*nnx + (ih-1)*numx*nnx;
 data2(n,1) = uts3(n,num,ih)/3600.;
 data2(n,2) = mlat3(n,num,ih);
 data2(n,3) = mlt3(n,num,ih);
 if (ifsm==1)
  data2(n,4) = aveflx(m);
 else
  data2(n,4) = eflx(n,num,ih);
 end
end
%data2 = [uthp(:,1) mlatp(:,1) mltp(:,1) abs(poyf(:,1))];	% works!
 if (num==1)
   [out,track,lines] = satelliteTrackPlot(data2,[0 cbarmax], cbarmax+5, PlotTitle, 1,cbarmax,{});
 else
   [out,track,lines] = satelliteTrackPlot(data2,[0 cbarmax], cbarmax+5, PlotTitle, 0,cbarmax,{});
 end
end		% for num=1:numh(ih)
hold off;
end		% for ih=1:2
end		% for if ifsm~=-1

% ------------------Finding Peaks and 1/e Fall Off in Eflux------------------

% Using the index of the maximum mlat (ixmlat) to find the maximum electron fluxes on either side of this index

% Define also the 1-D uts and mlat for later use

for ih=1:2 
 for num=1:numh(ih)
  nn = nnsav(num,ih);
  xmlat = minmlat;
  ixmlat = 0; 
  n1 = 1+(num-1)*nnx + (ih-1)*numx*nnx;
  nnn = nn+(num-1)*nnx + (ih-1)*numx*nnx;
   for n=1:nn
     m = n+(num-1)*nnx + (ih-1)*numx*nnx;
     mlat(m) = mlat3(n,num,ih);
     uts(m) = uts3(n,num,ih);
     mlt(m) = mlt3(n,num,ih);
     eflx0(m) = eflx(n,num,ih);
    if (abs(mlat3(n,num,ih)) > xmlat)
     xmlat = abs(mlat3(n,num,ih));
     ixmlat = m;
    end
   end

% Skip polar pass if ixmlat==0
ifskip = 0;
if (ixmlat==0)
 ifskip = 1;
end

if (ixmlat>0)
% Use the max funtion to find the two peaks
[peak1 ipeak1] = max(aveflx(n1:ixmlat));
ipeak1 = ipeak1 + n1 - 1;
% Skip polar pass if ipeak1=n1 or ixmlat
if (ipeak1==n1 | ipeak1==ixmlat)
 ifskip = 1;
end
[peak2 ipeak2] = max(aveflx(ixmlat:nnn));
ipeak2 = ipeak2 + ixmlat; 
% Skip polar pass if ipeak2=nn or ixmlat
if (ipeak2==nn | ipeak2==ixmlat)
 ifskip = 1;
end

peakarr = [peak1 peak2]
% Skip polar pass if peaks <eflxmin=2.718*fxmin
if (peak1<eflxmin)
'peak1<eflxmin'
 ifskip = 1;
end
if (peak2<eflxmin)
'peak2<eflxmin'
 ifskip = 1;
end
ipeakarr = [ipeak1 ipeak2];

np = length(peakarr);
end		% if ixmlat>0 (ie, mlat>minmlat (55deg))

% If skimming pass (ixmlat<minmlat) or if peak1<eflxmin or peak2<eflxmin then skip
%  NOTE:  fxmin=0.04 (very low), but eflxmin=0.5 (rel high), so reset eflxmin=fxmin*2.718
if (ifskip==0)

% Calculate 1/e fall off for each peak

peakfall = peakarr(1:np)/2.718;

% Finding the 1/e fall off value before and after each peak (4 total values) closest to equator and poles

for i = n1:ipeak1
	if aveflx(i) > peakfall(1) 
		ifall1 = i-1;  
		break
	end
end

for i = ixmlat:-1:ipeak1
	if aveflx(i) > peakfall(1) 
		ifall2 = i-1;  
		break
	end
end

for i = ixmlat:ipeak2
	if aveflx(i) > peakfall(2) 
		ifall3 = i-1;  
		break
	end
end

for i = nnn:-1:ipeak2
	if aveflx(i) > peakfall(2) 
		ifall4 = i;  
		break
	end
end

fall = [ifall1 ifall2 ifall3 ifall4];

% Finding endpoints for Gaussian curve (appx. where peak approaches zero)

for i = n1:nnn
	if i < fall(1) && i >= n1
		if aveflx(i) < peakarr(1)*.05 && aveflx(i) 
		istart(1) = i;
		start(1) = aveflx(i);
		end
	elseif i >= fall(2) && i < ixmlat
		[start(2) istart(2)] = min(aveflx(ifall2:ixmlat));
		istart(2) = istart(2) + fall(2);
	elseif i >= ixmlat && i < fall(3)
		[start(3) istart(3)] = min(aveflx(ixmlat:ifall3));
		istart(3) = istart(3) + ixmlat;
	elseif i >= fall(4) && i <= nnn
		if aveflx(i) < peakarr(2)*.05 
		istart(4) = i;
		start(4) = aveflx(i)
		break
		end
	end
end		

istart

% Calculates the Gaussian curves 
% Calculates the summation of each of the four half-widths from both the Gaussian and original EFlux curve
% Calculated the MLT delta h
%      m = n+(num-1)*nnx + (ih-1)*numx*nnx;

gaussum = [0 0 0 0];
nsmsum = [0 0 0 0];
origsum = [0 0 0 0];

	delh(1)      = abs(mlat(ipeakarr(1))-mlat(fall(1)));
for i = ifall1:ipeak1
	gaussflx1(i) = peakarr(1) * exp( -(abs(mlat(ipeakarr(1))-mlat(i)))^2 / delh(1)^2 );
	gaussum(1)   = gaussum(1) + gaussflx1(i);
	nsmsum(1)    = nsmsum(1) + aveflx(i);
        nv = i - (num-1)*nnx - (ih-1)*numx*nnx;
        origsum(1)   = origsum(1) + eflx(nv,num,ih);
end

	delh(2)      = abs(mlat(ipeakarr(1))-mlat(fall(2)));
for i = ipeak1:ifall2
	gaussflx2(i) = peakarr(1) * exp( -(abs(mlat(ipeakarr(1))-mlat(i)))^2 / delh(2)^2 );
	gaussum(2)   = gaussum(2) + gaussflx2(i);
	nsmsum(2)    = nsmsum(2) + aveflx(i);
        origsum(2)   = origsum(2) + eflx0(i);
end

	delh(3)      = abs(mlat(ipeakarr(2))-mlat(fall(3)));
for i = ifall3:ipeak2
	gaussflx3(i) = peakarr(2) * exp( -(abs(mlat(ipeakarr(2))-mlat(i)))^2 / delh(3)^2 );
	gaussum(3)   = gaussum(3) + gaussflx3(i);
	nsmsum(3)    = nsmsum(3) + aveflx(i);
        origsum(3)   = origsum(3) + eflx0(i);
end

	delh(4)      = abs(mlat(ipeakarr(2))-mlat(fall(4)));
for i = ipeak2:ifall4
	gaussflx4(i) = peakarr(2) * exp( -(abs(mlat(ipeakarr(2))-mlat(i)))^2 / delh(4)^2 );
	gaussum(4)   = gaussum(4) + gaussflx4(i);
	nsmsum(4)    = nsmsum(4) + aveflx(i);
        origsum(4)   = origsum(4) + eflx0(i);
end

delh
delmax=max(delh(1:4))+.5;

% Add additional flux from 1/e to where flux is zero
gaussum0 = gaussum;
origsum0 = origsum;
nsmsum0 = nsmsum;

for i = istart(1):ifall1-1	
	gaussflx1(i) = peakarr(1) * exp( -(abs(mlat(ipeakarr(1))-mlat(i)))^2 / delh(1)^2 );
	gaussum0(1)   = gaussum0(1) + gaussflx1(i);
	nsmsum0(1)    = nsmsum0(1) + aveflx(i);
        origsum0(1)   = origsum0(1) + eflx0(i);
end

for i = ifall2+1:istart(2)
	gaussflx2(i) = peakarr(1) * exp( -(abs(mlat(ipeakarr(1))-mlat(i)))^2 / delh(2)^2 );
	gaussum0(2)   = gaussum0(2) + gaussflx2(i);
	nsmsum0(2)    = nsmsum0(2) + aveflx(i);
        origsum0(2)   = origsum0(2) + eflx0(i);
end

for i = istart(3):ifall3-1
	gaussflx3(i) = peakarr(2) * exp( -(abs(mlat(ipeakarr(2))-mlat(i)))^2 / delh(3)^2 );
	gaussum0(3)   = gaussum0(3) + gaussflx3(i);
	nsmsum0(3)    = nsmsum0(3) + aveflx(i);
        origsum0(3)   = origsum0(3) + eflx0(i);
end

for i = ifall4+1:istart(4)
	gaussflx4(i) = peakarr(2) * exp( -(abs(mlat(ipeakarr(2))-mlat(i)))^2 / delh(4)^2 );
	gaussum0(4)   = gaussum0(4) + gaussflx4(i);
	nsmsum0(4)    = nsmsum0(4) + aveflx(i);
        origsum0(4)   = origsum0(4) + eflx0(i);
end

% Calculate Gaussian sum to origianl sum ratio
% 05002 47032UTsec NH 88.1mlat have peak and istart flux at ixmlat, but delh(2)=0.05 because 1/e flux set to 47031UT (should have been LATER UT than peak!!)
%  Thus, origsum,nsmsum,fitsum=0 (but suffix 0 goes out 1 so is non-zero)
%  ie, sometimes, cannot get good peak info on 1 side - need a flag? (good, bad?)

ratiof2o(1:4) = gaussum ./ origsum
ratiof2s(1:4) = gaussum ./ nsmsum
ratios2o(1:4) = nsmsum ./ origsum
ratiof2o0(1:4) = gaussum0 ./ origsum0
ratiof2s0(1:4) = gaussum0 ./ nsmsum0
ratios2o0(1:4) = nsmsum0 ./ origsum0

% Print out calculations for later plots vs MLT 
%  Print out first 6 flags (delh06eq peak06 delh06pol delh18pol peak18 peak18eq)
%    (0=good, 1=bad - default 0) and highest mlat info (must be >minmlat or 55deg) 
fprintf (CALC,'0 0 0 0 0 0 flags (0=good, 1=bad) for delh06eq peak06 delh06pol delh18pol peak18 peak18eq\n')
fprintf (CALC,'%2.2d %1.1d %5.0d %5.0d %6.2f %6.2f %6.2f num_orb ih(1=SH,2=NH) ixmlat,uts,mlat,mlt,aveflx \n',num,ih,ixmlat,uts(ixmlat),mlat(ixmlat),mlt(ixmlat),aveflx(ixmlat));
% Always list in order: (1) 06MLT eq start fall peak, (2) 06MLT pol start fall,
%  (3) 18MLT pol start fall, and (4) 18MLT eq start fall peak 
%  NOTE:  This won't be in UT order, even for the SH from 6-18MLT!
if (ih==1)
 i=1;
fprintf (CALC,'%6.3f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %5.0d %5.0d %6.2f %6.2f %6.2f %5.0d %5.0d %6.2f %6.2f %6.2f %5.0d %5.0d %6.2f %6.2f %6.2f %6.2f\n',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);
 i=2;
fprintf (CALC,'%6.3f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %5.0d %5.0d %6.2f %6.2f %6.2f %5.0d %5.0d %6.2f %6.2f %6.2f\n',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)));
 i=3;
fprintf (CALC,'%6.3f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %5.0d %5.0d %6.2f %6.2f %6.2f %5.0d %5.0d %6.2f %6.2f %6.2f\n',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)));
 i=4;
fprintf (CALC,'%6.3f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %5.0d %5.0d %6.2f %6.2f %6.2f %5.0d %5.0d %6.2f %6.2f %6.2f %5.0d %5.0d %6.2f %6.2f %6.2f %6.2f\n',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)),ipeak2,uts(ipeak2),mlat(ipeak2),mlt(ipeak2),aveflx(ipeak2),peak2);
else
 i=4;
fprintf (CALC,'%6.3f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %5.0d %5.0d %6.2f %6.2f %6.2f %5.0d %5.0d %6.2f %6.2f %6.2f %5.0d %5.0d %6.2f %6.2f %6.2f %6.2f\n',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)),ipeak2,uts(ipeak2),mlat(ipeak2),mlt(ipeak2),aveflx(ipeak2),peak2);
 i=3;
fprintf (CALC,'%6.3f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %5.0d %5.0d %6.2f %6.2f %6.2f %5.0d %5.0d %6.2f %6.2f %6.2f\n',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)));
 i=2;
fprintf (CALC,'%6.3f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %5.0d %5.0d %6.2f %6.2f %6.2f %5.0d %5.0d %6.2f %6.2f %6.2f\n',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)));
 i=1;
fprintf (CALC,'%6.3f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %5.0d %5.0d %6.2f %6.2f %6.2f %5.0d %5.0d %6.2f %6.2f %6.2f %5.0d %5.0d %6.2f %6.2f %6.2f %6.2f\n',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);
end

%----------------------- Store smoothed ekeV and eflux in peaks to fall-off -----

je06 = 0;
je18 = 0;
jp06 = 0;
jp18 = 0;

for i = ifall1:ipeak1
    if aveflx(i)>aveflx(ifall1)
	if mlt(i) >= 0 && mlt(i) < 12 
	   je06 = je06 + 1;
		utse6(je06)    = uts(i);
		avflxe6(je06)  = aveflx(i); 
		mlte6(je06)    = mlt(i);
		mlate6(je06)   = mlat(i);
		ekeve6(je06)   = mekev(i);	
	elseif mlt(i) >= 12 && mlt(i) <= 24
	   je18 = je18 + 1;
		utse18(je18)   = uts(i);
		avflxe18(je18) = aveflx(i);
		mlte18(je18)   = mlt(i);
		mlate18(je18)  = mlat(i);
		ekeve18(je18)  = mekev(i);
	end
    end
end

for i = ipeak1:ifall2
    if aveflx(i)>aveflx(ifall2)
	if mlt(i) >=0 && mlt(i) < 12
	   jp06 = jp06 + 1;
		utsp6(jp06)   = uts(i);
		avflxp6(jp06) = aveflx(i);
		mltp6(jp06)   = mlt(i);
		mlatp6(jp06)  = mlat(i);
		ekevp6(jp06)  = mekev(i);
	elseif mlt(i) >= 12 && mlt(i) <= 24
	   jp18 = jp18 + 1;
		utsp18(jp18)  = uts(i);
		avflxp18(jp18)= aveflx(i);
		mltp18(jp18)  = mlt(i);
		mlatp18(jp18) = mlat(i);
		ekevp18(jp18) = mekev(i);
	end
    end
end

for i = ifall3:ipeak2	
    if aveflx(i)>aveflx(ifall3)
	if mlt(i) >=0 && mlt(i) < 12
	   jp06 = jp06 + 1;
		utsp6(jp06)   = uts(i);
		avflxp6(jp06) = aveflx(i);
		mltp6(jp06)   = mlt(i);
		mlatp6(jp06)  = mlat(i);
		ekevp6(jp06)  = mekev(i);
	elseif mlt(i) >= 12 && mlt(i) <= 24
	   jp18 = jp18 + 1;
		utsp18(jp18)  = uts(i);
		avflxp18(jp18)= aveflx(i);
		mltp18(jp18)  = mlt(i);
		mlatp18(jp18) = mlat(i);
		ekevp18(jp18) = mekev(i);
	end
    end
end

for i = ipeak2:ifall4
    if aveflx(i)>aveflx(ifall4)
	if mlt(i) >=0 && mlt(i) < 12
	   je06 = je06 + 1;
		utse6(je06)    = uts(i);
		avflxe6(je06)  = aveflx(i); 
		mlte6(je06)    = mlt(i);
		mlate6(je06)   = mlat(i);
		ekeve6(je06)   = mekev(i);	
	elseif mlt(i) >= 12 && mlt(i) <= 24
	   je18 = je18 + 1;
		utse18(je18)   = uts(i);
		avflxe18(je18) = aveflx(i);
		mlte18(je18)   = mlt(i);
		mlate18(je18)  = mlat(i);
		ekeve18(je18)  = mekev(i);
	end
    end
end

for i=1:je06
 fprintf(HPOUT1, '%3.3d %5.0f %6.2f %6.2f %6.2f %6.2f\n',ndayN, utse6(i), avflxe6(i), mlte6(i), mlate6(i), ekeve6(i));
end
for i=1:je18
 fprintf(HPOUT2, '%3.3d %5.0f %6.2f %6.2f %6.2f %6.2f\n',ndayN, utse18(i), avflxe18(i), mlte18(i), mlate18(i), ekeve18(i));
end
for i=1:jp06
 fprintf(HPOUT3, '%3.3d %5.0f %6.2f %6.2f %6.2f %6.2f\n',ndayN, utsp6(i), avflxp6(i), mltp6(i), mlatp6(i), ekevp6(i));
end
for i=1:jp18
 fprintf(HPOUT4, '%3.3d %5.0f %6.2f %6.2f %6.2f %6.2f\n',ndayN, utsp18(i), avflxp18(i), mltp18(i), mlatp18(i), ekevp18(i));
end

%------------------------------ Plot each pass ----------------------------------
% top plot
% Plots E flux smoothed over 19 seconds (red), E flux unsmoothed raw data (cyan), and Gaussian curves (black) 

dark_green = [0 0.5 0];
fsize = 12;
fsiz = 11;
x1=uts(n1)/3600.;
x2=uts(nnn)/3600.;
  hemisphere = 'N';
  if (ih==1)
    hemisphere = 'S';
  end
numS = num2str(num);
PlotTitle = ([numS,' ',iyear, ' ', imonth,' ',iday,' ',hemisphere,'H DMSP-F', inst,' SSJ4'] );

	H1=figure;
		subplot ('Position',[0.1,0.70,0.83,0.20])
	hold on;
	plot(uts3(1:nn,num,ih)/3600.,eflx(1:nn,num,ih),'c');
	plot(uts(n1:nnn)/3600.,aveflx(n1:nnn),'m');
% istart is undefined function!
%	plot(uts(istart(1:2))/3600.,0.1*peakarr(1),'ks');
	plot(uts(ifall1:ipeak1)/3600.,gaussflx1(ifall1:ipeak1),'k', 'LineWidth',2);
	plot(uts(ipeak1:ifall2)/3600.,gaussflx2(ipeak1:ifall2),'k', 'LineWidth',2);


	plot(uts(ifall3:ipeak2)/3600.,gaussflx3(ifall3:ipeak2),'k', 'LineWidth',2);
	plot(uts(ipeak2:ifall4)/3600.,gaussflx4(ipeak2:ifall4),'k', 'LineWidth',2);
	text (x1+0.01,0.65*emax,'e-','Color','c','fontsize',fsize, 'fontweight','bold');
	text (x1+0.01,0.45*emax,'. = 1/e','color','b','fontsize',fsize, 'fontweight','bold');
	text (x1+0.01,0.25*emax,'\diamondsuit = Peak','color','r','fontsize',fsize, 'fontweight','bold');

% Plotting peaks and 1/e fall offs on top of eflx plot

for i = 1:np
		plot(uts(ipeakarr(i))/3600, peakarr(i), 'rd', 'MarkerFaceColor', 'r');	
end	
	plot(uts(fall(1:2))/3600, peakfall(1), 'm.');
	plot(uts(fall(3:4))/3600, peakfall(2), 'm.');

grid;
                        ylabel('E flux (mW/m^2)','fontsize',fsize,'fontweight','bold')
	                title (PlotTitle,'fontsize',fsize,'fontweight','bold');
	set (gca,'Xlim',[x1 x2],'Ylim',[0 20],'XTickLabel',{' ';});
		subplot ('Position',[0.1,0.50,0.83,0.20])
hold on;
	plot(uts3(1:nn,num,ih)/3600.,ekev(1:nn,num,ih),'g');
	plot(uts(n1:nnn)/3600.,mekev(n1:nnn),'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.20])
hold on;
	
	plot(uts3(1:nn,num,ih)/3600., condp(1:nn,num,ih),'r'); 
	plot(uts3(1:nn,num,ih)/3600., condp_Tot(1:nn,num,ih), 'k');
	plot(uts3(1:nn,num,ih)/3600., condUVp(1:nn,num,ih), 'g');
	plot(uts3(1:nn,num,ih)/3600., avecondp(1:nn,num,ih), 'm');

	text (x1+0.01,0.8*pmax,'\Sigma_P','Color','r','fontsize',fsize, 'fontweight','bold');
	text (x1+0.01,0.6*pmax,'\Sigma_PUV/EUV','color','g','fontsize',fsize, 'fontweight','bold');
	text (x1+0.01,0.4*pmax,'\Sigma_PTot','color','k','fontsize',fsize, 'fontweight','bold');
 	text (x1+0.01,0.2*pmax,'smooth','color','m','fontsize',fsize, 'fontweight','bold');
       
grid;
        ylabel('\Sigma_P (mho)', 'fontsize',fsiz,'fontweight','bold')
	set (gca,'Xlim',[x1 x2],'Ylim',[0 pmax],'XTickLabel',{' ';});

% bottom plots MLAT and MLT
% Plot MLT 
		subplot ('Position',[0.1,0.20,0.83,0.10])
	hold on;
	plot(uts3(1:nn,num,ih)/3600.,mlt3(1:nn,num,ih),'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(n1:nnn)/3600.,mlat(n1:nnn),'b');

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


hold off;
end		% if ifskip==0
 end		% for num=1:numh(ih)
end		% for ih=1:2
