Astronomical Computations

The Tracker Component Library provides astronomical and orbital mechanics functions, including JPL ephemeris queries, relativistic corrections, and orbital dynamics.

Note

The ephemeris functions require the astronomy extra (pip install nrl-tracker[astronomy]), which provides jplephem and astropy. The relativity, orbital mechanics, reference frame, time system, and Lambert functions need only the core install.

JPL Development Ephemeris

The ephemeris module provides access to JPL’s high-precision Development Ephemeris (DE) files for computing accurate positions and velocities of celestial bodies. Positions are returned in AU and velocities in AU/day, relative to the Solar System Barycenter (SSB) unless stated otherwise.

Basic Usage

from pytcl.astronomical.ephemerides import DEEphemeris
import numpy as np

# Create ephemeris object (auto-downloads kernel if needed)
eph = DEEphemeris(version='DE440')

# Query Sun position at J2000.0.  The Sun orbits the SSB within
# a couple of solar radii, so its SSB distance is small.
jd = 2451545.0  # Julian Date (TT)
r_sun, v_sun = eph.sun_position(jd)

print(f"Sun distance from SSB: {np.linalg.norm(r_sun):.6f} AU")
print(f"Sun speed: {np.linalg.norm(v_sun):.9f} AU/day")

Output:

Sun distance from SSB: 0.007668 AU
Sun speed: 0.000009154 AU/day

Available Bodies

The ephemeris supports queries for:

  • Sun - Position relative to Solar System Barycenter (~0.007 AU offset)

  • Moon - Position relative to SSB or Earth-centered

  • Planets - Mercury, Venus, Mars, Jupiter, Saturn, Uranus, Neptune

  • Barycenters - For any supported body

# Query planet positions
r_mars, v_mars = eph.planet_position('mars', jd)

# Moon geocentric position
r_moon_ec, v_moon_ec = eph.moon_position(jd, frame='earth_centered')

# Any body relative to SSB
r_body, v_body = eph.barycenter_position('mars', jd)

Ephemeris Versions

Supported DE versions with coverage:

  • DE440 (latest, 2020) - Covers 1550-2650, highest precision

  • DE432s (2013) - Covers 1350-3000, high precision for long-term

  • DE430 (2013) - Covers 1550-2650

  • DE405 (1998) - Covers 1600-2200, compact size

# Use a specific ephemeris version
eph_430 = DEEphemeris(version='DE430')
r, v = eph_430.sun_position(jd)

Frame Support

Positions can be returned in different frames:

# ICRF (default) - International Celestial Reference Frame
r_icrf, v_icrf = eph.sun_position(jd, frame='icrf')

# Ecliptic - J2000.0 ecliptic plane
r_ecliptic, v_ecliptic = eph.sun_position(jd, frame='ecliptic')

# Earth-centered (Moon only)
r_ec, v_ec = eph.moon_position(jd, frame='earth_centered')

Kernel Files

The library automatically downloads JPL ephemeris kernels (~100-120 MB) on first use. They are cached in ~/.jplephem/ for subsequent accesses.

To manually download kernels:

# Download a specific kernel
mkdir -p ~/.jplephem
cd ~/.jplephem
wget https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/de440.bsp

Module-Level Functions

Convenience functions are available for quick queries:

from pytcl.astronomical.ephemerides import (
    sun_position,
    moon_position,
    planet_position,
    barycenter_position,
)

# Use default DE440 ephemeris
r_sun, v_sun = sun_position(jd)
r_mars, v_mars = planet_position('mars', jd)

Relativistic Corrections

The relativity module provides functions for computing relativistic effects in orbital mechanics. The functions take a gravitational parameter \(GM\) (in m^3/s^2) rather than a mass; GM_EARTH and GM_SUN constants are provided, and GM_EARTH is the default for most functions.

Basic Schwarzschild Metrics

from pytcl.astronomical.relativity import (
    schwarzschild_radius,
    gravitational_time_dilation,
)

# Schwarzschild radius (event horizon) takes a mass in kg
M_sun = 1.989e30  # kg
Rs = schwarzschild_radius(M_sun)
print(f"Sun's Schwarzschild radius: {Rs:.2f} m")

