irfpy.util.keplernumeric

Kepler’s trajectory using numeric solvers

Kepler’s trajectory is determined by parameters together with the initial conditions.

A generic form of the equation of motion is

\[\frac{d\mathbf{v}}{dt} = -GM\frac{\hat{\mathbf{r}}}{r^2}\]

where G is the constant of gravity and M is the mass of the center body (or center of mass).

Summary reports the details about two-body problems

Kepler2RungeKutta is the most reliable one for non-reltivistic Kepler motion with centrifugal force supported.

\[\begin{split}\frac{dx}{dt} = v_x \\ \frac{dy}{dt} = v_y \\ \frac{dz}{dt} = v_z \\ \frac{dv_x}{dt} = -Gm_0 x / r^3 + \omega^2 x + 2v_y\omega_z - 2 v_z\omega_y\\ \frac{dv_y}{dt} = -Gm_0 y / r^3 + \omega^2 y + 2v_z\omega_x - 2 v_x\omega_z\\ \frac{dv_z}{dt} = -Gm_0 z / r^3 + \omega^2 z + 2v_x\omega_y - 2 v_y\omega_x\end{split}\]

here r is \(\sqrt{x^2+y^2+z^2}\).

Kepler2RungeKutta(mass0, mass1, xinit, vinit)	Kepler motion seen in a rotation frame.
Kepler(mass0, mass1, xinit, vinit, dt[, ...])	Simulate Kepler motion using Leap Flog method.
Kepler_SunEarth([dt])	Return Kepler class for Sun-Earth system
Summary(mass0, mass1, xinit, vinit[, ...])	Report summary of the two body problem.

class irfpy.util.keplernumeric.Kepler(mass0, mass1, xinit, vinit, dt, gravity_constant=6.674e-11)

Bases: object

Simulate Kepler motion using Leap Flog method.

Initalizer

Parameters:

• mass0 – Mass of the central body (or equivalent mass in the center of mass frame) MKSA system is assumed, unless you specify gravity_constant.

