Numerical Orbit Propagation - Part 2

Continuing from the previous tutorial on Earth-bound numerical orbit propagation, this tutorial introduces two concepts:

1. Propagating an orbit around the Moon

2. Propagating an orbit around the Saturn

While propagation itself does not differ from the Earth-bound case, each of them is slightly different in the coordinate frames used in the propagation.

Propagating an Orbit around the Moon

For a propagation around the Moon, the procedure is very similar to the Earth-bound case. The first step is to initialise the orbit. Note that the orbit is initialised in the Lunar inertial coordinate frame MoonCRS. This frame is provided by SatMAD.

from astropy import units as u
from astropy.coordinates import CartesianDifferential, CartesianRepresentation
from astropy.time import Time
from satmad.coordinates.frames import init_pvt, MoonCRS, CelestialBodyCRS
from satmad.core.celestial_bodies_lib import MOON
from satmad.core.celestial_body import CelestialBody

from satmad.propagation.numerical_propagators import NumericalPropagator
from satmad.utils.timeinterval import TimeInterval

# Initialises a Moon low-altitude satellite.
time_moon = Time("2020-01-01T11:00:00.000", scale="utc")

v_moon_crs = CartesianDifferential([1, -1, 0.6], unit=u.km / u.s)
r_moon_crs = CartesianRepresentation([1000, 1000, 2000], unit=u.km)
rv_moon_crs = init_pvt(MoonCRS, time_moon, r_moon_crs, v_moon_crs)

print(rv_moon_crs)

<SkyCoord (MoonCRS: obstime=2020-01-01T11:00:00.000): (x, y, z) in km
    (1000., 1000., 2000.)
 (v_x, v_y, v_z) in km / s
    (1., -1., 0.6)>

The second step is to initialise the numerical propagator around Moon. Note that, the central body is given explicitly as MOON. This is a “CelestialBody” object and a number of default Celestial Bodies are already provided (such as EARTH and MOON) in the “celestial_bodies” module.

# propagation config params
output_stepsize = 60 * u.s

# init propagator
prop_moon = NumericalPropagator(
    output_stepsize,
    rtol=1e-12,
    atol=1e-14,
    central_body=MOON
)

The final step is to propagate the orbit. This is not any different from propagation around the Earth.

#propagation run params
prop_duration = 3 * u.day
prop_start = rv_moon_crs.obstime + 0.5 * u.day

# run propagation
trajectory_moon = prop_moon.gen_trajectory(rv_moon_crs, TimeInterval(prop_start, prop_duration))

print(trajectory_moon)

print(f"Coords at time {prop_start + 0.12345 * u.day}:\n{trajectory_moon(prop_start + 0.12345 * u.day)}")

Trajectory from 2020-01-01T23:00:00.000 to 2020-01-04T23:00:00.000 in frame mooncrs. (Interpolators initialised: False)
Coords at time 2020-01-02T01:57:46.080:
<SkyCoord (MoonCRS: obstime=2020-01-02T01:57:46.080): (x, y, z) in km
    (-1176.05184748, -996.4101513, -2226.35450763)
 (v_x, v_y, v_z) in km / s
    (-0.48262601, 1.29170032, 0.27677641)>

Propagating an Orbit around the Saturn

The last example is doing the same procedure for Saturn. However, Saturn is not one of the default objects made available by SatMAD. So we will have to create the SATURN object and the SaturnCRS, the Saturn Celestial Reference Frame.

class SaturnCRS(CelestialBodyCRS):
    """Saturn Celestial Reference System. This is simply the ICRS shifted to the
    centre of the Saturn."""

    body = "Saturn"
    ephemeris_type = "builtin"

SATURN = CelestialBody(
    "Saturn",
    "User Generated Saturn Model.",
    3.7931187E16 * u.m ** 3 / u.s ** 2,
    inert_coord=SaturnCRS,
)

It is important to note that, the name “Saturn” defined in the Celestial Reference Frame (CRS) definition can only be a Solar System Body, given in the list solar_system_ephemeris.bodies. This ensures that the planet is recognised, its coordinates are computed and coordinate transformations from another (e.g. Heliocentric) coordinate system can be properly computed. Also note that the ephemeris_type is set to builtin. This uses the builtin models found in Astropy, which are accurate enough for most applications. It is possible to set this to jpl for the default internal implementation of the current JPL model (e.g. DE430), but this may download a very large file the first time it is called.

Now that coordinate systems are initialised, we can generate the initial conditions using the SaturnCRS coordinate system.

# Initialises a Saturn-based satellite.
time_saturn = Time("2018-01-01T11:00:00.000", scale="utc")

v_saturn_crs = CartesianDifferential([7, 26, 0.5], unit=u.km / u.s)
r_saturn_crs = CartesianRepresentation([2000, 10000, 70000], unit=u.km)
rv_saturn_crs = init_pvt(SaturnCRS, time_saturn, r_saturn_crs, v_saturn_crs)

print(rv_saturn_crs)

<SkyCoord (SaturnCRS: obstime=2018-01-01T11:00:00.000): (x, y, z) in km
    (2000., 10000., 70000.)
 (v_x, v_y, v_z) in km / s
    (7., 26., 0.5)>

The next step is to initialise the numerical propagator around Saturn. Note that, the central body is given explicitly as SATURN, which is the CelestialBody object defined above. This time we initialise with higher tolerance values, but feel free to experiment.

# propagation config params
output_stepsize = 120 * u.s

# init propagator
prop_saturn = NumericalPropagator(
    output_stepsize,
    rtol=1e-10,
    atol=1e-11,
    central_body=SATURN
)

The final step is to propagate the orbit. Once again, this is not any different from propagation around the Saturn.

#propagation run params
prop_duration = 1 * u.day
prop_start = rv_saturn_crs.obstime + 0.5 * u.day

# run propagation
trajectory_saturn = prop_saturn.gen_trajectory(rv_saturn_crs, TimeInterval(prop_start, prop_duration))

print(trajectory_saturn)

print(f"Coords at time {prop_start + 0.45321 * u.day}:\n{trajectory_saturn(prop_start + 0.45321 * u.day)}")

Trajectory from 2018-01-01T23:00:00.000 to 2018-01-02T23:00:00.000 in frame saturncrs. (Interpolators initialised: False)
Coords at time 2018-01-02T09:52:37.344:
<SkyCoord (SaturnCRS: obstime=2018-01-02T09:52:37.344): (x, y, z) in km
    (16560.67147392, 63478.96515132, 54644.17965665)
 (v_x, v_y, v_z) in km / s
    (5.2815623, 19.15792542, -12.10055772)>