close all

fxmin = 0.25;
mlt24 = 24.;
zero = [0 0];

% BIG Poy NH: 00183 all over; 195,211,261,262,279,331 70-80; 210,218 cap; 311 60-80
% Minor Poy NH: 176, 178, 240, 250, 263, 309, 311?, 313
% 00262 Sep 18:  2UT ~100 75mlat noon; 4UT ~60 75mlat 13mlt; 16UT ~50mW/m2 80mlat noon

% Read SSJ4 .txt file for a day 
hemisphere = 'N';
hemisphere = 'S';

hemisphere = 'S';
% SSJ4 in ~emery/reu11 old print format where txt_dmspssj/year/ have 'real <E>'
%filname = '/home/liamk/data/cenergy/cenergy_all_f15_20050514.txt';
%filname = '/home/emery/reu11/cenergy_all_f15_20050121.txt';
%filname = '/hao/aim3/dmsp/dmsp_rich/txt_dmspssj/2005/15_20050121.txt';
%idmfile = '/hao/aim3/dmsp/dmsp_idm_4s/2005/newl1505021a.txt';
iyear = 2005;
iday = 21;
nskp = 70801;
nrd = 72587-nskp;
nskpp = 71354;
nrdp = 73025-nskpp;
nskpi = 18667;
nrdi = 19086 - nskpi;
x1 = 19.8;
x2 = 20.4;
emax = 17.;
emax = 10.;
poyn = -180.;
dswmax = 60.;		% Dsw ~40-50!

%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;

iday = 262;
%filname = '/home/emery/reu11/cenergy_all_f15_20000918.txt';

nskp = 8116;
nrd = 9627-nskp;
nskpp = 7618;
nrdp = 9232-nskpp;
x1 = 2.1;
x2 = 2.7;

nskp = 12426;
nrd = 13937-nskp;
nskpp = 13768;
nrdp = 15354-nskpp;
x1 = 3.8;
x2 = 4.4;

emax = 11.;
emax = 10.;
imax = 2.;
poyn = -70.;
kmax = 9.;
kmax = 19.;
fxmin = 0.03;
fxmin = 0.15;

iday = 176;	% 15(6),17(8),19(time of FUV, but poor Poy ~4),21(*15mW/m2Poy)
iday = 263;	% Poy 9UT and <4UT right place but FUV 747UT
iday = 195;	% FUV 16UT, big Poy 13,15UT wrong place
iday = 211;     % FUV 1650UT, Poy wrong place
iday = 261;     % FUV 1748UT
fxmin = 0.01;
%filname = '/home/emery/reu11/cenergy_all_f15_20000917.txt';
filname = '/hao/aim3/dmsp/dmsp_rich/txt_dmspssj/2000/15_20000917.txt';
idmfile = '/hao/aim3/dmsp/dmsp_idm_4s/0009/l1500261a.txt';
nskp = 63195;
nrd = 64822-nskp;
nskpp = 63741;
nrdp = 65262-nskpp;
nskpi = 16421.;
nrdi = 16802-nskpi;
x1 = 17.7;
x2 = 18.3;
dswmax = 19.;

poyname = sprintf('/hao/aim3/dmsp/poynting_flux/pfomni_00%3.3d.txt',iday);
%poyname = sprintf('/hao/aim3/dmsp/poynting_flux/pfomni_05%3.3d.txt',iday);

x12 = [x1 x2];

yearS = int2str(iyear);
daynumS = int2str(iday);
PlotTitle = ([yearS,' ',daynumS,' ',hemisphere,'H  DMSP-F15 SSJ4'] );
PoyTitle = ([yearS,' ',daynumS,' ',hemisphere,'H  DMSP-F15 Poynting Flux'] );
xpoy = abs(poyn);
%xpoy = 80.;


%  Changing to array(jj,nj) makes the loop over >~950 go VERY slowly - WHY???

signhem = 1;
if hemisphere=='S'
 signhem = -1;
end

% Read Poynting day 1 sec file with OMNI data
 datpf=dlmread(poyname,'');
 jrdp = length(datpf(:,1));
if (iyear == 2004) 
 numpf = 7;
 numby = 10;
else
 numpf = 6;
 numby = 7;