# Time dilation factor at distance r from Earth's center
r = 7e6  # meters
factor = gravitational_time_dilation(r)  # gm defaults to GM_EARTH
print(f"Time dilation factor: {factor:.12f}")

Output:

Sun's Schwarzschild radius: 2954.13 m
Time dilation factor: 0.999999999366

Orbital Precession

Relativistic perihelion precession calculations:

from pytcl.astronomical.relativity import (
    GM_SUN,
    schwarzschild_precession_per_orbit,
)

# Mercury around the Sun
a = 5.79e10  # semi-major axis (m)
e = 0.2056   # eccentricity

precession = schwarzschild_precession_per_orbit(a, e, gm=GM_SUN)
print(f"Precession: {precession:.6e} rad/orbit")

Output:

Precession: 5.019383e-07 rad/orbit

For Mercury this accumulates to the famous ~43 arcsec/century.

GPS Time Effects

proper_time_rate(v, r) combines special relativistic (velocity) and general relativistic (gravity) time dilation, relative to a distant observer at rest. The familiar GPS number is the difference between the satellite clock and a clock on the ground:

from pytcl.astronomical.relativity import proper_time_rate

# GPS satellite: ~3.87 km/s at r ~ 26,600 km
rate_sat = proper_time_rate(3874.0, 2.66e7)

# Ground clock: equatorial rotation speed at Earth's surface
rate_ground = proper_time_rate(465.0, 6.371e6)

seconds_per_day = 86400
offset = (rate_sat / rate_ground - 1.0) * seconds_per_day
print(f"GPS clock gain: {offset*1e6:.1f} microseconds/day")

Output:

GPS clock gain: 38.6 microseconds/day

(The satellite clock runs fast by ~45 us/day from gravity and slow by ~7 us/day from velocity, for a net gain of ~38.6 us/day.)

Light Propagation Effects

Shapiro delay - the extra light travel time in a gravitational field. The function takes observer, source, and gravitating-body positions:

from pytcl.astronomical.relativity import shapiro_delay

# Earth-to-Mars signal passing near the Sun (superior conjunction)
AU = 1.495978707e11
observer = np.array([AU, 1e9, 0.0])         # Earth (m)
source = np.array([-1.52 * AU, -1e9, 0.0])  # Mars, opposite side (m)
sun = np.zeros(3)

delay = shapiro_delay(observer, source, sun)  # gm defaults to GM_SUN
print(f"Shapiro delay: {delay*1e6:.1f} microseconds")

Output:

Shapiro delay: 147.5 microseconds

Post-Newtonian Effects

post_newtonian_acceleration returns the total acceleration (Newtonian plus the first post-Newtonian Schwarzschild correction):

from pytcl.astronomical.relativity import post_newtonian_acceleration

r_vec = np.array([1e11, 0.0, 0.0])     # heliocentric position (m)
v_vec = np.array([0.0, 30000.0, 0.0])  # velocity (m/s)

a_total = post_newtonian_acceleration(r_vec, v_vec, gm=GM_SUN)
a_newt = -GM_SUN / np.linalg.norm(r_vec) ** 3 * r_vec

print(f"Total acceleration: {a_total} m/s^2")
print(f"1PN correction magnitude: {np.linalg.norm(a_total - a_newt):.3e} m/s^2")

Output:

Total acceleration: [-0.01327124  0.          0.        ] m/s^2
1PN correction magnitude: 6.510e-10 m/s^2

Geodetic Precession (De Sitter Effect)

Precession per orbit of a gyroscope transported around a central mass:

from pytcl.astronomical.relativity import geodetic_precession

# GPS-like orbit around Earth
a_orbit = 2.66e7  # semi-major axis (m)
e_orbit = 0.01
inc = np.deg2rad(55.0)

dw = geodetic_precession(a_orbit, e_orbit, inc)  # gm defaults to GM_EARTH
print(f"Geodetic precession: {dw:.3e} rad/orbit")

Output:

Geodetic precession: 1.572e-09 rad/orbit

Frame-Dragging (Lense-Thirring Effect)

Nodal precession rate caused by the spin of the central body:

from pytcl.astronomical.relativity import lense_thirring_precession

