apps120911_iotorus.app11_ncparticle_track

'''Neutral cloud particle motion tracer.

Trace particles in
1. Inertial frame
2. Io fixed frame
'''

import os
import sys
import logging
logging.basicConfig()
import datetime
import math

import matplotlib.pyplot as plt
import numpy as np
import scipy as sp

from pyana.util.keplernumeric import Kepler2RungeKutta as Kepler

rj = 71500e3   # 71500 km.
rio = 5.9 * rj
vio = 17e3     # 17 km/s, roughly.
wio = vio / rio
mj = 1.8986e27   # Mass of Jupiter.
mio = 8.94e22    # Mass of Io
dt = 600.  # Time step for system
iopos = np.array([rio, 0., 0.])
iovel = np.array([0, vio, 0])

vcr = 74e3

def main():
    plt.figure(figsize=(14, 7))
    plt.subplot(121)

    for t in np.arange(0, 1.5e5, 1.08e4):

        phi = wio * t
        print(phi, np.rad2deg(phi), 'deg', t)

        # Position of Io (inertia)
        posio = np.array([np.cos(phi), np.sin(phi), 0]) * rio

        plt.plot([posio[0]], [posio[1]], 'ro')

        # Corotational velocity particle at Io (inertia)
        velcr = np.array([-np.sin(phi), np.cos(phi), 0]) * vcr
        
        plt.plot([posio[0], posio[0] + velcr[0] * 1e4],
            [posio[1], posio[1] + velcr[1] * 1e4], 'b:')

        # Trajectory
        trjctry = Kepler(mj, 1.0, posio, velcr)
        tlist = np.arange(0., 10800, 10)
        trjctry3d = trjctry.get_pos(tlist)

        plt.plot(trjctry3d[:, 0], trjctry3d[:, 1], 'b-')

    plt.gca().set_aspect('equal')

    plt.subplot(122)

    for vrel_km in np.arange(-30, 31, 10):
        vrel = vrel_km * 1e3

        phi = 0

        # Position of Io (inertia)
        posio = np.array([np.cos(phi), np.sin(phi), 0]) * rio

        plt.plot([posio[0]], [posio[1]], 'ro')

        # Corotational velocity particle at Io (inertia)
        velcr = np.array([-np.sin(phi), np.cos(phi), 0]) * (vcr + vrel)
        
        # Trajectory
        trjctry = Kepler(mj, 1.0, posio, velcr)
        tlist = np.arange(0., 10800 * 2.2, 10)
        trjctry3d = trjctry.get_pos(tlist)

        plt.plot(trjctry3d[:, 0], trjctry3d[:, 1])
        txt = '%.1f eV/amu' % (((vcr + vrel) / 13.8e3) ** 2)
        plt.text(trjctry3d[-1, 0], trjctry3d[-1, 1], txt, horizontalalignment='right')

    plt.plot(np.cos(np.linspace(0, 2 * np.pi, 60)) * 5.9 * rj,
                np.sin(np.linspace(0, 2 * np.pi, 60)) * 5.9 * rj, 'b:')

    plt.plot(np.cos(np.linspace(0, 2 * np.pi, 60)) * 15.9 * rj,
                np.sin(np.linspace(0, 2 * np.pi, 60)) * 15.9 * rj, 'b:')


    plt.xlim(-1e9, 1e9)
    plt.ylim(-0e9, 35 * rj)
    plt.gca().set_aspect('equal')
    plt.savefig('app11_ncparticle_track.png')


if __name__ == "__main__":
    retval = main()