# 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, inipos, '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()