''' From the S/C, backtrace the particle. Vsc from spice assumed.
As a model data, 2031-08-20 to 08-21 was used.
It is at the ganymede orbit with inclination of -23.5 deg.
17.5 Rj.
'''
import datetime
import numpy as np
import matplotlib.pyplot as plt
#from matplotlib import rc
#rc('text', usetex=True)
#rc('font', family='serif')
from pyana.juice import jspice as js
import pyana.util.keplernumeric as kepler
def main():
'''Main script'''
# Initialize spice.
js.init()
# Define parameters related to time.
dts = 600
dt = datetime.timedelta(seconds=dts)
t0 = datetime.datetime(2031, 8, 21, 15, 0, 0)
rj = 71500
fig = plt.figure(figsize=(14, 6))
axxy = fig.add_subplot(121)
axxz = fig.add_subplot(122)
### Plot of the s/c moon positions
for it in range(48): # 8 hour. = 10 min * 48
# Time to backtrace.
t = t0 - it * dt
print(t)
# First, Moons position and SC position calculated in JEqS frame.
io = js.get_position(t, frame='JEQS', target='IO') / rj
# eu = js.get_position(t, frame='JEQS', target='EUROPA') / rj
# ga = js.get_position(t, frame='JEQS', target='GANYMEDE') / rj
# ca = js.get_position(t, frame='JEQS', target='CALLISTO') / rj
# Then, s/c position is calculated in JEqS frame
sc = js.get_position(t, frame='JEQS') / rj
if t == t0:
l = 'o'
axxy.text(io[0], io[1], 'I')
# axxy.text(eu[0], eu[1], 'E')
# axxy.text(ga[0], ga[1], 'G')
# axxy.text(ca[0], ca[1], 'C')
axxz.text(io[0], io[2], 'I')
# axxz.text(eu[0], eu[2], 'E')
# axxz.text(ga[0], ga[2], 'G')
# axxz.text(ca[0], ca[2], 'C')
else:
l = '.'
axxy.plot(io[0], io[1], 'b' + l)
# axxy.plot(eu[0], eu[1], 'b' + l)
# axxy.plot(ga[0], ga[1], 'b' + l)
# axxy.plot(ca[0], ca[1], 'b' + l)
axxy.plot(sc[0], sc[1], 'r' + l)
axxz.plot(io[0], io[2], 'b' + l)
# axxz.plot(eu[0], eu[2], 'b' + l)
# axxz.plot(ga[0], ga[2], 'b' + l)
# axxz.plot(ca[0], ca[2], 'b' + l)
axxz.plot(sc[0], sc[2], 'r' + l)
# Plot of Jupiter body
axxy.plot(1.0 * np.cos(np.linspace(0, 2 * np.pi, 120)),
1.0 * np.sin(np.linspace(0, 2 * np.pi, 120)), 'k')
# Plot of Io orbit
axxy.plot(5.9 * np.cos(np.linspace(0, 2 * np.pi, 120)),
5.9 * np.sin(np.linspace(0, 2 * np.pi, 120)), 'k:')
### Looking vector with velocity of 74 km/s
scpos_t0 = js.get_position(t0, frame='JEQS') / rj
sc2io = js.get_position(t0, frame='JEQS', target='IO', origin='JGO')
enavel_insc = -sc2io / np.sqrt((sc2io * sc2io).sum()) * 74. # in km/s
scvel_t0 = js.get_velocity(t0, frame='JEQS')
enavel = enavel_insc - scvel_t0
print('VEL', enavel, enavel_insc, scvel_t0)
axxy.plot(scpos_t0[0] + [0, enavel[0] / 10.], scpos_t0[1] + [0, enavel[1] / 10.], 'k')
axxy.plot(scpos_t0[0] + [0, enavel_insc[0] / 10.], scpos_t0[1] + [0, enavel_insc[1] / 10.], 'k')
axxy.plot(scpos_t0[0] + [0, scvel_t0[0] / 10.], scpos_t0[1] + [0, scvel_t0[1] / 10.], 'k')
enavel = enavel / rj
enavel_insc = enavel_insc / rj
for it in range(48):
time = -it * dts
print(time, time*enavel)
enapos = scpos_t0 + time * enavel
axxy.plot(enapos[0], enapos[1], 'g.')
axxz.plot(enapos[0], enapos[2], 'g.')
enapos = scpos_t0 + time * enavel_insc
axxy.plot(enapos[0], enapos[1], 'g+')
axxz.plot(enapos[0], enapos[2], 'g+')
### Kepler motion. in MKSA
mj = 1.89813e27
scpos_t0 = js.get_position(t0, frame='JEQS') * 1e3 # in m
sc2io = js.get_position(t0, frame='JEQS', target='IO', origin='JGO') * 1e3 # in m
# enavel = -sc2io / np.sqrt((sc2io * sc2io).sum()) * 74e3 # in m/s
enavel = enavel * rj * 1e3
kep = kepler.Kepler2RungeKutta(mj, 1., scpos_t0, enavel)
tlist = np.arange(0, -600 * 48 - 1, -600)
print(tlist)
pos = kep.get_pos(tlist) / 1e3 / rj
axxy.plot(pos[:, 0], pos[:, 1], 'c.')
axxz.plot(pos[:, 0], pos[:, 2], 'c.')
axxy.set_xlim(-25, 10)
axxy.set_ylim(-10, 15)
axxy.set_aspect('equal')
axxz.set_aspect('equal')
fig.savefig('app16_fromscram_enabtrace.png')
fig.savefig('app16_fromscram_enabtrace.eps')
return
if __name__ == "__main__":
main()