apps121025_test_liouville.app02_setupΒΆ

A simple setup. Refactoring of app01

Set a massive star in the center. Set a particle (reference) at the arbitrary position. Set a particle velocity (reference) at

To eliminate the problem into 2D x 2D for the first step, the particle position is on x-axis and the velocity has only vx, vy components.

''' A simple setup. Refactoring of app01

Set a massive star in the center.
Set a particle (reference) at the arbitrary position.
Set a particle velocity (reference) at 

To eliminate the problem into 2D x 2D for the first step,
the particle position is on x-axis and the velocity
has only vx, vy components.
'''

import numpy as np
import matplotlib.pyplot as plt

import irfpy.util.keplernumeric as kepler

def kep_particle(pos, vel):
    ''' Kepler's particle simulation for a specific particle
    '''
    ### Set the mass (center).  Assume Jupiter (roughly).
    ms = 2e27 # kg

    ### Instance the kepler class
#    kep = kepler.Kepler2RungeKutta(ms, 1., pos, vel)
    kep = kepler.KeplerRungeKutta(ms, 1., pos, vel)

    ### Trace the particle along the time.
    tlist = [0, 3e4]

    posvel = kep.get_posvel(tlist)
    return posvel[-1, :3], posvel[-1, 3:]

def main():
    ### Set the particle position. Assume Io position.
    inipos = np.array([4.5e8, 0, 0])  # m

    ### Set the particle velocity. Assume 10 km/s at 60 deg inward from x.
    inivel = np.array([-10e3 * 0.5, 10e3 * np.sqrt(3./4), 0])
    
    pos1, vel1 = kep_particle(inipos, inivel)

    plt.figure()
    plt.plot(inipos[0], inipos[1], 'o')
    plt.plot(pos1[0], pos1[1], 'o')
    plt.gca().set_aspect('equal')

    plt.savefig('app02_setup.png')

if __name__ == "__main__":
    main()