end
numbz = numby+1;
numvsw = numbz + 1;
numdsw = numvsw + 1;
jj = 0;
nj = 0;
for n=nskpp+1:nskpp+nrdp
% for n=1:jrdp
 if (sign(datpf(n,2))==signhem & abs(datpf(n,2))>45.)
  jj = jj + 1;
  if (jj==1)
   nj = nj + 1;
   jbeg(nj) = n;
  end
  poyf(jj) = datpf(n,numpf);
  by(jj) = datpf(n,numby);
  bz(jj) = datpf(n,numbz);
  vsw(jj) = datpf(n,numvsw);
  dsw(jj) = datpf(n,numdsw);
 uthp(jj) = (datpf(n,1)-iday)*24.;
 mlatp(jj) = datpf(n,2);
 mltp(jj) = datpf(n,3);
 jend(nj) = n;
 else
  if (jj > 0)
   jend(nj) = n-1;
  end
  jj = 0;
 end
end
nj
jbeg
jend
jj = jend(nj)-jbeg(nj) + 1;
uthlen=length(uthp)
utsp = uthp(jj)*3600.
maxpoy = max(poyf);
minpoy = min(poyf);
minmlatp = min(mlatp);
maxmlatp = max(mlatp);

% 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

% Plot all Poynting flux on 1 day for 1 hemisphere
for mp=1:nj
% will label start of pass in hours 
mp
jj = jend(mp) - jbeg(mp) + 1
for n=1:jj
m = jbeg(mp) + n - 1;
data2(n,1) = (datpf(m,1)-iday)*24.;
data2(n,2) = datpf(m,2);
data2(n,3) = datpf(m,3);
data2(n,4) = abs(datpf(m,numpf));
%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!
% if (mp==1)			% use this to plot ALL orbits on 1 plot
 if (mp>0)			% use this to plot each orbit separately
   [out,track,lines] = satelliteTrackPlot(data2,[0 xpoy], 15, PoyTitle, 1,0,{});
 else
   hold on;
   [out,track,lines] = satelliteTrackPlot(data2,[0 xpoy], 15, PoyTitle, 0,0,{});
 end
end

% 6 values (n5) for data5 do not plot correctly in satelliteTrackPlot
%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

%stop

% Read IDM 4s data
 idmdat = dlmread(idmfile,'',0,0);
 ird = length(idmdat(:,1));
ni = 0;
for n=nskpi+1:nskpi+nrdi
 if (sign(idmdat(n,9))==signhem & abs(idmdat(n,9))>45. & idmdat(n,2)<999.)
  ni = ni + 1;
  utsidm(ni) = idmdat(n,1);
  mlatidm(ni) = idmdat(n,9);
  mltidm(ni) = idmdat(n,10);
  vyidm(ni) = idmdat(n,2);
 end
end
ni
minidm = min(vyidm(1:ni))
maxidm = max(vyidm(1:ni))

% Read SSJ4 data - NOTE! have to skip first 2 line description before using dlmread!
 datk=dlmread(filname,'',2,0);
 krd = length(datk(:,1));

nn = 0;
nk = 0;
%SSJ=fopen(filname,'r');
%for n=1:nskp
% topssj = fgetl(SSJ);
%end
%for n=1:nrd
% lindat = fgetl(SSJ);
% weidat = sscanf(lindat,'%f');
%for n=2:krd
for n=nskp+1:nskp+nrd
 if (sign(datk(n,8))==signhem & abs(datk(n,8))>45.)
  nn = nn + 1;
  if (nn==1)
   nk = nk + 1;
   kbeg(nk) = n;
  end
% uts(nn) = weidat(2);
% mlat(nn) = 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;
  uts(nn) = datk(n,2);
  mlat(nn) = datk(n,8);
  mlt(nn) = datk(n,10);
  eflx(nn) = datk(n,11);
  ekev(nn) = datk(n,13)*0.001;
  iflx(nn) = datk(n,14);
  ikev(nn) = datk(n,16)*0.001; 
% new format with 'mean energy' for elec>32eV and ions>68eV
  ekev2(nn) = datk(n,17)*0.001; 
  ikev2(nn) = datk(n,18)*0.001; 
  if (eflx(nn)<fxmin)
   ekev(nn) = 0.;
   ekev2(nn) = 0.;
  end
  if (iflx(nn)<fxmin)
   ikev(nn) = 0.;
   ikev2(nn) = 0.;
  end
 kend(nk) = n;
 else
  if (nn > 0)
   kend(nk) = n-1;
  end
  nn = 0;
 end
end
nk
kbeg
kend
nn = kend(nk)-kbeg(nk) + 1;
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))
  
