;;;;;;;;;;;;;;;;;;;;;;
; omni-event-sw.pro ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
; Convert Solar Wind data from the CDAWeb (http://cdaweb.gsfc.nasa.gov/istp_public/)  
; for the "OMNI" 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 January 25, 2008
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; MANUALLY ENTER FILENAME HERE ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
omnifile='data/OMNI_HRO_1MIN_28057.txt.new'

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; You should not need to edit below this line ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
a=' ' 

close,10
openr,10,omnifile
readf,10,rows,cols
mat=fltarr(cols,rows)
readf,10,mat
close,10

jd=fltarr(rows)
ut=fltarr(rows)
bxgse = ut & bygse = ut & bzgse = ut & bt = ut
vxgse = ut & vygse = ut & vzgse = ut 
den = ut & temp = ut
bxgsm=0.0*bxgse
bygsm=0.0*bygse
bzgsm=0.0*bzgse
bxsm=0.0*bxgse
bysm=0.0*bygse
bzsm=0.0*bzgse
vxgsm=0.0*vxgse
vygsm=0.0*vygse
vzgsm=0.0*vzgse
vxsm=0.0*vxgse
vysm=0.0*vygse
vzsm=0.0*vzgse

startTime=julday( mat(1,0), mat(0,0), mat(2,0), $
		 mat(3,0), mat(4,0), mat(5,0))
endTime=julday( mat(1,rows-1), mat(0,rows-1), mat(2,rows-1), $
	       mat(3,rows-1), mat(4,rows-1), mat(5,rows-1))
jd(*)=julday( mat(1,*),mat(0,*),mat(2,*),mat(3,*),mat(4,*),mat(5,*))-startTime

ut(*)=(mat(3,*))*60.0+mat(4,*)+mat(5,*)/60.0
bxgse(*)=mat(6,*)
bygse(*)=mat(7,*)
bzgse(*)=mat(8,*)
vxgse(*)=mat(9,*)
vygse(*)=mat(10,*)
vzgse(*)=mat(11,*)
den(*)=mat(12,*)
temp(*)=mat(13,*)
bt=sqrt(bxgsm^2+bygsm^2+bzgsm^2)
vt=sqrt(vxgse^2+vygse^2+vzgse^2)
cs=sqrt(2*1.38e-23*temp/(den*1.673e-27))/1000.0

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

tiltangle=0.0*bxgsm
for i = 0L,rows-1 do begin
  iDay = ymd2dn(mat(2,i),mat(1,i),mat(0,i))
  recalc,mat(2,i),iDay,mat(3,i),mat(4,i),mat(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)', mat(1,0), mat(0,0), mat(2,0), $
                                 mat(3,0), mat(4,0), mat(5,0)
   print,'Last timestamp (MM DD YYYY  HH MM SS): '
   print, format='(i,i,i,i,i,f)', mat(1,rows-1), mat(0,rows-1), mat(2,rows-1), $
                                  mat(3,rows-1), mat(4,rows-1), mat(5,rows-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( mat(1,0), mat(0,0), mat(2,0), $
			  mat(3,0)+starthr, mat(4,0), mat(5,0))
   if (tmpStartTime lt startTime) then begin
     starthr = mat(3,0)
     print,'Warning:  You entered an invalid time.  Defaulting to', starthr
   endif

   tmpEndTime = julday( mat(1,0), mat(0,0), mat(2,0), $
		        mat(3,0)+starthr+endhr, mat(4,0), mat(5,0))

   if (tmpEndTime gt endTime) then begin     
     endhr = (mat(0,rows-1)-mat(0,0))*24.0 + $
	     (mat(3,rows-1)-mat(3,0)) + $
	     (mat(4,rows-1)-mat(4,0))/60.0 + $
	     (mat(5,rows-1)-mat(5,0))/(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
   finaltime_jd = julday( mat(1,0),mat(0,0),mat(2,0), $
			  mat(3,0),finaltime,mat(5,0) )-startTime

   finalvxgsm=spline(jd,vxgsm,finaltime_jd)
   finalvygsm=spline(jd,vygsm,finaltime_jd)
   finalvzgsm=spline(jd,vzgsm,finaltime_jd)
   finalvxsm=spline(jd,vxsm,finaltime_jd)
   finalvysm=spline(jd,vysm,finaltime_jd)
   finalvzsm=spline(jd,vzsm,finaltime_jd)
   finalden=spline(jd,den,finaltime_jd)
   finalcs=spline(jd,cs,finaltime_jd)
   finalbxgsm=spline(jd,bxgsm,finaltime_jd)
   finalbygsm=spline(jd,bygsm,finaltime_jd)
   finalbzgsm=spline(jd,bzgsm,finaltime_jd)
   finalbxsm=spline(jd,bxsm,finaltime_jd)
   finalbysm=spline(jd,bysm,finaltime_jd)
   finalbzsm=spline(jd,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,tiltangle,finaltime_jd)

   ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   ; now plot the SW Data in the standard form ;
   ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   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(mat(2,0),mat(1,0),mat(0,0))
    close,11
    openw,11,smfile
    printf,10,format='(i5,i4,2i3)',$
           mat(2,0),iDay,floor(StartHr),(starthr-floor(starthr))*60
    printf,11,format='(i5,i4,2i3)',$
           mat(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(mat(1,0), mat(0,0), mat(2,0), starthr,0,0)     
     actualStopTime = julday(mat(1,0),mat(0,0), mat(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