;Edward Burin des Roziers
;ASEN5050

;THEMIS.PRO is the main program for simulating the Lunar perturbations on the Themis probes.

pro themis

Re = 6378.14d

t = 0.0                                        ;Initial time sec
h = 60.0                                ;Time step sec

;Initial conditions of satellite
posvel_Sat = sat(t)

;set up numerical integration
tbody0 = three_body(t,posvel_Sat)

t_length = 43200
time = dblarr(t_length)
x = dblarr(t_length)
y = dblarr(t_length)
z = dblarr(t_length)
xdot = dblarr(t_length)
ydot = dblarr(t_length)
zdot = dblarr(t_length)
a_Sat = dblarr(t_length)
inc_Sat = dblarr(t_length)

for i=0l,t_length-1 do begin
        print, t_length - i
        posvel_Sat = rk4(posvel_Sat, tbody0, t, h, 'three_body', /double)
        time[i] = t
        aeiomwnu = c2kr(posvel_Sat[0:2],posvel_Sat[3:5])
        a_Sat[i] = aeiomwnu[0]
        inc_Sat[i] = aeiomwnu[2]
        x[i] = posvel_Sat[0]
        y[i] = posvel_Sat[1]
        z[i] = posvel_Sat[2]
        xdot[i] = posvel_Sat[3]
        ydot[i] = posvel_Sat[4]
        zdot[i] = posvel_Sat[5]
        t = t+h
        tbody0 = three_body(t,posvel_Sat)

endfor

time_min = time/60.
time_hrs = time_min/60.
time_dys = time_hrs/24.
inc_Sat = inc_Sat * !radeg                        ;deg

!P.multi = [0,1,3]
!P.charsize = 2

set_plot,'eps'
device, file='tbody.ps', xsize=26.0, ysize=20.1, xoffset=1.0, yoffset=1.0,$
        /landscape

plot, time_dys, x, title='Satellite X-Position',$
        xrange=[0,30], xstyle=1, ytitle='km',$
        xtitle='Days'
plot, time_dys, y, title='Satellite Y-Position',$
        xrange=[0,30], xstyle=1, ytitle='km',$
        xtitle='Days'
plot, time_dys, z, title='Satellite Z-Position',$
        xrange=[0,30], xstyle=1, ytitle='km',$
        xtitle='Days'

device, /close
set_plot, 'x'

end