PRO Dst_empirical, inputfile, angle_gsm_file

; Read inputfile (simulation output in non-dimensional form)

t_sim = fltarr(1)
n_sim = fltarr(1)
vx_sim = fltarr(1)
vy_sim = fltarr(1)
vz_sim = fltarr(1)
by_sim = fltarr(1)
bz_sim = fltarr(1)
;angle_gsm = fltarr(1)
v = fltarr(1)

restore,angle_gsm_file
angle_gsm=gamma_SIMtoGSM ; [deg]

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, inputfile, /get_lun

; Loop until EOF is found
while ~ eof(lun) do begin
   ; Read a line of text
   readf, lun, tsim, nsim, vxsim, vysim, vzsim, bysim, bzsim
   t_sim = [t_sim,tsim]
   n_sim = [n_sim,nsim]
   vx_sim = [vx_sim,vxsim]
   vy_sim = [vy_sim,vysim]
   vz_sim = [vz_sim,vzsim]
   v = [v,sqrt(vxsim^2.+vysim^2.+vzsim^2.)]
   by_sim = [by_sim,bysim]
   bz_sim = [bz_sim,bzsim]
endwhile

; Close the file
close, lun

t_sim = t_sim(1:*)
n_sim = n_sim(1:*)
vx_sim = vx_sim(1:*)
vy_sim = vy_sim(1:*)
vz_sim = vz_sim(1:*)
v = v(1:*)
by_sim = by_sim(1:*)
bz_sim = bz_sim(1:*)

; Transform Bz to GSM coordinates from simulation coordinates
angle_gsm = angle_gsm*!PI/180. ; rad

bz_gsm=by_sim*sin(angle_gsm)+bz_sim*cos(angle_gsm)

;pram=0.5*n_sim*v^2.
pram=((n_sim*rhodim)*(v*vdim)^2.)*1.e-9 ; [nPa]

; Dst empirical relation
f = 7.26 ; [nT/sqrt(nPa)]
g = 11.0 ; [nT]

; Calculate VBs at each solution time
vbs_gsm = fltarr(1)

for i = 0, n_elements(bz_gsm)-1 do begin
  if (bz_gsm[i] lt 0.) then begin
    vbs_gsm = [vbs_gsm,abs(v[i]*vdim*bz_gsm[i]*bdim)*1000.] ; VBs in [mV/m]
  endif else begin
    vbs_gsm = [vbs_gsm,0.]
  endelse
endfor

vbs_gsm=vbs_gsm(1:*)

print,'Calculating the Dst index (O Brien and McPherron 2000).... '
print,'via Runge-Kutta 4th order method'

dst = fltarr(1)

print,'Specify initial conditions [0,0]:'

xi = 0.0 ; initial Dst=0 (assumed to be in quiet conditions)
ti = (t_sim(0)*tdim)/3600. ; [hours]

print,'ti:', ti
print,'Dst(ti):', xi

; Getting Dst*(ti)
xi = xi - f * sqrt(pram(0)) + g
dst1 = xi + f * sqrt(pram(0)) - g
dst = [dst,dst1]

for i = 1, n_elements(t_sim)-1 do begin

   tf = (t_sim(i)*tdim)/3600. ; [hours]

   h = tf - ti
   P1 = vbs_gsm(i)

   ; Calling rk4 function
   xf = rk4_pro(ti, xi, P1, h)

   ; Calculating the dst
   dst1 = xf + f * sqrt(pram(i)) - g

   ; Save into arrays
   dst = [dst,dst1]

   ti = tf
   xi = xf

endfor

   dst = dst(1:*)

   print, min(dst)
   print, max(dst)

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

   save,filename='Dst_empirical.sav',t_sim,dst

END