; 

f=dblarr(91, 793)

openr,1, 'whi_08080_08113'
readf,1,f
close,1

year = reform( f(0,*))
doy =reform(f(1,*))
hour=reform(f(2,*))
hour=hour/100.

time=year + (doy-1)/365. + hour/24./365.
day = doy + hour/24.

kp=reform(f(3,*))
ap=reform(f(4,*))


bx= reform(f(8,*))
by= reform(f(9,*))
bz= reform(f(10,*))

bmag=sqrt(bx*bx+by*by+bz*bz)

den = reform(f(11,*))
vel = reform(f(12,*))


hps = reform(f(20,*))
zeros=where(hps eq 0)

hpn =  reform(f(34,*))
zeron=where(hpn eq 0)

hp=hpn+hps
hp(zeros)=2*hpn(zeros)
hp(zeron)=2*hps(zeron)


;flux=f(5,*)
;ssn=f(7,*)

openr,1, 'solar_whi.dat'
ff=fltarr(5,28)
readf,1,ff
close,1
ssn=ff(3,*)
flux=ff(4,*)
days=106+findgen(28)

; add radbelt hourly data

h=' '
f=fltarr(6, 10321-7)
openr,1, '2007-2008GOES12_1hour.txt'
for i=0, 6 do readf,1, h
readf,1,f
close,1

ryear =reform(f(0,*))
rmonth =reform(f(1,*))
rday =reform(f(2,*))
rhour =reform(f(3,*))
rminute =reform(f(4,*))

radbelt = reform(f(5,*))
missing = where( radbelt gt 1e6)
radbelt(missing) =-9999

; this works only for leap !!!!!1

mday_leap = [0,31,60,91,121,152,182,213,244,274,305,335]   ;leap year

rdoy = mday_leap(rmonth-1) + rday + double(rhour)/24.

ri=where(ryear eq 2008 and rdoy ge 80 and rdoy le 107+4)


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; whi period = CR 2068 = Mar 20 to April 16 2008 = doy 80 - 107
; but we want to add on four days to include wind travel time
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

i=where(doy ge 80 and doy le 107+4)

colors,/standard


; THIS IS THE PAPER PLOT FOR RADBELT

set_plot, 'ps'
device, bits=16, /color, /encapsulated, xsize=18, ysize=12, $
/bold, /helvetica, filename='whi_vel_rad.eps'

!p.font=0
!p.multi=[0,1,1]
!p.charsize=1.2
!p.symsize=.1

plot_io, rdoy(ri), radbelt(ri), xstyle=1, ystyle=5, min =0., xtitle='Day of year 2008', xmargin=[8,8],yrange=[1,100000],/nodata
;,title='Whole Heliosphere Interval'
axis,yaxis=1,ystyle=1,ytitle='Radiation belt(electrons/cm^2/s/sr)',color=16
oplot,rdoy(ri),radbelt(ri),color=16,min=0.,thick=8

plot, day(i), vel(i), ytitle='Solar Wind Speed (km/sec)', xstyle=5,ystyle=9,xmargin=[8,8],yrange=[300,750],/noerase,min=1,color=4
;oplot,day(i),-1*vel(i)*bz(i)*0.1,color=16,thick=4,max=600

device,/close

;; ENDS HERE

; THIS IS THE PAPER PLOT FOR AURORA

set_plot, 'ps'
device, bits=16, /color, /encapsulated, xsize=18, ysize=12, $
/bold, /helvetica, filename='whi_velbz_hemp.eps'

!p.font=0
!p.multi=[0,1,1]
!p.charsize=1.2
!p.symsize=.1

plot, day(i), hp(i), xstyle=1, ystyle=5, min =1, xtitle='Day of year 2008', xmargin=[8,8],yrange=[0,200],/nodata,title='Whole Heliosphere Interval'
axis,yaxis=1,ystyle=1,ytitle='Auroral Electron Power (GW)',color=2
oplot,day(i),hp(i),color=2,min=1,thick=8

plot, day(i), vel(i), ytitle='Solar Wind Speed (km/sec)', xstyle=5,ystyle=9,xmargin=[8,8],yrange=[0,750],/noerase,min=1,thick=8,color=4
oplot,day(i),-1*vel(i)*bz(i)*0.1,thick=4,max=600
;oplot,day(i),vel(i),color=4,thick=8,min=1,ytitle='Solar Wind Speed (km/sec)'

device,/close

;; ENDS HERE

set_plot, 'ps'
device, bits=16, /color, /encapsulated, xsize=18, ysize=12, $
/bold, /helvetica, filename='whi_velhemp.eps'

