Astronomical & Celestial Mechanics

Overview

The Tracker Component Library provides orbital mechanics, ephemeris, and celestial reference frame functionality for satellite operations, space vehicle navigation, and astronomical calculations.

Key Modules:

  • astronomical.orbital_mechanics - Keplerian elements, anomaly conversions, propagation

  • astronomical.reference_frames - ECEF/ECI, GCRF/ITRF/TEME, precession, nutation, polar motion

  • astronomical.ephemerides - Sun, Moon, planet positions (JPL DE ephemeris via jplephem)

  • astronomical.sgp4 - SGP4/SDP4 propagator (wraps the sgp4 package)

  • astronomical.tle - Two-Line Element parsing and formatting

  • astronomical.lambert - Lambert’s problem and orbit transfers

  • astronomical.relativity - Relativistic corrections

  • astronomical.time_systems - Julian dates, GMST, UTC/TAI/TT/GPS conversions

All functions are re-exported at the pytcl.astronomical package level.

Note

The orbital mechanics, Lambert, and SGP4 APIs work in kilometers and km/s (the default gravitational parameter is mu = 398600.4418 km^3/s^2), unlike the coordinate-conversion modules, which work in meters.

Orbital Elements and Representations

Keplerian Elements (6 parameters describing an orbit) are held in the OrbitalElements named tuple:

import numpy as np
from pytcl.astronomical import OrbitalElements, orbital_elements_to_state

elements = OrbitalElements(
    a=6700.0,                # semi-major axis (km)
    e=0.001,                 # eccentricity (0 = circular)
    i=np.radians(98.2),      # inclination (rad, sun-sync here)
    raan=np.radians(45.0),   # right ascension of ascending node (rad)
    omega=np.radians(30.0),  # argument of perigee (rad)
    nu=0.0,                  # true anomaly (rad, 0 = perigee)
)

# Convert to a Cartesian state (StateVector with .r and .v)
state = orbital_elements_to_state(elements)
print(f"Position (km):   {state.r}")
print(f"Velocity (km/s): {state.v}")

Output:

Position (km):   [4436.31508215 3761.26976554 3312.43462317]
Velocity (km/s): [-2.05538115 -3.40409261  6.61810164]

Anomaly Conversions

In orbital mechanics, there are three ways to measure position along the orbit:

  1. True Anomaly (nu): Actual angle from focus

  2. Eccentric Anomaly (E): Auxiliary angle (used in equations)

  3. Mean Anomaly (M): Time-based angle (increases uniformly with time)

from pytcl.astronomical import (
    true_to_eccentric_anomaly, eccentric_to_true_anomaly,
    eccentric_to_mean_anomaly, mean_to_eccentric_anomaly,
)

nu = 0.1   # true anomaly (radians)
e = 0.001  # eccentricity

# True -> Eccentric
E = true_to_eccentric_anomaly(nu, e)

# Eccentric -> Mean
M = eccentric_to_mean_anomaly(E, e)

# Mean -> Eccentric (Newton-Raphson iteration)
E_recovered = mean_to_eccentric_anomaly(M, e)

# Eccentric -> True
nu_recovered = eccentric_to_true_anomaly(E_recovered, e)

print(f"Original nu={nu:.6f}, recovered nu={nu_recovered:.6f}")

Hyperbolic and parabolic counterparts exist as well (mean_to_hyperbolic_anomaly, true_anomaly_to_parabolic_anomaly, etc.).

State Conversions

from pytcl.astronomical import StateVector, state_to_orbital_elements

# Cartesian state to Keplerian elements
state = StateVector(
    r=np.array([6600.0, 0.0, 0.0]),  # km
    v=np.array([0.0, 7.8, 0.0]),     # km/s
)
elements = state_to_orbital_elements(state)
print(f"a={elements.a:.1f} km, e={elements.e:.4f}")

# Keplerian back to Cartesian
state_recovered = orbital_elements_to_state(elements)

Keplerian Propagation

Basic Circular Orbit

from pytcl.astronomical import kepler_propagate, orbital_period

# Initial orbital elements: LEO polar orbit
elements = OrbitalElements(a=6700.0, e=0.0, i=np.pi/2,
                           raan=0.0, omega=0.0, nu=0.0)

period = orbital_period(elements.a)  # seconds
print(f"Orbital period: {period / 60:.1f} minutes")

# Propagate 100 seconds (two-body motion, no perturbations)
elements_prop = kepler_propagate(elements, 100.0)

