apps120911_iotorus.app15_fromsc_enabtrace

''' From the S/C, backtrace the particle.  Vsc=0 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 pyana.juice import jspice as js
import pyana.util.keplernumeric as kepler

def main():
    '''Main script'''

    js.init()

    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)

    for it in range(48):   # 8 hour.

        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)

    ### 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 = -sc2io / np.sqrt((sc2io * sc2io).sum()) * 74. / 71500.   # in Rj/s
    print(enavel)
    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.')

    ### 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
    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_aspect('equal')
    axxz.set_aspect('equal')
    fig.savefig('app15_fromsc_enabtrace.png')

    return

if __name__ == "__main__":
    main()