apps121025_test_liouville.app04_nrm_setup

A normalized setup. Adopted from apps121025_test_liouville.app02_setup

Normalization is done as:

\[\begin{split}X_n &= x / (GM)^{0.5} \\ V_n &= v / (GM)^{0.25} \\ T_n &= t / (GM)^{0.25}\end{split}\]

Then, you may get the simplest equation of motion.

\[\begin{split}dX_n/dT_n &= dV_n \\ dV_n/dT_n &= -\frac{X_n}{|X_n|^3}\end{split}\]

where

\[\begin{split}GM &= 6.674674\times 10^{-11} * 2\times 10^{27} = 1.3348\times 10^{17} \\ (GM)^{0.5} &= 3.6535\times 10^8 \\ (GM)^{0.25} &= 1.9114\times 10^4\end{split}\]

Initial position, originall x=4.5e8 becomes 1.2317. Initial velocity, originally v=10e3 becomes 0.52318. Time integration, originally t=3e4 becomes 1.56953.

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.

r''' A normalized setup. Adopted from :mod:`apps121025_test_liouville.app02_setup`

Normalization is done as:

.. math::

    X_n &= x / (GM)^{0.5}   \\
    V_n &= v / (GM)^{0.25}   \\
    T_n &= t / (GM)^{0.25}

Then, you may get the simplest equation of motion.

.. math::

    dX_n/dT_n &= dV_n  \\
    dV_n/dT_n &= -\frac{X_n}{|X_n|^3}

where

.. math::

    GM &= 6.674674\times 10^{-11} * 2\times 10^{27} = 1.3348\times 10^{17}  \\
    (GM)^{0.5} &= 3.6535\times 10^8  \\
    (GM)^{0.25} &= 1.9114\times 10^4

Initial position,  originall x=4.5e8 becomes 1.2317.
Initial velocity, originally v=10e3 becomes 0.52318.
Time integration, originally t=3e4 becomes 1.56953.

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
    '''
    kep = kepler.KeplerRungeKutta(1., 1., pos, vel, gravity_constant=1.)

    ### Trace the particle along the time.
    tlist = np.arange(0, 1.51, 0.01)

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

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

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

    plt.figure()
    plt.subplot(211)
    plt.plot(inipos[0], inipos[1], 'o')
    plt.plot(pos1[:, 0], pos1[:, 1], 'o')
    plt.subplot(212)
    plt.plot(vel1[:, 0], vel1[:, 1], 'o')
    plt.gca().set_aspect('equal')

    plt.savefig('app04_nrm_setup.png')

if __name__ == "__main__":
    main()