# Check: semi-major axis is unchanged
print(f"a unchanged: {np.isclose(elements.a, elements_prop.a)}")

Output:

Orbital period: 91.0 minutes
a unchanged: True

kepler_propagate_state does the same starting from a StateVector.

With Mean Anomaly Propagation

from pytcl.astronomical import mean_motion

# For efficient long-term propagation, advance the mean anomaly
a = 6700.0
e = 0.001
n = mean_motion(a)  # rad/s

# Initial state
nu0 = 0.0
E0 = true_to_eccentric_anomaly(nu0, e)
M0 = eccentric_to_mean_anomaly(E0, e)

# Propagate mean anomaly
t_prop = 3600.0  # 1 hour
M_prop = M0 + n * t_prop

# Convert back to true anomaly
E_prop = mean_to_eccentric_anomaly(M_prop, e)
nu_prop = eccentric_to_true_anomaly(E_prop, e)

print(f"True anomaly after 1 hour: {np.degrees(nu_prop):.2f} deg")

Output:

True anomaly after 1 hour: 237.36 deg

Perturbations

The Keplerian propagator is a pure two-body model. For perturbed motion the library offers SGP4 (below), which models J2-J4, drag, and resonance effects for TLE-based satellites. The secular J2 rates themselves are a one-line formula if you need them directly, as in the sun-synchronous design example at the end of this page.

Atmospheric Drag (for low Earth orbits)

pytcl.atmosphere.simplified_thermosphere provides upper-atmosphere density. It is a simplified barometric model with solar-flux and geomagnetic coupling – not NRLMSISE-00 – and is only meaningful above roughly 200 km.

from pytcl.atmosphere import simplified_thermosphere

# Density at 400 km altitude (angles in radians, altitude in meters)
state_atm = simplified_thermosphere(
    latitude=np.radians(40.0),
    longitude=np.radians(-74.0),
    altitude=400e3,
    year=2026,
    day_of_year=57,
    seconds_in_day=43200.0,
    f107=150.0,   # 10.7 cm solar flux (SFU)
    ap=4.0,       # geomagnetic index
)
print(f"Density at 400 km: {state_atm.density:.3e} kg/m^3")

def drag_acceleration(rho, vel, area, mass, cd=2.2):
    """Cannonball drag: a = -0.5 (cd A / m) rho |v| v."""
    v_mag = np.linalg.norm(vel)
    return -0.5 * (cd * area / mass) * rho * v_mag * vel

vel = np.array([7500.0, 0.0, 0.0])  # m/s
a_drag = drag_acceleration(state_atm.density, vel, area=10.0, mass=500.0)

Output:

Density at 400 km: 2.886e-12 kg/m^3

Reference Frames: ECEF and ECI

Earth-Centered Inertial (ECI) vs Earth-Centered Earth-Fixed (ECEF)

  • ECI: Fixed to distant stars (inertial). Doesn’t rotate with Earth.

  • ECEF: Rotates with Earth. Fixed to ground stations.

The simple conversion is a rotation about the z-axis by the Greenwich Mean Sidereal Time (GMST):

import numpy as np
from pytcl.astronomical import cal_to_jd, gmst
from pytcl.astronomical.reference_frames import ecef_to_eci, eci_to_ecef

# Satellite position in ECI (km); positions only, shape (3,)
sat_eci = np.array([6600.0, 0.0, 0.0])

# Time of observation -> Julian date -> GMST (radians)
jd = cal_to_jd(2026, 2, 26, 12, 0, 0.0)
theta = gmst(jd)

# Transform to ECEF (what ground stations see)
sat_ecef = eci_to_ecef(sat_eci, theta)

print(f"ECI position:  {sat_eci}")
print(f"ECEF position: {sat_ecef}")

# Transform back
sat_eci_recovered = ecef_to_eci(sat_ecef, theta)

Output:

ECI position:  [6600.    0.    0.]
ECEF position: [6045.68286597 2647.58733268    0.        ]

Components of the Full Transformation

A high-precision ECI/ECEF transformation includes:

  1. Precession - Long-term wobble of Earth’s spin axis (~26,000 year period)

  2. Nutation - Short-term oscillation (~18.6 year period)

  3. Earth Rotation (GMST/GAST) - Daily rotation

  4. Polar Motion - Shift of the rotation axis (~1 meter, from EOP data)

Each component is available as a matrix:

from pytcl.astronomical import (
    precession_matrix_iau76, nutation_matrix, polar_motion_matrix,
)

jd = cal_to_jd(2026, 2, 26, 12, 0, 0.0)

