PRO Dst_vs_NumMagStormIndex_fig8, inputfileTotalKE, dst_measured_file, dst_empirical_file 

; Read inputfileTotalKE (simulation output)
; Read dst_measured_file (measured hourly equatorial Dst (Final) index)
; Read dst_empirical_file (empirical Dst index calculated from the simulation output at L1 point from O'Brien and McPherron, 2000)

t_simTotalKE = fltarr(1)
totalke_sim = fltarr(1)

restore,dst_measured_file

rhodim=5.*1.66054e-27*1.e6 ; [kg/m3]
vdim=3.5e-9/sqrt(4.e-7*!PI*rhodim) ; [m/s]
bdim=3.5e-9 ; [T]
tdim=6.378e6/vdim ; [s]

openr, lun, inputfileTotalKE, /get_lun

; Loop until EOF is found
while ~ eof(lun) do begin
   ; Read a line of text
   readf, lun, tsimTotalKE, totalkesim
   t_simTotalKE = [t_simTotalKE,tsimTotalKE]
   totalke_sim = [totalke_sim,totalkesim]
endwhile

; Close the file
close, lun

t_simTotalKE = t_simTotalKE(1:*)
totalke_sim = totalke_sim(1:*)

t_simTotalKE = t_simTotalKE + (16.*60.+47.0208)/tdim

wdef,0,800,400

utplot,t*3600.,Dst_Final,background=255,color=0,yrange=[-100.,10.],title='Dst & Numerical Magnetic Storm Index vs Time (05 April 2010 Storm)',ytitle='Dst (nT), total(delB_2.5RE_10RE(t)) (nT)'

; Dessler-Parker relation
U = 4.*!dpi*8.e15*3.e-5/(3.*4.e-7*!dpi)
delB_TotalKE = -2.*totalke_sim*90.*3.e4/(3.*U) ; PI*(100*RE^2.-6.25*RE^2.)*(0.3RE) for which the height of the ring current slab is a variable, in this case 0.3RE.

oplot,t_simTotalKE*tdim,gauss_smooth(delB_TotalKE,2),color=0,linestyle=3

restore,dst_empirical_file
oplot,t_sim*tdim,dst,color=0,linestyle=2

legend,['Dst (measured)','Dst (empirical)','Numerical Magnetic Storm Index'],linestyle=[0,2,3],/top,/right,textcolor=0,color=0,box=0

write_png,'Dst_05042010_Storm_Overlap.png',tvrd(/true)

END