;;;;;;;;;;;;;;;;;;;;;;
; ace-event-sw.pro ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
; Convert Solar Wind data from the CDAWeb (http://cdaweb.gsfc.nasa.gov/istp_public/)  
; for the "ACE" data source into a format that the LFM can work with.  
; 
; NOTE:  You must first convert the CDA web file using "fixCDAweb.pl".
;
; See documentation on the wiki at:
; https://wiki.ucar.edu/display/LTR/Solar+Wind
;
; Last modified by Peter Schmitt (schmitt<at>ucar.edu) on February 11, 2008
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; MANUALLY ENTER FILENAMES HERE ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Magnetosphere file:
ace_mag_file='AC_H0_MFI_24259.txt.new'
; Solar wind file:
ace_sw_file='AC_H0_SWE_24259.txt.new'

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; You should not need to edit below this line ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

a=' '
; Read magnetosphere data from AC_H0_MFI (16 second time cadence)
close,10
openr,10,ace_mag_file
readf,10,rows_mag,cols_mag
data_mag=fltarr(cols_mag,rows_mag)
readf,10,data_mag
close,10

jd_mag=fltarr(rows_mag)
ut_mag=fltarr(rows_mag)

;
; Read Solar Wind data from AC_H0_SWE (64 second time cadence)
;
close,10
openr,10,ace_sw_file
readf,10,rows_sw,cols_sw
data_sw=fltarr(cols_sw,rows_sw)
readf,10,data_sw
close,10
jd_sw=fltarr(rows_sw)
ut_sw=fltarr(rows_sw)

;
; initialize variables for storing data:
bxgse = ut_sw & bygse = ut_sw & bzgse = ut_sw 
bt = ut_sw
vxgse = ut_sw & vygse = ut_sw & vzgse = ut_sw
acexgse = ut_sw & aceygse = ut_sw & acezgse = ut_sw
den = ut_sw & temp = ut_sw
; B (gsm)
bxgsm=0.0*bxgse & bygsm=0.0*bygse & bzgsm=0.0*bzgse
; B (sm)
bxsm=0.0*bxgse & bysm=0.0*bygse & bzsm=0.0*bzgse
; V (gsm)
vxgsm=0.0*vxgse & vygsm=0.0*vygse & vzgsm=0.0*vzgse
; V (sm)
vxsm=0.0*vxgse & vysm=0.0*vygse & vzsm=0.0*vzgse
; ACE (gsm)
acexgsm=0.0*acexgse & aceygsm=0.0*aceygse & acezgsm=0.0*acezgse
; ACE (sm)
acexsm=0.0*acexgse & aceysm=0.0*aceygse & acezsm=0.0*acezgse


;
; store data
;
startTime_mag=julday( data_mag(1,0), data_mag(0,0), data_mag(2,0), $
                      data_mag(3,0), data_mag(4,0), data_mag(5,0))
endTime_mag=julday( data_mag(1,rows_mag-1), data_mag(0,rows_mag-1), data_mag(2,rows_mag-1), $
                    data_mag(3,rows_mag-1), data_mag(4,rows_mag-1), data_mag(5,rows_mag-1))
jd_mag(*)=julday( data_mag(1,*),data_mag(0,*),data_mag(2,*),data_mag(3,*),data_mag(4,*),data_mag(5,*))-startTime_mag
ut_mag(*)=(data_mag(3,*))*60.0+data_mag(4,*)+data_mag(5,*)/60.0

startTime_sw=julday( data_sw(1,0), data_sw(0,0), data_sw(2,0), $
                     data_sw(3,0), data_sw(4,0), data_sw(5,0))
endTime_sw=julday( data_sw(1,rows_sw-1), data_sw(0,rows_sw-1), data_sw(2,rows_sw-1), $
                   data_sw(3,rows_sw-1), data_sw(4,rows_sw-1), data_sw(5,rows_sw-1))
jd_sw(*)=julday( data_sw(1,*),data_sw(0,*),data_sw(2,*),data_sw(3,*),data_sw(4,*),data_sw(5,*))-startTime_sw
ut_sw(*)=(data_sw(3,*))*60.0+data_sw(4,*)+data_sw(5,*)/60.0

;
; The magnetosphere data (Bx, By, Bz) is stored at a 4-minute cadence
; whereas the solar wind data (rho, T, Vx,Vy,Vz, Px,Py,Pz) is stored
; at a 64-second cadence.
; 
; interpolate the 4-minute data to 64-seconds:
;

bxgse=spline(jd_mag, data_mag(6,*), jd_sw)
bygse=spline(jd_mag, data_mag(7,*), jd_sw)
bzgse=spline(jd_mag, data_mag(8,*), jd_sw)
bt=sqrt(bxgse^2+bygse^2+bzgse^2)

den(*)=data_sw(6,*)
temp(*)=data_sw(7,*)
vxgse(*)=data_sw(8,*)
vygse(*)=data_sw(9,*)
vzgse(*)=data_sw(10,*)
acexgse(*)=data_sw(11,*)
aceygse(*)=data_sw(12,*)
acezgse(*)=data_sw(13,*)

vt=sqrt(vxgse^2+vygse^2+vzgse^2)
cs=sqrt(2*1.38e-23*temp/(den*1.673e-27))/1000.0

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Compute Average Arrival Time for event ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; We will shift the time series by the average arrival time
; avg arrival time = {average position of ACE} / {average velocity}
; |t| = |x| / |v|
; Units are seconds:
;    seconds = { km } / { km/s }

avgX=sqrt( (MEAN(acexgse))^2 + $
	   (MEAN(aceygse))^2 + $
	   (MEAN(acezgse))^2 )

avgV=sqrt( (MEAN(vxgse))^2 + $
	   (MEAN(vygse))^2 + $
	   (MEAN(vzgse))^2 )

avgArrivalTime = avgX/avgV

;;;;;;;;;;;;;;;;;;;;;;;;;
;convert from gse to sm ;
;;;;;;;;;;;;;;;;;;;;;;;;;

tiltangle=0.0*bxgsm
for i = 0L,rows_sw-1 do begin
  iDay = ymd2dn(data_sw(2,i),data_sw(1,i),data_sw(0,i))
  recalc,data_sw(2,i),iDay,data_sw(3,i),data_sw(4,i),data_sw(5,i)
  k=-1
  gsmx=0
  gsmy=0
  gsmz=0
  gsmgse,gsmx,gsmy,gsmz,vxgse(i),vygse(i),vzgse(i),k
  vxgsm(i)=gsmx
  vygsm(i)=gsmy
  vzgsm(i)=gsmz
  k=-1
  smgsm,smx,smy,smz,vxgsm(i),vygsm(i),vzgsm(i),k
  vxsm(i)=smx
  vysm(i)=smy
  vzsm(i)=smz
  k=-1
  gsmx=0
  gsmy=0
  gsmz=0
  gsmgse,gsmx,gsmy,gsmz,bxgse(i),bygse(i),bzgse(i),k
  bxgsm(i)=gsmx
  bygsm(i)=gsmy
  bzgsm(i)=gsmz
  k=-1
  smgsm,smx,smy,smz,bxgsm(i),bygsm(i),bzgsm(i),k
  bxsm(i)=smx
  bysm(i)=smy
  bzsm(i)=smz
  k=1
  smx=0
  smy=0
  smz=1
  smgsm,smx,smy,smz,gsmx,gsmy,gsmz,k
  tiltangle(i)=atan(gsmx,gsmz)
endfor

;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Get start and stop times ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;
print,'First timestamp (MM DD YYYY  HH MM SS): '
print, format='(i,i,i,i,i,f)', data_sw(1,0), data_sw(0,0), data_sw(2,0), $
                              data_sw(3,0), data_sw(4,0), data_sw(5,0)
print,'Last timestamp (MM DD YYYY  HH MM SS): '
print, format='(i,i,i,i,i,f)', data_sw(1,rows_sw-1), data_sw(0,rows_sw-1), data_sw(2,rows_sw-1), $
                               data_sw(3,rows_sw-1), data_sw(4,rows_sw-1), data_sw(5,rows_sw-1)
read,'Enter starting time (in hours, relative to first timestamp): ',starthr
read,'Enter ending time (in hours, relative to first timestamp): ',endhr

; Make sure starthr and endhr are valid
tmpStartTime = julday( data_sw(1,0), data_sw(0,0), data_sw(2,0), mat(3,0)+starthr, data_sw(4,0), data_sw(5,0))
if (tmpStartTime lt startTime_sw) then begin
   starthr = data_sw(3,0)
   print,'Warning:  You entered an invalid time.  Defaulting to', starthr
endif

; Add on the average arrival time (avgArrivalTime), in seconds.
tmpEndTime = julday( data_sw(1,0), data_sw(0,0), data_sw(2,0), mat(3,0)+starthr+endhr, data_sw(4,0), data_sw(5,0)+avgArrivalTime)

if (tmpEndTime gt endTime_sw) then begin     
  ; make room for "avgArrivalTime" offset:
  endhr = (data_sw(0,rows_sw-1)-data_sw(0,0))*24.0 + $
	     (data_sw(3,rows_sw-1)-data_sw(3,0)) + $
	     (data_sw(4,rows_sw-1)-data_sw(4,0))/60.0 + $
	     (data_sw(5,rows_sw-1)-data_sw(5,0)-avgArrivalTime)/(60.0*60.0)
  print,'Warning:  You entered an invalid time.  Defaulting to', endhr
endif   

   ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;  
   ; Get both plasma and magnetic data at 1 min intervals ;
   ; (i.e. interpolate)                                   ;
   ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;  
   print,"Interpolating to 1 min interval data"

   finaltime=findgen((endhr-starthr)*60.0)+starthr*60.0+1.0
   ; Set the final julian date time.  
   ; offset by avgArrivalTime
   finaltime_jd = julday( data_sw(1,0),data_sw(0,0),data_sw(2,0), $
			  data_sw(3,0),finaltime,data_sw(5,0)+avgArrivalTime)-startTime_sw

   finalvxgsm=spline(jd_sw,vxgsm,finaltime_jd)
   finalvygsm=spline(jd_sw,vygsm,finaltime_jd)
   finalvzgsm=spline(jd_sw,vzgsm,finaltime_jd)
   finalvxsm=spline(jd_sw,vxsm,finaltime_jd)
   finalvysm=spline(jd_sw,vysm,finaltime_jd)
   finalvzsm=spline(jd_sw,vzsm,finaltime_jd)
   finalden=spline(jd_sw,den,finaltime_jd)
   finalcs=spline(jd_sw,cs,finaltime_jd)
   finalbxgsm=spline(jd_sw,bxgsm,finaltime_jd)
   finalbygsm=spline(jd_sw,bygsm,finaltime_jd)
   finalbzgsm=spline(jd_sw,bzgsm,finaltime_jd)
   finalbxsm=spline(jd_sw,bxsm,finaltime_jd)
   finalbysm=spline(jd_sw,bysm,finaltime_jd)
   finalbzsm=spline(jd_sw,bzsm,finaltime_jd)
   finalBt=sqrt(finalbxgsm^2+finalbygsm^2+finalbzgsm^2)
   finalVt=sqrt(finalvxgsm^2+finalvygsm^2+finalvzgsm^2)
   finalVa=1.0e-5*2.2e6*finalBt/sqrt(finalden)   
   finaltilt=spline(jd_sw,tiltangle,finaltime_jd)

   ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   ; now plot the SW Data in the standard form ;
   ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   print,'Now plotting information (press enter)'
   read,a
   ; plot a visual failsafe check on data
   print,"Time interpolation"
   plot,finaltime_jd,xtitle='Time',ytitle='Time (days)',$
       background=255,color=0
   read,a
   print,"X velocity (km/s)"
   plot,finaltime_jd,finalvxgsm,title='X velo (km/s); black=gsm; red=sm',$
       xtitle='Time (days)',ytitle='Vx (km/s)',$
       background=255,color=0
   oplot,finaltime_jd,finalvxsm,color=1
   read,a
   print,"Y velocity (km/s)"
   plot,finaltime_jd,finalvygsm,title='Y velo (km/s); black=gsm; red=sm',$
       xtitle='Time (days)',ytitle='Vy (km/s)',$
       background=255,color=0
   oplot,finaltime_jd,finalvysm,color=1
   read,a
   print,"Z velocity (km/s)"
   plot,finaltime_jd,finalvzgsm,title='Z velo (km/s); black=gsm; red=sm',$
       xtitle='Time (days)',ytitle='Vz (km/s)',$
       background=255,color=0
   oplot,finaltime_jd,finalvzsm,color=1
   read,a
   print,"Density (#/cm^3)"
   plot,finaltime_jd,finalden,title='Density (#/cm^3)',$
       xtitle='Time (days)',ytitle='Density (#/cm^3)',$
       background=255,color=0
   read,a
   print,"Thermal speed (km/s)"
   plot,finaltime_jd,finalcs,title='Thermal Speed (km/s)',$
       xtitle='Time (days)',ytitle='Thermal Speed (km/s)',$
       background=255,color=0
   read,a
   print,"B_x field (nT)"
   plot,finaltime_jd,finalbxgsm,title='X B field (nT); black=gsm; red=sm',$
       xtitle='Time (days)',ytitle='Bx (nT)',$
       background=255,color=0
   oplot,finaltime_jd,finalbxsm,color=1
   read,a
   print,"B_y field (nT)"
   plot,finaltime_jd,finalbygsm,title='Y B field (nT); black=gsm; red=sm',$
       xtitle='Time (days)',ytitle='By (nT)',$
       background=255,color=0
   oplot,finaltime_jd,finalbysm,color=1
   read,a
   print,"B_z field (nT)"
   plot,finaltime_jd,finalbzgsm,title='Z B field (nT); black=gsm; red=sm',$
       xtitle='Time (days)',ytitle='Bz (nT)',$
       background=255,color=0
   oplot,finaltime_jd,finalbzsm,color=1
   read,a
   print,"Total B field (nT)"
   plot,finaltime_jd,finalbt,title='Total B field (nT)',$
       xtitle='Time (days)',ytitle='|B| (nT)',$
       background=255,color=0
   read,a
   print,"Tilt angle"
   plot,finaltime_jd,finaltilt,title='Tilt Angle',$
       background=255,color=0
   read,a

  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  ; attempt to use linear regression fiting to determine Bx(By,Bz) ;
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  Nterms = 2
  Npoints = N_Elements(finalbxgsm)
  x = fltarr(Nterms,Npoints)
  x(0,*) = finalbygsm
  x(1,*) = finalbzgsm
  w = replicate(1.0,Npoints)
  coeffs = regress(x,finalbxgsm,w,newfinalbxgsm,A0,/relative_weight, CHISQ=chi_squared)
  print, '[Bx, By, Bz] fit'
  print, A0, coeffs(0), coeffs(1)
  plot,finaltime_jd,finalbxgsm,title='Bx (black) and Interpolated Bx (green) field',$
       xtitle='Time (days)',ytitle='Bx (nT)',$
       background=255,color=0
  oplot,finaltime_jd,newfinalbxgsm,color=2
  read,a

  print,'Regression:  chi^2 = ', chi_squared

; output data to file

  read,'OK to create SW data file?(y/n) ',a
  if (a eq 'y') then begin
    read, 'Enter file name: ',a
    close,10
    openw,10,a
    smfile=strcompress(a+'.sm',/remove_all)
    length=size(finaltime)
    iDay = ymd2dn(data_sw(2,0),data_sw(1,0),data_sw(0,0))
    close,11
    openw,11,smfile
    printf,10,format='(i5,i4,2i3)',$
           data_sw(2,0),iDay,floor(StartHr),(starthr-floor(starthr))*60
    printf,11,format='(i5,i4,2i3)',$
           data_sw(2,0),iDay,floor(StartHr),(starthr-floor(starthr))*60
    printf,10,format='(i5,3i)',length(1),10 
    printf,11,format='(i5,3i)',length(1),11 
    printf,10," DATA: "
    printf,11," DATA: "
    for i=0,(length(1)-1) do begin
      printf,10,format='(10(f13.2,2x),f6.1,2x,3f6.1,f8.2,4f9.4)', finaltime(i), $
       finalden(i), finalvxgsm(i), finalvygsm(i), finalvzgsm(i), finalcs(i), $
       finalbxgsm(i), finalbygsm(i), finalbzgsm(i), finalbt(i)
      printf,11,format='(10(f13.2,2x),f6.1,2x,3f10.4,f8.2,5f9.4)', finaltime(i), $
       finalden(i), finalvxsm(i), finalvysm(i), finalvzsm(i), finalcs(i), $
       finalbxsm(i), finalbysm(i), finalbzsm(i), finalbt(i),finaltilt(i)
    endfor
    close,10
    close,11

    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    ; plot everything, including regression: ;
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     t0 = min(finaltime_jd)
     t1 = max(finaltime_jd)
     ;
     !p.multi=[0,1,8]
     !x.range=[t0,t1]
     !x.style=1
     !x.minor=6
     !p.charsize=2
     zero1=[t0,t1]
     zero2=[0,0]
     zero3=[180,180]
    
     yMin = min(finalden)
     yMax = max(finalden)

     actualStartTime = julday(data_sw(1,0), data_sw(0,0), data_sw(2,0), starthr,0,0)     
     actualStopTime = julday(data_sw(1,0),data_sw(0,0), data_sw(2,0), starthr+endhr,0,0)

     caldat,actualStartTime, start_month, start_day, start_year
     caldat,actualStopTime, stop_month, stop_day, stop_year

     titleString = string(format='(%"CMIT Wind (GSM) Data for [%d-%d-%d] to [%d-%d-%d]")', $
    			start_month,start_day,start_year, $
    		     	stop_month, stop_day, stop_year)
     plot,zero1,[0,50],ytitle='!17rho', $
       title=titleString,$
        xtitle=' ',xtickname=replicate(' ',8), yticks=3,$
        pos=[0.1,0.8,0.9,0.9],ystyle=1,yrange=[yMin,yMax],/nodata,$
	background=255,color=0
     oplot,finaltime_jd,finalden,linestyle=0,thick=1,color=0
     ;
     yMin=min(finalvxgsm)
     yMin2=min(finalvxsm)
     yMin=min(yMin,yMin2)
     yMax=max(finalvxgsm)
     yMax2=max(finalvxsm)
     yMax=max(yMax,yMax2)
     plot,zero1,[-50,50],ytitle='!17Vx',$
        xtitle=' ',xtickname=replicate(' ',8), yticks=2,$
        pos=[0.1,0.7,0.9,0.795],ystyle=1,yrange=[yMin,yMax],/nodata,$
	color=0
     oplot,finaltime_jd,finalvxgsm,linestyle=0,thick=1,color=0
     ;oplot,finaltime_jd,finalvxsm,linestyle=0,color=1,thick=1
     oplot,zero1,zero2,linestyle=5,color=8
     ;
     yMin=min(finalvygsm)
     yMin2=min(finalvysm)
     yMin=min(yMin,yMin2)
     yMax=max(finalvygsm)
     yMax2=max(finalvysm)
     yMax=max(yMax,yMax2)
     plot,zero1,[-50,50],ytitle='!17Vy',$
        xtitle=' ',xtickname=replicate(' ',8), yticks=2,$
        pos=[0.1,0.6,0.9,0.695],ystyle=1,yrange=[yMin,yMax],/nodata,$
	color=0
     oplot,finaltime_jd,finalvygsm,linestyle=0,thick=1,color=0
     ;oplot,finaltime_jd,finalvysm,linestyle=0,color=1,thick=1
     oplot,zero1,zero2,linestyle=5,color=8
     ;
     yMin=min(finalvzgsm)
     yMin2=min(finalvzsm)
     yMin=min(yMin,yMin2)
     yMax = max(finalvzsm)
     yMax2 = max(finalvzgsm)
     yMax = max(yMax, yMax2)
     plot,zero1,[-50,50],ytitle='!17Vz',$
        xtitle=' ',xtickname=replicate(' ',8), yticks=2,$
        pos=[0.1,0.5,0.9,0.595],ystyle=1,yrange=[yMin,yMax],/nodata,$
	color=0
     oplot,finaltime_jd,finalvzgsm,linestyle=0,thick=1,color=0
     ;oplot,finaltime_jd,finalvzsm,linestyle=0,color=1,thick=1
     oplot,zero1,zero2,linestyle=5,color=8
     ;
     yMin=min(finalcs)
     yMax=max(finalcs)
     plot,zero1,[0,50],ytitle='!17TS', $
        xtitle=' ',xtickname=replicate(' ',8), yticks=2,$
        pos=[0.1,0.4,0.9,0.495],ystyle=1,yrange=[yMin,yMax],/nodata,$
	color=0
     oplot,finaltime_jd,finalcs,linestyle=0,thick=1,color=0
     ;
     yMin=min(finalbxgsm)
     yMin2=min(finalbxsm)
     yMin=min(yMin,yMin2)
     yMax=max(finalbxgsm)
     yMax2=max(finalbxsm)
     yMax=max(yMax,yMax2)
     plot,zero1,[-50,50],ytitle='!17Bx',$
        xtitle=' ',xtickname=replicate(' ',8), yticks=2,$
        pos=[0.1,0.3,0.9,0.395],ystyle=1,yrange=[yMin,yMax],/nodata,$
	color=0
     oplot,finaltime_jd,finalbxgsm,linestyle=0,thick=1,color=0
     oplot,finaltime_jd,newfinalbxgsm,linestyle=0,color=2,thick=1
     ;oplot,finaltime_jd,finalbxsm,linestyle=0,color=1,thick=1
     oplot,zero1,zero2,linestyle=5,color=8
     ;
     yMin=min(finalbygsm)
     yMin2=min(finalbysm)
     yMin=min(yMin,yMin2)
     yMax=max(finalbygsm)
     yMax2=max(finalbysm)
     yMax=max(yMax,yMax2)
     plot,zero1,[-20,20],ytitle='!17By',$
        xtitle=' ',xtickname=replicate(' ',8), yticks=2,$
        pos=[0.1,0.2,0.9,0.295],ystyle=1,yrange=[yMin, yMax],/nodata,$
	color=0
     oplot,finaltime_jd,finalbygsm,linestyle=0,thick=1,color=0
     ;oplot,finaltime_jd,finalbysm,linestyle=0,color=1,thick=1
     oplot,zero1,zero2,linestyle=5,color=8
     ;
     yMin=min(finalbzgsm)
     yMin2=min(finalbzsm)
     yMin=min(yMin,yMin2)
     yMax=max(finalbzgsm)
     yMax2=max(finalbzsm)
     yMax=max(yMax,yMax2)
     plot,zero1,[-10,30],ytitle='!17Bz',xtitle='!17Time',$
        pos=[0.1,0.1,0.9,0.195],ystyle=1,yrange=[yMin, yMax],$
	yticks=2,/nodata,$
	color=0
     oplot,finaltime_jd,finalbzgsm,linestyle=0,thick=1,color=0
     ;oplot,finaltime_jd,finalbzsm,linestyle=0,color=1,thick=1
     oplot,zero1,zero2,linestyle=5,color=8
    
     xyouts,.5,.95,a+'.jpg', alignment=0.5, CHARSIZE=2, /NORMAL
    
     new_Bstr = string(format='(%"<Bx, By, Bz> = <%f, %f, %f>")', $
    		     A0, coeffs(0), coeffs(1))
    
     xyouts,.5,.01, new_Bstr, alignment=0.5, CHARSIZE=1.4, /NORMAL	

     ;;;;;;;;;
     ; FIXME ;
     ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     ; Give the option of saving as a postscript ;
     ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     write_jpeg,a+'.jpg',tvrd(true=3),true=3,quality=100
    
     !p.multi=0
     !x.range=0
     !p.charsize=1
  endif
end