;set_plot,'x'

!p.font=0
!p.multi=[0,1,1]
!p.charsize=1.2
!p.symsize=.1

plot, day(i), hp(i), xstyle=1, ystyle=5, min =1, xtitle='Day of year 2008', xmargin=[8,8],yrange=[-10,110],/nodata
axis,yaxis=1,ystyle=1,ytitle='Auroral Electron Power (orange), Bz (purple)'
oplot,day(i),-1*bz(i)*10.,color=13,min=0,psym=2
oplot,day(i),hp(i),color=2,min=1

plot, day(i), vel(i), ytitle='Solar Wind Speed (blue)', xstyle=5,ystyle=9,xmargin=[8,8],yrange=[300,750],/noerase,title='Whole Heliosphere Interval',min=1,/nodata
oplot,day(i),vel(i),color=4,min=1,thick=8

device,/close

set_plot, 'ps'
device, bits=16, /color, /encapsulated, xsize=18, ysize=12, $
/bold, /helvetica, filename='whi_velmag.eps'

;set_plot,'x'

!p.font=0
!p.multi=[0,1,1]
!p.charsize=1.2

plot, day(i), bmag(i), xstyle=1, ystyle=5, min =1, xtitle='Day of year 2008', xmargin=[8,8],color=13,yrange=[0,10]
axis,yaxis=1,ystyle=1,ytitle='Solar Wind Magnetic Field',color=13

plot, day(i), vel(i), ytitle='Solar Wind Speed (km/s)', xstyle=5,ystyle=9,xmargin=[8,8],yrange=[300,750],/noerase,title='Whole Heliosphere Interval',min=1,color=4

device,/close

set_plot, 'ps'
device, bits=16, /color, /encapsulated, xsize=18, ysize=12, $
/bold, /helvetica, filename='whi_velmagz.eps'

;set_plot,'x'

!p.font=0
!p.multi=[0,1,1]
!p.charsize=1.2

plot, day(i), bz(i), xstyle=1, ystyle=5, xtitle='Day of year 2008', xmargin=[8,8],yrange=[-10,10],max=0,psym=2,color=13
axis,yaxis=1,ystyle=1,ytitle='Solar Wind Magnetic Field (Bz)',color=13

plot, day(i), vel(i), ytitle='Solar Wind Speed (km/s)', xstyle=5,ystyle=9,xmargin=[8,8],yrange=[300,750],/noerase,title='Whole Heliosphere Interval',min=1,color=4

device,/close
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

f=dblarr(91, 793)

openr,1, 'wsm_96223_96256'
readf,1,f
close,1

year = reform(f(0,*))
doy = reform(f(1,*))
hour= reform(f(2,*))
hour=hour/100.

time=year + (doy-1)/365. + hour/24./365.
day = doy + hour/24.

kp = reform(f(3,*))
ap = reform(f(4,*))

bx = reform(f(8,*))
by = reform(f(9,*))
bz = reform(f(10,*))

bmag=sqrt(bx*bx+by*by+bz*bz)

den = reform(f(11,*))
vel = reform(f(12,*))


hps = reform(f(20,*))
zeros=where(hps eq 0)

hpn =  reform(f(34,*))
zeron=where(hpn eq 0)

hp=hpn+hps
hp(zeros)=2*hpn(zeros)
hp(zeron)=2*hps(zeron)

; not sure if this is a  bad data point 
; hp(459)=0

openr,1, 'solar_wsm.dat'
ff=fltarr(5,30)
readf,1,ff
close,1
ssn=ff(3,*)
flux=ff(4,*)
days=223+findgen(30)



; add radbelt hourly data

h=' '
f=fltarr(6, 10231-7)
openr,1, '1996-1997GOES8_1hour.txt'
for i=0, 6 do readf,1, h
readf,1,f
close,1

ryear =reform(f(0,*))
rmonth =reform(f(1,*))
rday =reform(f(2,*))
rhour =reform(f(3,*))
rminute =reform(f(4,*))

radbelt = reform(f(5,*))
missing = where( radbelt gt 1e6)
radbelt(missing) =-9999

; this works only for leap !!!!!1

mday_leap = [0,31,60,91,121,152,182,213,244,274,305,335]   ;leap year

rdoy = mday_leap(rmonth-1) + rday + double(rhour)/24.

ri=where(ryear eq 1996 and rdoy ge 223 and rdoy le 250+4)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; wsm = CR 1912-1913 =  Aug 10 to Sep 8, 1996 = doy 223 - 252
; NOTE this is not a 27-day interval but a 30-day interval !!!!
; lets do 32 days starting from WSM start day of Aug 10
; so 223 -- 254

