'''Particle motion tracer.
Trace particles in
1. Inertial frame
2. Io fixed frame
Neutral particle trajectories are traced under two frames.
'''
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])
def io_kepler():
''' A function to calculate single particle in inertia frame.
'''
kep_io = Kepler(mj, mio, iopos, iovel)
return kep_io
def singleparticle(vx_io, vy_io, vz_io):
''' A single particle. Parmeters are initial velocity in Io fixed frame.
'''
# Start from io position with velocity +1 km/s in y.
parti = Kepler(mj, mio, iopos, iovel + np.array([vx_io, vy_io, vz_io]))
return parti
def calinertia_dispinertia():
plt.subplot(221)
tlist = np.linspace(0, 36000, 30)
kep_io = io_kepler()
pos_io = kep_io.get_pos(tlist)
plt.plot(pos_io[:, 0], pos_io[:, 1], 'b:')
kep_part = singleparticle(-57e3, 0, 0)
pos_part = kep_part.get_pos(tlist)
plt.plot(pos_part[:, 0], pos_part[:, 1], 'r-o', label='-vx')
kep_part = singleparticle(57e3, 0, 0)
pos_part = kep_part.get_pos(tlist)
plt.plot(pos_part[:, 0], pos_part[:, 1], 'r-*', label='+vx')
kep_part = singleparticle(0, -57e3, 0)
pos_part = kep_part.get_pos(tlist)
plt.plot(pos_part[:, 0], pos_part[:, 1], 'r-x', label='-vy')
kep_part = singleparticle(0, 57e3, 0)
pos_part = kep_part.get_pos(tlist)
plt.plot(pos_part[:, 0], pos_part[:, 1], 'r-^', label='+vy')
plt.legend(loc='lower left')
plt.gca().set_aspect('equal')
plt.xlim(-.2e10, .2e10)
plt.ylim(-.2e10, .2e10)
plt.title('Calc inertia. Disp inertia.')
def inertia2io(t, vec, inverse=False):
''' Inerital to io frame, omega = wio
'''
theta = wio * t
if inverse:
theta = -theta
xax_io = np.array([np.cos(theta), np.sin(theta), 0])
yax_io = np.array([-np.sin(theta), np.cos(theta), 0])
vec_io = np.array([xax_io.dot(vec), yax_io.dot(vec), vec[2]])
return vec_io
def calinertia_dispioframe():
plt.subplot(222)
tlist = np.linspace(0, 36000, 30)
kep_io = io_kepler()
pos_io = kep_io.get_pos(tlist).T
pos_io_io = np.array([inertia2io(tlist[i], pos_io[:, i]) for i in range(len(tlist))])
plt.plot(pos_io_io[:, 0], pos_io_io[:, 1], 'b:')
kep_part = singleparticle(-57e3, 0, 0)
pos_part = kep_part.get_pos(tlist).T
pos_part_io = np.array([inertia2io(tlist[i], pos_part[:, i]) for i in range(len(tlist))])
plt.plot(pos_part_io[:, 0], pos_part_io[:, 1], 'r-o', label='-vx')
kep_part = singleparticle(57e3, 0, 0)
pos_part = kep_part.get_pos(tlist).T
pos_part_io = np.array([inertia2io(tlist[i], pos_part[:, i]) for i in range(len(tlist))])
plt.plot(pos_part_io[:, 0], pos_part_io[:, 1], 'r-*', label='+vx')
kep_part = singleparticle(0, -57e3, 0)
pos_part = kep_part.get_pos(tlist).T
pos_part_io = np.array([inertia2io(tlist[i], pos_part[:, i]) for i in range(len(tlist))])
plt.plot(pos_part_io[:, 0], pos_part_io[:, 1], 'r-x', label='-vy')
kep_part = singleparticle(0, 57e3, 0)
pos_part = kep_part.get_pos(tlist).T
pos_part_io = np.array([inertia2io(tlist[i], pos_part[:, i]) for i in range(len(tlist))])
plt.plot(pos_part_io[:, 0], pos_part_io[:, 1], 'r-^', label='+vy')
plt.legend(loc='lower left')
plt.gca().set_aspect('equal')
plt.xlim(-.2e10, .2e10)
plt.ylim(-.2e10, .2e10)
plt.title('Calc inertia. Disp ioframe(cnv).')
def calioframe_dispinertia():
''' Calculation in IO /display in inertia
'''
plt.subplot(223)
tlist = np.linspace(0, 36000, 30)
kep_io = Kepler(mj, mio, iopos, np.array([0, 0, 0]), omega=np.array([0, 0, wio]))
pos_io = kep_io.get_pos(tlist).T
pos_io = np.array([inertia2io(tlist[i], pos_io[:, i], inverse=True) for i in range(len(tlist))])
plt.plot(pos_io[:, 0], pos_io[:, 1], 'b:')
part_io = Kepler(mj, mio, iopos, np.array([-57e3, 0, 0]), omega=np.array([0, 0, wio]))
part_io = part_io.get_pos(tlist).T
part_io = np.array([inertia2io(tlist[i], part_io[:, i], inverse=True) for i in range(len(tlist))])
plt.plot(part_io[:, 0], part_io[:, 1], 'ro')
part_io = Kepler(mj, mio, iopos, np.array([57e3, 0, 0]), omega=np.array([0, 0, wio]))
part_io = part_io.get_pos(tlist).T
part_io = np.array([inertia2io(tlist[i], part_io[:, i], inverse=True) for i in range(len(tlist))])
plt.plot(part_io[:, 0], part_io[:, 1], 'r*')
part_io = Kepler(mj, mio, iopos, np.array([0, -57e3, 0]), omega=np.array([0, 0, wio]))
part_io = part_io.get_pos(tlist).T
part_io = np.array([inertia2io(tlist[i], part_io[:, i], inverse=True) for i in range(len(tlist))])
plt.plot(part_io[:, 0], part_io[:, 1], 'rx')
part_io = Kepler(mj, mio, iopos, np.array([0, 57e3, 0]), omega=np.array([0, 0, wio]))
part_io = part_io.get_pos(tlist).T
part_io = np.array([inertia2io(tlist[i], part_io[:, i], inverse=True) for i in range(len(tlist))])
plt.plot(part_io[:, 0], part_io[:, 1], 'r^')
plt.legend(loc='lower left')
plt.gca().set_aspect('equal')
plt.xlim(-.2e10, .2e10)
plt.ylim(-.2e10, .2e10)
plt.title('Calc ioframe. Disp inertia(cnv).')
def calioframe_dispioframe():
''' Calculation/display both done in IO frame
'''
plt.subplot(224)
tlist = np.linspace(0, 36000, 30)
kep_io = Kepler(mj, mio, iopos, np.array([0, 0, 0]), omega=np.array([0, 0, wio]))
pos_io = kep_io.get_pos(tlist)
plt.plot(pos_io[:, 0], pos_io[:, 1], 'b:')
part_io = Kepler(mj, mio, iopos, np.array([-57e3, 0, 0]), omega=np.array([0, 0, wio]))
part_io = part_io.get_pos(tlist)
plt.plot(part_io[:, 0], part_io[:, 1], 'ro')
part_io = Kepler(mj, mio, iopos, np.array([57e3, 0, 0]), omega=np.array([0, 0, wio]))
part_io = part_io.get_pos(tlist)
plt.plot(part_io[:, 0], part_io[:, 1], 'r*')
part_io = Kepler(mj, mio, iopos, np.array([0, -57e3, 0]), omega=np.array([0, 0, wio]))
part_io = part_io.get_pos(tlist)
plt.plot(part_io[:, 0], part_io[:, 1], 'rx')
part_io = Kepler(mj, mio, iopos, np.array([0, 57e3, 0]), omega=np.array([0, 0, wio]))
part_io = part_io.get_pos(tlist)
plt.plot(part_io[:, 0], part_io[:, 1], 'r^')
plt.legend(loc='lower left')
plt.gca().set_aspect('equal')
plt.xlim(-.2e10, .2e10)
plt.ylim(-.2e10, .2e10)
plt.title('Calc ioframe. Disp ioframe.')
def main():
plt.figure(figsize=(10,10))
calinertia_dispinertia()
calinertia_dispioframe()
calioframe_dispinertia()
calioframe_dispioframe()
plt.plot([15.9 * rj], [0], 'bx')
plt.savefig('app10_motions.png')
if __name__ == "__main__":
retval = main()