• mass1 – Mass of the body of interest. It does not affect to the motion, but kept it by historical reason :(

• xinit (np.array with (3,)) – Initial position of mass1 at t=0

• vinit (np.array with (3,)) – Initial velocity of mass1 at t=0

• dt – Time resolution. Negative dt will be backtracing.

• gravity_constant – Gravity constant. The MKSA value for default.

USAGE:

Let’s see the Io motion around Jupiter. The mass of Jupiter is 1.89686e27 kg, and Io is 8.932e22 (which does not affect the calculation). Put the Io at 4.217e5 km with velocity of 17.334 km/s at t=0. Time resolution is set 3600 s (1 hour).

>>> kepler = Kepler(1.89686e27, 8.932e22,
...                 np.array((4.217e8, 0, 0)),
...                 np.array((0, 1.7334e4, 0)),
...                 3600)

Then, you can proceed the calculation by 100 steps. This means that the calculation was done until 3.6e5 s (100 hours).

>>> kepler.proceed(100)    # 100 step forward.

You may find the position at t=7.6e4 s, which is about the half of the orbital period of Io (full orbital period is 1.52853e5 s). Thus, the position should be close to [-4.217e8, 0, 0] with velocity of (0, -1.7334e4, 0).

>>> pos = kepler.get_pos(7.6e4)
>>> print(pos)      
[-4.26590115e+08   1.92925628e+07   0.00000000e+00]

>>> vel = kepler.get_vel(7.6e4)
>>> print(vel)    
[  -778.52664693 -17145.37459877      0.        ]

proceed(niter)

Proceed n steps.

step()

Proceed 1 step.

get_xlist()

Return the list of positions.

Returns:

Tuple. (xtlist, xlist), where xtlist is (N,) shape and xlist is (N, 3) shape.

get_vlist()

Return the list of velocity.

Returns:

Tuple. (vtlist, vlist), where vtlist is (N,) shape and vlist is (N, 3) shape.

get_pos(t)

Return the position, interpolated.

Parameters:

t – Scalar or tuple of the time. Should be interpolated.

Returns:

np.array with (3,) or (3, N) shape.

>>> kep = Kepler_SunEarth(dt=86400)
>>> kep.proceed(100)  # 100 day move.
>>> tx, xyz = kep.get_xlist()   # Return results in list format.
>>> print(tx[38])
3283200
>>> x, y, z = xyz[38, :]  # Position of 38th data.
>>> print('%.5e %.5e' % (x, y))
1.18759e+11 9.09623e+10
>>> x, y, z = kep.get_pos(3300000)
>>> print('%.5e %.5e' % (x, y))
1.18454e+11 9.13590e+10
>>> kep.proceed(100)  # 100 day more.
>>> x, y, z = kep.get_pos(86400 * 180.38)  # At 180.38 days
>>> print('%.5e %.5e' % (x, y))
-1.49434e+11 5.65757e+09

get_vel(t)

Return the velocity.

irfpy.util.keplernumeric.Kepler_SunEarth(dt=3600)

Return Kepler class for Sun-Earth system

>>> kse = Kepler_SunEarth()
>>> kse.proceed(24 * 365)  # About 1 year.
>>> tx, x = kse.get_xlist()
>>> print('%.2e %.2e %.2e' % (x[-1, 0], x[-1, 1], x[-1, 2]))  # Last data
1.50e+11 -3.11e+08 0.00e+00

irfpy.util.keplernumeric.Kepler_EarthMoon()

Return Kepler class for Earth-moon system

class irfpy.util.keplernumeric.KeplerRungeKutta(mass0, mass1, xinit, vinit, gravity_constant=6.674e-11)

Bases: object

Kepler motion solver using Runge-Kutta.

Similar functionality to Kepler, but using classical RK solver. No leap flog used.

Equation of motion is

\[\begin{split}\frac{dx}{dt} = v_x \\ \frac{dy}{dt} = v_y \\ \frac{dz}{dt} = v_z \\ \frac{dv_x}{dt} = -Gm_0 x / r^3 \\ \frac{dv_y}{dt} = -Gm_0 y / r^3 \\ \frac{dv_z}{dt} = -Gm_0 z / r^3\end{split}\]

here r is \(\sqrt{x^2+y^2+z^2}\).

>>> rkkep = KeplerRungeKutta(1e27, 1.5, (5e8, 0, 0), (0, 20e3, 0))
>>> x, y, z, vx, vy, vz = rkkep.get_posvel(np.arange(0, 1000, 30)).T
>>> print('%.4e %.4e %.4e' % (x[-1], y[-1], z[-1]))
4.9987e+08 1.9798e+07 0.0000e+00

Compare with Kepler class. >>> kep = Kepler(1e27, 1.5, np.array((5e8, 0, 0)), … np.array((0, 20e3, 0)), 10) >>> kep.proceed(150) >>> x, y, z = kep.get_pos(990) >>> print(‘%.4e %.4e %.4e’ % (x, y, z)) 4.9987e+08 1.9798e+07 0.0000e+00

Parameters:

• mass0 – The center object.

• mass1 – The flying object. It does not affect the results.

• xinit – (3,) shaped numpy.array to specify the initial position of mass1.

• vinit – (3,) shpaed numpy.array to specify the initial velocity of mass1.

• gravity_constant – The gravity constant. You may change it if the equation and variables are normalized.

get_posvel(tlist)

Return the pos, vel

get_pos(tlist)

get_vel(tlist)

class irfpy.util.keplernumeric.Kepler2RungeKutta(mass0, mass1, xinit, vinit, omega=array([0., 0., 0.]), gravity_constant=6.674e-11)

Bases: object

Kepler motion seen in a rotation frame.

Kepler motion in a rotation frame. The center of the rotation MUST be at the origin, where massive object (mass0) is located. Note that the \(\omega=0\) case, the results should return exactly the same as KeplerRungeKutta.

Probably it can handle the “center of gravity” frame, but not yet checked…

Equation of motion is:

\[\begin{split}\frac{dx}{dt} = v_x \\ \frac{dy}{dt} = v_y \\ \frac{dz}{dt} = v_z \\ \frac{dv_x}{dt} = -Gm_0 x / r^3 + \omega^2 x + 2v_y\omega_z - 2 v_z\omega_y\\ \frac{dv_y}{dt} = -Gm_0 y / r^3 + \omega^2 y + 2v_z\omega_x - 2 v_x\omega_z\\ \frac{dv_z}{dt} = -Gm_0 z / r^3 + \omega^2 z + 2v_x\omega_y - 2 v_y\omega_x\end{split}\]

here r is \(\sqrt{x^2+y^2+z^2}\).

>>> rkkep = Kepler2RungeKutta(1e27, 1.5, (5e8, 0, 0), (0, 20e3, 0),
... omega=np.array([0.0, 0.0, 0.0]))

The simplest way of obtaining the position and velocity is to give a list of the time to integrate. Usually, it starts from 0 and increases monotonically with a specific time step. Note that the integral starts with the first element of the time step.

>>> x, y, z, vx, vy, vz = rkkep.get_posvel(np.arange(0, 30001, 30)).T
>>> print('%.4e %.4e %.4e' % (x[0], y[0], z[0]))
5.0000e+08 0.0000e+00 0.0000e+00
>>> print('%.4e %.4e %.4e' % (x[-1], y[-1], z[-1]))  # Data on t=30000.
4.0260e+08 5.6589e+08 0.0000e+00

If you switch on the rotation, you can instance as

>>> rkkep_rott = Kepler2RungeKutta(1e27, 1.5, (5e8, 0, 0), (0, 15e3, 0),
... omega=np.array([0, 0, 5e3 / 5e8]))

Here the rotation frame is rotating with angular frequency of 5e3 / 5e8 = 1e-5 rad/s. This means that the initial velocity of the particle, measured in rotation frame, is 5e3 less than that seen in inertia. This corresponds to 15e3.

>>> x,  y, z, vx, vy, vz = rkkep_rott.get_posvel(np.arange(0, 30001, 30)).T
>>> print('%.4e %.4e %.4e' % (x[0], y[0], z[0]))
5.0000e+08 0.0000e+00 0.0000e+00
>>> print('%.4e %.4e %.4e' % (x[-1], y[-1], z[-1]))  # Data on t=30000.
5.5185e+08 4.2164e+08 0.0000e+00

The final coordinates at t=30000 in rotation frame can be (have been) validated from the rotational conversion using \(\omega t=0.3\) rad.

Parameters:

• mass0 – The center object.

• mass1 – The flying object. It does not affect the results.

• xinit – (3,) shaped numpy.array to specify the initial position of mass1. Given in the rotation frame.

• vinit – (3,) shpaed numpy.array to specify the initial velocity of mass1. Given in the rotation frame.

• omega – (3,) shpaed the angular velocity vector of the rotation frame seen from inertia frame.

• gravity_constant – The gravity constant. You may change it if the equation and variables are normalized.

get_posvel(tlist)

Return the pos, vel

Get the posistion and velocity as (N, 3) shape where N is the shape of tlist.

Note that the tlist should start from t=0. If tlist[0] is not 0, the initial conditions are considered at the time of tlist[0]. Thus, the results will be “shifted”.

get_pos(tlist)

get_vel(tlist)

class irfpy.util.keplernumeric.Summary(mass0, mass1, xinit, vinit, gravity_constant=6.674e-11)

Bases: object

Report summary of the two body problem.

See http://www.mns.kyutech.ac.jp/~kishine/notes/sentan_note101130e.pdf for summary of conservables.

>>> msun = 1.9884e30  # kg
>>> mearth = 5.972e24 # kg
>>> rearth = 149597870700.   #m
>>> vearth = 29780.  # m/s

>>> summary = Summary(msun, mearth, [rearth, 0, 0], [0, vearth, 0])
>>> # summary.print_summary() # will show the summary

print_summary()

Print summary of the two-body system.

get_eccentricity_vector()

Return the eccentricity vector.

Use get_eccentricity() to get scalar data.

get_eccentricity()

Return the eccentricity.

If eccentricity, e, is 0, the motion is circular. For e between 0 and 1, the motion is eliptic. For e=1, the motion is parabolic. For e>1, the motion is hyperbolic.

get_total_energy()

Return the total energy

Usually the unit is in Joule, if you instance this object in MKSA system.

get_angular_momentum()

Return the angular momentum.

If you specify the initial values in MKSA system, the angular momentum vector is in the unit of \(kg m^2/s\).

get_areal_velocity()

Return the areal velcoity vector

get_laplace_vector()

Return the Laplace vector

irfpy.util.keplernumeric.doctests()