; CR 1912 ends on Aug 22 = doy 235 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


i=where(doy ge 223 and doy le 250+4)



; THIS IS THE PAPER PLOT FOR RADBELT

set_plot, 'ps'
device, bits=16, /color, /encapsulated, xsize=18, ysize=12, $
/bold, /helvetica, filename='wsm_vel_rad.eps'

!p.font=0
!p.multi=[0,1,1]
!p.charsize=1.2
!p.symsize=.1

plot_io, rdoy(ri), radbelt(ri), xstyle=1, ystyle=5, min =0., xtitle='Day of year 1996', xmargin=[8,8],yrange=[1,100000],/nodata
;,title='Whole Sun Month'
axis,yaxis=1,ystyle=1,ytitle='Radiation belt(electrons/cm^2/s/sr)',color=16
oplot,rdoy(ri),radbelt(ri),color=16,min=0.,thick=8

plot, day(i), vel(i), ytitle='Solar Wind Speed (km/sec)', xstyle=5,ystyle=9,xmargin=[8,8],yrange=[300,750],/noerase,min=1,color=4
;oplot,day(i),-1*vel(i)*bz(i)*0.1,color=16,thick=4,max=600

device,/close

; THIS IS THE END

; THIS IS THE PAPER PLOT FOR AURORA
set_plot, 'ps'
device, bits=16, /color, /encapsulated, xsize=18, ysize=12, $
/bold, /helvetica, filename='wsm_velbz_hemp.eps'

!p.font=0
!p.multi=[0,1,1]
!p.charsize=1.2
!p.symsize=.1

plot, day(i), hp(i), xstyle=1, ystyle=5, min =1, xtitle='Day of year 1996', xmargin=[8,8],yrange=[0,200],/nodata,title='Whole Sun Month'
axis,yaxis=1,ystyle=1,ytitle='Auroral Electron Power (GW)',color=2
oplot,day(i),hp(i),color=2,min=1,thick=8

plot, day(i), vel(i), ytitle='Solar Wind Speed (km/sec)', xstyle=5,ystyle=9,xmargin=[8,8],yrange=[0,750],/noerase,min=1,thick=8,color=4
oplot,day(i),-1*vel(i)*bz(i)*0.1,thick=4,max=600
;oplot,day(i),vel(i),color=4,thick=8,min=1,ytitle='Solar Wind Speed (km/sec)'

device,/close

; THIS IS THE END

device, bits=16, /color, /encapsulated, xsize=18, ysize=12, $
/bold, /helvetica, filename='wsm_velhemp.eps'

plot, day(i), hp(i), xstyle=1, ystyle=5, xtitle='Day of year 1996', xmargin=[8,8],yrange=[-10,110],min=1,/nodata
axis,yaxis=1,ystyle=1,ytitle='Auroral Electron Power (orange), Bz (purple)'
oplot,day(i),-1*bz(i)*10.,color=13,min=0,psym=2
oplot,day(i),hp(i),color=2,min=1

plot, day(i), vel(i), ytitle='Solar Wind Speed (blue)', xstyle=5,ystyle=9,xmargin=[8,8],yrange=[300,750],/noerase,title='Whole Sun Month',min=1,/nodata
oplot,day(i),vel(i),color=4,min=1,thick=8

device,/close

device, bits=16, /color, /encapsulated, xsize=18, ysize=12, $
/bold, /helvetica, filename='wsm_velmag.eps'

plot, day(i), bmag(i), xstyle=1, ystyle=5, xtitle='Day of year 1996', xmargin=[8,8],color=13,yrange=[0,10],min=1
axis,yaxis=1,ystyle=1,ytitle='Solar Wind Magnetic Field',color=13

plot, day(i), vel(i), ytitle='Solar Wind Speed (km/s)', xstyle=5,ystyle=9,xmargin=[8,8],yrange=[300,750],/noerase,title='Whole Sun Month',min=1,color=4


device,/close

device, bits=16, /color, /encapsulated, xsize=18, ysize=12, $
/bold, /helvetica, filename='wsm_velmagz.eps'

plot, day(i), bz(i), xstyle=1, ystyle=5, xtitle='Day of year 1996', xmargin=[8,8],color=13,yrange=[-10,10]
axis,yaxis=1,ystyle=1,ytitle='Solar Wind Magnetic Field (Bz)',color=13

plot, day(i), vel(i), ytitle='Solar Wind Speed (km/s)', xstyle=5,ystyle=9,xmargin=[8,8],yrange=[300,750],/noerase,title='Whole Sun Month',min=1,color=4

device,/close

stop
end