% Calculate 7sec and 15sec sliding averages for DMSP data in proton spot to cf w WIC/SI12
%  00261 ~18:01-18:04 64884-64996s of 63989-65615s total (#895-1007)
 clear uts7;
 clear avef7;
 clear avif7;
 clear avekev1;
 clear avekev2;
 clear avikev1;
 clear avikev2;
 miss15 = 0;
for nsl=898:1004
 m7 = nsl-897;
 uts7(m7) = uts(nsl);
 avef = 0.;
 avek1 = 0.;
 avek2 = 0.;
 avif = 0.;
 avik1 = 0.;
 avik2 = 0.;
 ne7 = 0;
 ni7 = 0;
 for n7=nsl-3:nsl+3
    avef = avef + eflx(n7);
    if (eflx(n7)>=fxmin)
	ne7 = ne7 + 1;
	avek1 = avek1 + ekev(n7);
        avek2 = avek2 + ekev2(n7);
    end 
 end
 for n7=nsl-7:nsl+7
    avif = avif + iflx(n7);
    if (iflx(n7)>=fxmin)
	ni7 = ni7 + 1;
	avik1 = avik1 + ikev(n7);
        avik2 = avik2 + ikev2(n7);
    end 
 end
 avef7(m7) = avef/7.;
 avif7(m7) = avif/15.;
 avekev1(m7) = avek1/ne7;
 avekev2(m7) = avek2/ne7;
 avkiv1(m7) = 3000.;
 avkiv2(m7) = 3000.;
 if (ni7>0)
  avikev1(m7) = avik1/ni7;
  avikev2(m7) = avik2/ni7;
 else
  miss15 = miss15 + 1;
 end
end
m7
miss15
avef7
avif7
maxfwic = max(avef7)
maxfsi12 = max(avif7)
avefwic = mean(avef7(23:96))
avifsi12 = mean(avif7(23:96))
% Find median from 64910-64983 (23-96) vs 64887-64993 (1-107)
medekev1 = median(avekev1(23:96))
medekev2 = median(avekev2(23:96))
medikev1 = median(avikev1(23:96))
medikev2 = median(avikev2(23:96))

% plot proton/electron spot
SpotTitle = sprintf('Proton/Electron Spot   %s',PlotTitle);
dark_green = [0 0.5 0];
fsize = 12;
fsiz = 11;
emx7 = 9;
kmx7 = 9.9;
x7 = uts7(1);
x8 = uts7(m7);
x78 = [x7 x8];
	H1=figure;
		subplot ('Position',[0.1,0.70,0.83,0.20])
	plot(uts7,avif7*10.,'c');
	hold on;
	plot(uts7,avef7,'g');
	text (x7+5,0.65*emx7,'e-','Color','g','fontsize',fsize, 'fontweight','bold');
	text (x7+5,0.45*emx7,'p+*10.','color','c','fontsize',fsize, 'fontweight','bold');
grid;
                        ylabel('Eflux mW/m^2','fontsize',fsize,'fontweight','bold')
	                title (SpotTitle,'fontsize',fsize,'fontweight','bold');
	set (gca,'Xlim',[x7 x8],'Ylim',[0 emx7],'XTickLabel',{' ';});
		subplot ('Position',[0.1,0.50,0.83,0.20])
	plot(uts7,avikev1,'c');
hold on;
	plot(uts7,avikev2,'Color','k','LineStyle',':');
	plot(uts7,avekev1*10,'g');
	plot(uts7,avekev2*10,'Color',dark_green,'LineStyle','--');
	text (x7+5,0.65*kmx7,'e-*10.','Color','g','fontsize',fsize, 'fontweight','bold');
	text (x7+5,0.45*kmx7,'p+','color','c','fontsize',fsize, 'fontweight','bold');
grid;
         ylabel('<E> keV','fontsize',fsize,'fontweight','bold')
	set (gca,'Xlim',[x7 x8],'Ylim',[0 kmx7],'XTickLabel',{' ';});
% Plot Poynting flux
		subplot ('Position',[0.1,0.40,0.83,0.10])
hold on;
	plot(uthp(1:jj)*3600.,poyf(1:jj),'r');
grid;
box;
         ylabel('S//','fontsize',fsize,'fontweight','bold')
	set (gca,'Xlim',[x7 x8],'Ylim',[poyn 0],'XTickLabel',{' ';});
% Plot 4s IDM Vy in m/s
		subplot ('Position',[0.1,0.30,0.83,0.10])
hold on;
	plot(utsidm(1:ni),vyidm(1:ni),'color',dark_green);
	plot(x78,zero,'--k');
grid;
box;
        ylabel('Vy km/s','fontsize',fsiz,'fontweight','bold')
        xlabel('UT (h)','fontsize',fsize,'fontweight','bold')
	set (gca,'Xlim',[x7 x8],'Ylim',[-0.9 2.9]);
hold off;


% Plot all electron and ion flux on 1 day for 1 hemisphere
EfxTitle = ([yearS,' ',daynumS,' ',hemisphere,'H  DMSP-F15 e- -10*p+ Energy Flux'] );
keVTitle = ([yearS,' ',daynumS,' ',hemisphere,'H  DMSP-F15 e- -p+ Mean KeV Energy'] );
elecTitle = ([yearS,' ',daynumS,' ',hemisphere,'H  DMSP-F15 e- EFlux,<E>keV'] );
protTitle = ([yearS,' ',daynumS,' ',hemisphere,'H  DMSP-F15 p+ EFlux,<E>keV/10'] );
for mp=1:nk
% will label start of pass in hours 
mp
nt = kend(mp) - kbeg(mp) + 1
for n=1:nt
m = kbeg(mp) + n - 1;
% Store energy flux (-10*iflx)
data1(n,1) = datk(m,2)/3600.;
data1(n,2) = datk(m,8);
data1(n,3) = datk(m,10);
data1(n,4) = datk(m,11);
data2(n,1) = datk(m,2)/3600.;
data2(n,2) = datk(m,8);
data2(n,3) = datk(m,10);
data2(n,4) = -10.*datk(m,14);
% Store mean keV energies (-1*ikev)
data3(n,1) = datk(m,2)/3600.;
data3(n,2) = datk(m,8);
data3(n,3) = datk(m,10);
data3(n,4) = datk(m,13)*0.001;
  if (datk(m,11)<fxmin)
   data3(n,4) = 0.;
  end
data4(n,1) = datk(m,2)/3600.;
data4(n,2) = datk(m,8);
data4(n,3) = datk(m,10);
data4(n,4) = -1.*datk(m,16)*0.001;
  if (datk(m,14)<fxmin)
   data4(n,4) = 0.;
  end
end

   [out,track,lines] = satelliteTrackPlot(data1,[0 emax], 15, EfxTitle, 1,emax,{});
   hold on;
   [out,track,lines] = satelliteTrackPlot(data2,[0 emax], 15, EfxTitle, 0,emax,{});
   [out,track,lines] = satelliteTrackPlot(data3,[0 kmax], 15, keVTitle, 1,kmax,{});
   hold on;
   [out,track,lines] = satelliteTrackPlot(data4,[0 kmax], 15, keVTitle, 0,kmax,{});

% Plot electons and protons separately
for n=1:nt
m = kbeg(mp) + n - 1;
% Store electrons
data1(n,1) = datk(m,2)/3600.;
data1(n,2) = datk(m,8);
data1(n,3) = datk(m,10);
data1(n,4) = datk(m,11);
data2(n,1) = datk(m,2)/3600.;
data2(n,2) = datk(m,8);
data2(n,3) = datk(m,10);
data2(n,4) = -1.*datk(m,13)*0.001;
  if (datk(m,11)<fxmin)
   data2(n,4) = 0.;
  end
% Store protons
data3(n,1) = datk(m,2)/3600.;
data3(n,2) = datk(m,8);
data3(n,3) = datk(m,10);
data3(n,4) = datk(m,14);
data4(n,1) = datk(m,2)/3600.;
data4(n,2) = datk(m,8);
data4(n,3) = datk(m,10);
data4(n,4) = -1.*datk(m,16)*0.0001;
  if (datk(m,14)<fxmin)
   data4(n,4) = 0.;
  end
end

   [out,track,lines] = satelliteTrackPlot(data1,[0 emax], 15, elecTitle, 1,emax,{});
   hold on;
   [out,track,lines] = satelliteTrackPlot(data2,[0 emax], 15, elecTitle, 0,emax,{});
   [out,track,lines] = satelliteTrackPlot(data3,[0 imax], 15, protTitle, 1,imax,{});
   hold on;
   [out,track,lines] = satelliteTrackPlot(data4,[0 imax], 15, protTitle, 0,imax,{});

end	% end mp=1:nk


% top plot
%  For 17 Sep 2000 only, plot times of 18:01:08 (64868s) and 18:05:13 (65113s)
image1=[64868 64868];
image2=[65113 65113];
PlotMLTTitle = sprintf('Start %4.1fMLT     %s     %4.1fMLT End',mlt(1),PlotTitle,mlt(nn) );
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 (x1+0.01,0.65*emax,'e-','Color','g','fontsize',fsize, 'fontweight','bold');
	text (x1+0.03,0.65*emax,'p+','color','c','fontsize',fsize, 'fontweight','bold');
%  For 17 Sep 2000 only, plot times of 18:01:08 (64868s) and 18:05:13 (65113s)
yall=[0 emax];
	plot (image1/3600.,yall,'--m');
	plot (image2/3600.,yall,'--m');
grid;
                        ylabel('Eflux mW/m^2','fontsize',fsize,'fontweight','bold')
	                title (PlotMLTTitle,'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.,ikev2(1:nn),'k','LineStyle',':');
	plot(uts(1:nn)/3600.,ekev(1:nn),'g');
	plot(uts(1:nn)/3600.,ekev2(1:nn),':k');
%  For 17 Sep 2000 only, plot times of 18:01:08 (64868s) and 18:05:13 (65113s)
yall=[0 kmax];
	plot (image1/3600.,yall,'--m');
	plot (image2/3600.,yall,'--m');
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(1:jj),poyf(1:jj),'r');
%  For 17 Sep 2000 only, plot times of 18:01:08 (64868s) and 18:05:13 (65113s)
yall=[poyn 0];
	plot (image1/3600.,yall,'--m');
	plot (image2/3600.,yall,'--m');
grid;
box;
         ylabel('S//','fontsize',fsize,'fontweight','bold')
	set (gca,'Xlim',[x1 x2],'Ylim',[poyn 0],'XTickLabel',{' ';});
% Plot 4s IDM Vy in m/s
		subplot ('Position',[0.1,0.30,0.83,0.10])
hold on;
	plot(utsidm(1:ni)/3600.,vyidm(1:ni),'color',dark_green);
	plot(x12,zero,'--k');
%  For 17 Sep 2000 only, plot times of 18:01:08 (64868s) and 18:05:13 (65113s)
yall=[-1.9 2.9];
	plot (image1/3600.,yall,'--m');
	plot (image2/3600.,yall,'--m');
grid;
box;
        ylabel('Vy km/s','fontsize',fsiz,'fontweight','bold')
	set (gca,'Xlim',[x1 x2],'Ylim',[-1.9 2.9],'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(1:jj),dsw(1:jj),'k');
%  For 17 Sep 2000 only, plot times of 18:01:08 (64868s) and 18:05:13 (65113s)
yall=[0 dswmax];
	plot (image1/3600.,yall,'--m');
	plot (image2/3600.,yall,'--m');
grid;
box;
        ylabel('Dsw','fontsize',fsiz,'fontweight','bold')
	set (gca,'Xlim',[x1 x2],'Ylim',[0 dswmax],'XTickLabel',{' ';});
% Plot 5-min IMF By and Bz (old MLAT plot)
		subplot ('Position',[0.1,0.10,0.83,0.10])
	hold on;
%	plot(uts(1:nn)/3600.,mlat(1:nn),'b');
%        plot(uthp(1:jj),mlatp(1:jj),'k');
%grid;
%        ylabel('MLAT','fontsize',fsiz,'fontweight','bold')
        xlabel('UT (h)','fontsize',fsize,'fontweight','bold')
%	set (gca,'Xlim',[x1 x2]);
	plot(uthp(1:jj),bz(1:jj),'r');
        plot(uthp(1:jj),by(1:jj),'color',dark_green);
	text (x1+0.01,10.,'y','Color',dark_green,'fontsize',fsize, 'fontweight','bold');
	text (x1+0.01,-10.,'z','Color','r','fontsize',fsize, 'fontweight','bold');
	plot(x12,zero,'--k');
%  For 17 Sep 2000 only, plot times of 18:01:08 (64868s) and 18:05:13 (65113s)
yall=[-29 29];
	plot (image1/3600.,yall,'--m');
	plot (image2/3600.,yall,'--m');
grid;
box;
        ylabel('By,z','fontsize',fsiz,'fontweight','bold')
	set (gca,'Xlim',[x1 x2],'Ylim',[-29 29]);

hold off;