theta = gmst(jd)
print(f"GMST: {theta:.4f} rad = {np.degrees(theta):.2f} deg")

P = precession_matrix_iau76(jd)      # IAU 1976 precession
N = nutation_matrix(jd)              # IAU 1980 nutation
W = polar_motion_matrix(xp=1e-6, yp=2e-6)  # pole coordinates (rad)

Output:

GMST: 5.8704 rad = 336.35 deg

High-Precision Transformations (GCRF/ITRF)

The full chain (precession + nutation + sidereal rotation + polar motion) is packaged as gcrf_to_itrf / itrf_to_gcrf. Supply UT1 and TT Julian dates and, if available, the pole coordinates from IERS EOP data:

from pytcl.astronomical import gcrf_to_itrf, itrf_to_gcrf

r_gcrf = np.array([6600.0, 0.0, 0.0])  # km
jd_ut1 = cal_to_jd(2026, 2, 26, 12, 0, 0.0)
jd_tt = jd_ut1 + 69.184 / 86400  # TT - UTC as of 2026

r_itrf = gcrf_to_itrf(r_gcrf, jd_ut1=jd_ut1, jd_tt=jd_tt)
r_back = itrf_to_gcrf(r_itrf, jd_ut1=jd_ut1, jd_tt=jd_tt)

teme_to_itrf / itrf_to_teme convert between the TEME frame that SGP4 outputs and the Earth-fixed ITRF frame.

Ephemeris Functions

The ephemeris functions use a JPL Development Ephemeris kernel through the jplephem package (pip install nrl-tracker[astronomy]). Positions are in AU and, in the 'icrf' frame, relative to the solar system barycenter (SSB), not to Earth.

Sun and Moon Positions

from pytcl.astronomical import cal_to_jd, sun_position, moon_position

jd = cal_to_jd(2026, 2, 26, 12, 0, 0.0)  # Julian date (TT)

# Sun position relative to the solar system barycenter (AU, ICRF)
sun_pos, sun_vel = sun_position(jd)
print(f"Sun position (AU, barycentric ICRF): {sun_pos}")

# Moon position relative to Earth's center
moon_pos, moon_vel = moon_position(jd, frame='earth_centered')
AU_KM = 149597870.7
print(f"Earth-Moon distance: {np.linalg.norm(moon_pos) * AU_KM:.0f} km")

Output:

Sun position (AU, barycentric ICRF): [-0.00266542 -0.00510309 -0.00207743]
Earth-Moon distance: 366214 km

Planet Positions

from pytcl.astronomical import planet_position

# Available: mercury, venus, mars, jupiter, saturn, uranus, neptune
for planet in ['mars', 'jupiter']:
    pos, vel = planet_position(planet, jd)
    print(f"{planet.capitalize()}: {np.linalg.norm(pos):.2f} AU from the SSB")

Output:

Mars: 1.39 AU from the SSB
Jupiter: 5.23 AU from the SSB

For repeated queries, instantiate DEEphemeris once and call its methods; the convenience functions above create and cache one internally.

SGP4 Propagator (TLE-Based)

From Two-Line Element (TLE)

parse_tle returns a TLE named tuple (angles in radians, mean motion in radians/minute); sgp4_propagate takes the TLE and the time since epoch in minutes and returns position/velocity in the TEME frame:

from pytcl.astronomical import (
    parse_tle, sgp4_propagate, tle_epoch_to_datetime,
)

# Example TLE (International Space Station)
tle_line1 = "1 25544U 98067A   24001.50000000  .00016717  00000-0  10270-3 0  9997"
tle_line2 = "2 25544  51.6400 247.4627 0006703 130.5360 325.0288 15.49815350479003"

tle = parse_tle(tle_line1, tle_line2, name="ISS (ZARYA)")

print(f"Catalog number: {tle.catalog_number}")
print(f"Epoch: {tle_epoch_to_datetime(tle)}")
rev_per_day = tle.mean_motion * 1440 / (2 * np.pi)
print(f"Mean motion: {rev_per_day:.4f} revolutions/day")

# Propagate 10 minutes past epoch
state = sgp4_propagate(tle, 10.0)
print(f"Position (TEME, km):   {state.r}")
print(f"Velocity (TEME, km/s): {state.v}")

Output:

Catalog number: 25544
Epoch: 2024-01-01 12:00:00+00:00
Mean motion: 15.4982 revolutions/day
Position (TEME, km):   [4602.28257877 3229.53092379 3800.56614424]
Velocity (TEME, km/s): [-0.96603157  6.34040043 -4.2050613 ]