J_earth = 5.86e33  # Earth's spin angular momentum (kg m^2/s)
rate = lense_thirring_precession(a_orbit, e_orbit, inc, J_earth)
print(f"Lense-Thirring nodal rate: {rate:.3e} rad/s")

Output:

Lense-Thirring nodal rate: 4.625e-16 rad/s

Orbital Mechanics

The orbital mechanics module works with StateVector (position/velocity) and OrbitalElements named tuples. Distances are in km, velocities in km/s, and the gravitational parameter defaults to Earth’s (398600.4418 km^3/s^2):

from pytcl.astronomical.orbital_mechanics import (
    StateVector,
    state_to_orbital_elements,
    orbital_elements_to_state,
    kepler_propagate_state,
)

r0 = np.array([7000.0, 0.0, 0.0])  # position (km)
v0 = np.array([0.0, 7.546, 0.0])   # velocity (km/s)
state = StateVector(r=r0, v=v0)

# Convert state to classical orbital elements (a, e, i, raan, omega, nu)
elements = state_to_orbital_elements(state)
print(f"a = {elements.a:.1f} km, e = {elements.e:.4f}")

# Two-body propagation via Kepler's equation
state_1h = kepler_propagate_state(state, 3600.0)
print(f"Position after 1 h: {state_1h.r.round(1)} km")

Output:

a = 6999.9 km, e = 0.0000
Position after 1 h: [-5172.3 -4716.5     0. ] km

kepler_propagate propagates OrbitalElements directly, and orbital_elements_to_state converts elements back to a state vector.

Reference Frame Transformations

The reference frames module provides coordinate system conversions:

from pytcl.astronomical.reference_frames import (
    equatorial_to_ecliptic,
    mean_obliquity_iau80,
    precession_matrix_iau76,
)

# Equatorial to ecliptic requires the obliquity of the ecliptic
jd = 2451545.0
eps = mean_obliquity_iau80(jd)
r_equatorial = np.array([1.0, 0.0, 0.0])
r_ecliptic = equatorial_to_ecliptic(r_equatorial, eps)

# IAU 1976 precession from J2000.0 to a target epoch
jd_future = jd + 36525.0  # 100 years later
P = precession_matrix_iau76(jd_future)
r_precessed = P @ r_equatorial

Additional transformations include gcrf_to_itrf / itrf_to_gcrf, teme_to_itrf, eci_to_ecef / ecef_to_eci, nutation and polar motion matrices, and sidereal time functions (gmst_iau82, gast_iau82).

Time Systems

Astronomical time conversions and calculations:

from pytcl.astronomical.time_systems import (
    jd_to_mjd,
    mjd_to_jd,
    cal_to_jd,
    jd_to_cal,
)

# Julian Date conversions
jd = 2451545.0  # J2000.0 epoch
mjd = jd_to_mjd(jd)
print(f"MJD: {mjd}")

# Calendar conversions
jd2 = cal_to_jd(2000, 1, 1, 12, 0, 0.0)
print(f"JD for 2000-01-01 12:00: {jd2}")
print(f"Back to calendar: {jd_to_cal(jd2)}")

Output:

MJD: 51544.5
JD for 2000-01-01 12:00: 2451545.0
Back to calendar: (2000, 1, 1, 12, 0, 0.0)

Time scale conversions between UTC, TAI, TT, GPS, and Unix time are also available (utc_to_tai, tai_to_tt, utc_to_gps, unix_to_jd, and friends), along with sidereal time (gmst, gast).

Lambert Problem

Solve for orbits connecting two positions. Positions are in km, the time of flight in seconds, and the solvers return a LambertSolution named tuple with fields v1, v2, a, e, and tof:

from pytcl.astronomical.lambert import lambert_universal, lambert_izzo

r1 = np.array([7000.0, 0.0, 0.0])  # km
r2 = np.array([0.0, 7000.0, 0.0])  # km
tof = 2700.0                       # seconds

sol = lambert_universal(r1, r2, tof)  # mu defaults to Earth's, km^3/s^2
print(f"Initial velocity: {sol.v1.round(3)} km/s")
print(f"Final velocity: {sol.v2.round(3)} km/s")

Output:

Initial velocity: [3.585 5.964 0.   ] km/s
Final velocity: [-5.964 -3.585 -0.   ] km/s

See Also