Note

parse_tle verifies the TLE line checksums by default, so hand-edited lines are rejected (pass verify_checksum=False to skip the check).

Batch Propagation

from pytcl.astronomical import sgp4_propagate_batch

times = np.arange(0.0, 90.0, 10.0)  # minutes since epoch
positions, velocities = sgp4_propagate_batch(tle, times)
# positions: (9, 3) km, velocities: (9, 3) km/s

Lambert’s Problem

Orbit Transfers (from one orbit to another)

from pytcl.astronomical import lambert_universal, hohmann_transfer

# Initial orbit: LEO (km); final orbit: GEO
r1 = 6700.0
r2 = 42164.0

# Hohmann transfer delta-v and time of flight
dv1, dv2, tof = hohmann_transfer(r1, r2)
print(f"Hohmann: dv1={dv1:.3f} km/s, dv2={dv2:.3f} km/s, "
      f"tof={tof / 3600:.2f} hours")

# Solve Lambert's problem for a 90-degree transfer in the same time
r_initial = np.array([r1, 0.0, 0.0])
r_final = np.array([0.0, r2, 0.0])

sol = lambert_universal(r_initial, r_final, tof)
print(f"Departure velocity: {sol.v1} km/s")
print(f"Arrival velocity:   {sol.v2} km/s")

Output:

Hohmann: dv1=2.420 km/s, dv2=1.465 km/s, tof=5.28 hours
Departure velocity: [6.73474339 7.50479386 0.        ] km/s
Arrival velocity:   [-1.19253673 -0.42248626 -0.        ] km/s

lambert_izzo implements the same interface with Izzo’s algorithm, and bi_elliptic_transfer / minimum_energy_transfer cover the other classical transfer computations.

Relativistic Corrections

Clock Rate (Proper Time) (main relativistic correction for GPS/timing)

from pytcl.astronomical.relativity import proper_time_rate

# GPS satellite: orbital speed (m/s) and geocentric radius (m)
rate_sat = proper_time_rate(v=3874.0, r=26560e3)

# Ground clock on the rotating equator
rate_ground = proper_time_rate(v=465.0, r=6378e3)

# GPS satellite clocks run ~38 microseconds/day faster than ground clocks
delta = rate_sat - rate_ground
print(f"Fractional rate difference: {delta:.2e}")
print(f"Clock offset: {delta * 86400 * 1e6:.1f} microseconds/day")

Output:

Fractional rate difference: 4.46e-10
Clock offset: 38.5 microseconds/day

Shapiro Delay (photon traveling through spacetime curvature)

from pytcl.astronomical.relativity import shapiro_delay, GM_EARTH

# Ground station position (ECEF, equator)
station = np.array([6378e3, 0.0, 0.0])  # meters

# GPS satellite position
sat = np.array([26560e3, 0.0, 0.0])  # meters

# Signal delay due to Earth's spacetime curvature
# (observer, light source, gravitating body position, GM)
delay = shapiro_delay(station, sat, np.zeros(3), GM_EARTH)

print(f"Shapiro delay: {delay * 1e12:.1f} picoseconds")

Output:

Shapiro delay: 42.2 picoseconds

Relativistic Range Correction (Shapiro delay as an equivalent range)

from pytcl.astronomical.relativity import relativistic_range_correction

r1 = 6378e3          # Station geocentric radius (m)
r2 = 26560e3         # Satellite geocentric radius (m)
rho = r2 - r1        # Station-satellite range (m)

dr = relativistic_range_correction(r1, r2, rho, GM_EARTH)
print(f"Range correction: {dr * 1000:.1f} mm")

Output:

Range correction: 12.7 mm

Complete Satellite Tracking Example

Track a Satellite from a Ground Station

import datetime
import numpy as np
from pytcl.astronomical import (
    cal_to_jd, parse_tle, sgp4_propagate, teme_to_itrf,
    tle_epoch_to_datetime,
)
from pytcl.coordinate_systems.conversions import ecef2enu, geodetic2ecef

class GroundTracker:
    def __init__(self, tle, station_lat, station_lon, station_alt):
        """
        Args:
            tle: parsed TLE (from parse_tle)
            station_lat, station_lon: ground station (radians)
            station_alt: ground station altitude (meters)
        """
        self.tle = tle
        self.epoch = tle_epoch_to_datetime(tle)
        self.lat = station_lat
        self.lon = station_lon
        self.station_ecef = geodetic2ecef(station_lat, station_lon,
                                          station_alt)

    def get_topocentric(self, time):
        """Satellite azimuth/elevation/range as seen from the station."""
        minutes = (time - self.epoch).total_seconds() / 60.0
        state = sgp4_propagate(self.tle, minutes)

        # SGP4 output is TEME; rotate to Earth-fixed and convert to m
        jd = cal_to_jd(time.year, time.month, time.day,
                       time.hour, time.minute, time.second)
        pos_ecef = teme_to_itrf(state.r, jd) * 1000.0

        enu = ecef2enu(pos_ecef, self.lat, self.lon, self.station_ecef)
        east, north, up = enu
        range_m = np.linalg.norm(enu)
        azimuth = np.arctan2(east, north)  # from north
        elevation = np.arctan2(up, np.hypot(east, north))
        return np.degrees(azimuth), np.degrees(elevation), range_m

# Example: Track the ISS from NYC
tle1 = "1 25544U 98067A   24001.50000000  .00016717  00000-0  10270-3 0  9997"
tle2 = "2 25544  51.6400 247.4627 0006703 130.5360 325.0288 15.49815350479003"
tle = parse_tle(tle1, tle2)

tracker = GroundTracker(tle, np.radians(40.7128), np.radians(-74.0060),
                        10.0)

# Sample one orbit starting at the TLE epoch
t_start = tle_epoch_to_datetime(tle)
for i in range(0, 90, 15):
    t = t_start + datetime.timedelta(minutes=i)
    az, el, rng = tracker.get_topocentric(t)
    print(f"{i:3d} min: Az={az:7.1f} deg, El={el:6.1f} deg, "
          f"Range={rng / 1000:8.1f} km")

Output (negative elevation means the satellite is below the horizon):

 0 min: Az=   24.1 deg, El= -38.6 deg, Range=  8608.8 km
15 min: Az=  -21.0 deg, El= -56.9 deg, Range= 11188.4 km
30 min: Az=  -87.5 deg, El= -65.3 deg, Range= 12041.4 km
45 min: Az= -145.7 deg, El= -51.8 deg, Range= 10528.0 km
60 min: Az=  175.5 deg, El= -29.8 deg, Range=  7100.4 km
75 min: Az=  109.9 deg, El= -14.1 deg, Range=  4373.8 km

Performance Considerations

Propagation Speed

import time

# Keplerian propagation: very fast (two-body only)
elements = OrbitalElements(a=6700.0, e=0.001, i=1.0,
                           raan=0.0, omega=0.0, nu=0.0)
t0 = time.perf_counter()
for _ in range(10000):
    kepler_propagate(elements, 100.0)
print(f"Keplerian: {time.perf_counter() - t0:.3f} s for 10k propagations")

# SGP4: slower but includes perturbations
t0 = time.perf_counter()
for _ in range(10000):
    sgp4_propagate(tle, 1.0)
print(f"SGP4: {time.perf_counter() - t0:.3f} s for 10k propagations")

For many time steps of a single satellite, sgp4_propagate_batch amortizes the per-call overhead.

Accuracy vs Speed Tradeoff

Method              Accuracy     Speed      Best For

Keplerian           ~1 km        Fastest    Circular orbits, quick estimates
SGP4                ~1-3 km      Medium     TLE-based (public data)
Numerical + drag    Best         Slow       Precise long-term (GPS, research)

Common Orbital Scenarios

Sun-Synchronous Orbit Design

For a sun-synchronous orbit the J2-driven RAAN precession must equal 360 degrees per year. Solve the secular J2 rate equation for inclination:

from scipy.optimize import fsolve

J2 = 1.08262668e-3
Re = 6378.137        # km
mu = 398600.4418     # km^3/s^2

def design_sun_sync_orbit(altitude_km):
    """Find the inclination whose J2 nodal precession tracks the Sun."""
    a = Re + altitude_km
    n = np.sqrt(mu / a**3)
    target = 2 * np.pi / (365.25 * 86400)  # rad/s

    def sun_sync_condition(i):
        # Secular RAAN rate for a circular orbit:
        # dRAAN/dt = -1.5 n J2 (Re/a)^2 cos(i)
        return -1.5 * n * J2 * (Re / a) ** 2 * np.cos(i) - target

    return np.degrees(fsolve(sun_sync_condition, np.radians(98.0))[0])

print(f"Sun-sync inclination at 700 km: {design_sun_sync_orbit(700.0):.2f} deg")

Output:

Sun-sync inclination at 700 km: 98.19 deg

See Also