Astronomical

Orbital mechanics and astronomical computations.

Astronomical calculations for target tracking.

This module provides time system conversions, orbital mechanics, Lambert problem solvers, reference frame transformations, and high-precision ephemerides for celestial bodies.

Examples

>>> from pytcl.astronomical import kepler_propagate, OrbitalElements
>>> import numpy as np

>>> # Propagate an orbit forward in time
>>> elements = OrbitalElements(a=7000, e=0.01, i=0.5, raan=0, omega=0, nu=0)
>>> new_elements = kepler_propagate(elements, 3600)  # 1 hour

>>> # Solve Lambert's problem
>>> from pytcl.astronomical import lambert_universal
>>> r1 = np.array([5000, 10000, 2100])
>>> r2 = np.array([-14600, 2500, 7000])
>>> sol = lambert_universal(r1, r2, 3600)

>>> # Query Sun position with high precision
>>> from pytcl.astronomical import sun_position
>>> r_sun, v_sun = sun_position(2451545.0)  # J2000.0

pytcl.astronomical.cal_to_jd(year, month, day, hour=0, minute=0, second=0.0)

Convert calendar date to Julian Date.

Parameters:

• year (int) – Year (negative for BCE).

• month (int) – Month (1-12).

• day (int) – Day of month (1-31).

• hour (int, optional) – Hour (0-23).

• minute (int, optional) – Minute (0-59).

• second (float, optional) – Second (0-59.999…).

Returns:

Julian Date.

Return type:

float

Examples

>>> cal_to_jd(2000, 1, 1, 12, 0, 0)  # J2000.0
2451545.0
>>> cal_to_jd(1970, 1, 1, 0, 0, 0)  # Unix epoch
2440587.5

Notes

Uses the algorithm from Meeus, “Astronomical Algorithms”, 2nd ed. Valid for dates after October 15, 1582 (Gregorian calendar).

pytcl.astronomical.jd_to_cal(jd)

Convert Julian Date to calendar date.

Parameters:

jd (float) – Julian Date.

Returns:

(year, month, day, hour, minute, second).

Return type:

tuple

Examples

>>> jd_to_cal(2451545.0)  # J2000.0
(2000, 1, 1, 12, 0, 0.0)

pytcl.astronomical.mjd_to_jd(mjd)

Convert Modified Julian Date to Julian Date.

Parameters:

mjd (float) – Modified Julian Date.

Returns:

Julian Date.

Return type:

float

Examples

>>> mjd = 44239.0  # 1980-01-01
>>> jd = mjd_to_jd(mjd)
>>> jd
2444239.5

Notes

MJD = JD - 2400000.5

pytcl.astronomical.jd_to_mjd(jd)

Convert Julian Date to Modified Julian Date.

Parameters:

jd (float) – Julian Date.

Returns:

Modified Julian Date.

Return type:

float

Examples

>>> jd = 2444239.5  # 1980-01-01
>>> mjd = jd_to_mjd(jd)
>>> mjd
44239.0

pytcl.astronomical.utc_to_tai(year, month, day, hour=0, minute=0, second=0.0)

Convert UTC to TAI (Julian Date).

Parameters:

• year (int) – Calendar date.

• month (int) – Calendar date.

• day (int) – Calendar date.

• hour (int, optional) – Time of day.

• minute (int, optional) – Time of day.

• second (float, optional) – Seconds (including fractional).

Returns:

TAI as Julian Date.

Return type:

float

Notes

TAI = UTC + leap_seconds

pytcl.astronomical.tai_to_utc(jd_tai)

Convert TAI (Julian Date) to UTC.

Parameters:

jd_tai (float) – TAI as Julian Date.

Returns:

• jd_utc (float) – UTC as Julian Date.

• leap_seconds (int) – Number of leap seconds applied.

Return type:

Tuple[float, int]

Notes

This is an approximate conversion that may have small errors near leap second boundaries.

pytcl.astronomical.tai_to_tt(jd_tai)

Convert TAI (Julian Date) to TT (Terrestrial Time).

Parameters:

jd_tai (float) – TAI as Julian Date.

Returns:

TT as Julian Date.

Return type:

float

Notes

TT = TAI + 32.184 seconds (constant offset)

pytcl.astronomical.tt_to_tai(jd_tt)

Convert TT (Terrestrial Time) to TAI (Julian Date).

Parameters:

jd_tt (float) – TT as Julian Date.

Returns:

TAI as Julian Date.

Return type:

float

pytcl.astronomical.utc_to_tt(year, month, day, hour=0, minute=0, second=0.0)

Convert UTC to TT (Terrestrial Time).

Parameters:

• year (int) – Calendar date.

• month (int) – Calendar date.

• day (int) – Calendar date.

• hour (int, optional) – Time of day.

• minute (int, optional) – Time of day.

• second (float, optional) – Seconds.

Returns:

TT as Julian Date.

Return type:

float

Notes

TT = UTC + leap_seconds + 32.184

pytcl.astronomical.tt_to_utc(jd_tt)

Convert TT to UTC.

Parameters:

jd_tt (float) – TT as Julian Date.

Returns:

• jd_utc (float) – UTC as Julian Date.

• leap_seconds (int) – Number of leap seconds applied.

Return type:

Tuple[float, int]

pytcl.astronomical.tai_to_gps(jd_tai)

Convert TAI to GPS time.

Parameters:

jd_tai (float) – TAI as Julian Date.

Returns:

GPS time as Julian Date.

Return type:

float

Notes

GPS time = TAI - 19 seconds (offset at GPS epoch) GPS time does not include leap seconds after 1980.

pytcl.astronomical.gps_to_tai(jd_gps)

Convert GPS time to TAI.

Parameters:

jd_gps (float) – GPS time as Julian Date.

Returns:

TAI as Julian Date.

Return type:

float

pytcl.astronomical.utc_to_gps(year, month, day, hour=0, minute=0, second=0.0)

Convert UTC to GPS time.

Parameters:

• year (int) – Calendar date.

• month (int) – Calendar date.

• day (int) – Calendar date.

• hour (int, optional) – Time of day.

• minute (int, optional) – Time of day.

• second (float, optional) – Seconds.

Returns:

GPS time as Julian Date.

Return type:

float

pytcl.astronomical.gps_to_utc(jd_gps)

Convert GPS time to UTC.

Parameters:

jd_gps (float) – GPS time as Julian Date.

Returns:

• jd_utc (float) – UTC as Julian Date.

• leap_seconds (int) – Number of leap seconds applied.

Return type:

Tuple[float, int]

pytcl.astronomical.unix_to_jd(unix_time)

Convert Unix timestamp to Julian Date.

Parameters:

unix_time (float) – Unix timestamp (seconds since 1970-01-01 00:00:00 UTC).

Returns:

Julian Date.

Return type:

float

Examples

>>> unix_to_jd(0.0)  # Unix epoch
2440587.5

pytcl.astronomical.jd_to_unix(jd)

Convert Julian Date to Unix timestamp.

Parameters:

jd (float) – Julian Date.

Returns:

Unix timestamp.

Return type:

float

Examples

>>> jd = 2440587.5  # 1970-01-01 00:00:00 UTC
>>> unix_to_jd = jd_to_unix(jd)
>>> unix_to_jd
0.0

pytcl.astronomical.gps_week_seconds(jd_gps)

Convert GPS Julian Date to GPS week and seconds of week.

Parameters:

jd_gps (float) – GPS time as Julian Date.

Returns:

• week (int) – GPS week number.

• seconds (float) – Seconds into the week.

Return type:

Tuple[int, float]

Examples

>>> gps_week_seconds(JD_GPS_EPOCH)
(0, 0.0)

pytcl.astronomical.gps_week_to_utc(week, seconds)

Convert GPS week and seconds to UTC.

Parameters:

• week (int) – GPS week number.

• seconds (float) – Seconds into the week.

Returns:

• jd_utc (float) – UTC as Julian Date.

• leap_seconds (int) – Number of leap seconds applied.

Return type:

Tuple[float, int]

pytcl.astronomical.gmst(jd_ut1)

Compute Greenwich Mean Sidereal Time.

Parameters:

jd_ut1 (float) – UT1 as Julian Date.

Returns:

GMST in radians.

Return type:

float

Notes

Uses the IAU 1982 expression for GMST.

References

[1]

Explanatory Supplement to the Astronomical Almanac, 3rd ed.

pytcl.astronomical.gast(jd_ut1, dpsi=0.0, eps=0.0)

Compute Greenwich Apparent Sidereal Time.

Parameters:

• jd_ut1 (float) – UT1 as Julian Date.

• dpsi (float, optional) – Nutation in longitude (radians). Default 0.

• eps (float, optional) – Mean obliquity of ecliptic (radians). Default uses approximate value.

Returns:

GAST in radians.

Return type:

float

Notes

GAST = GMST + equation of equinoxes The equation of equinoxes = dpsi * cos(eps)

For high precision applications, nutation parameters should be computed from the IAU 2006/2000A precession-nutation model.

pytcl.astronomical.get_leap_seconds(year, month, day)

Get the number of leap seconds (TAI-UTC offset) for a given date.

Parameters:

• year (int) – Year.

• month (int) – Month (1-12).

• day (int) – Day of month.

Returns:

TAI-UTC offset in seconds.

Return type:

int

Examples

>>> get_leap_seconds(2020, 1, 1)
37
>>> get_leap_seconds(1980, 1, 6)
19

class pytcl.astronomical.LeapSecondTable

Bases: object

Table of leap seconds for UTC-TAI conversion.

The leap second table contains the dates when leap seconds were added and the cumulative TAI-UTC offset.

entries

List of (year, month, day, tai_utc_offset) tuples.

Type:

list of tuple

Notes

This table must be updated when new leap seconds are announced. Last update: 2017-01-01 (37 seconds).

As of 2024, no new leap seconds have been added since 2017.

__init__()

get_offset(year, month, day)

Get TAI-UTC offset for a given date.

Parameters:

• year (int) – Year.

• month (int) – Month (1-12).

• day (int) – Day of month.

Returns:

TAI-UTC offset in seconds.

Return type:

int

class pytcl.astronomical.OrbitalElements(a, e, i, raan, omega, nu)

Bases: NamedTuple

Classical (Keplerian) orbital elements.

a

Semi-major axis (km). Negative for hyperbolic orbits.

Type:

float

e

Eccentricity. 0=circular, 0<e<1=elliptic, e=1=parabolic, e>1=hyperbolic.

Type:

float

i

Inclination (radians), 0 to pi.

Type:

float

raan

Right ascension of ascending node (radians), 0 to 2*pi.

Type:

float

omega

Argument of periapsis (radians), 0 to 2*pi.

Type:

float

nu

True anomaly (radians), 0 to 2*pi.

Type:

float

a: float

Alias for field number 0

e: float

Alias for field number 1

i: float

Alias for field number 2

raan: float

Alias for field number 3

omega: float

Alias for field number 4

nu: float

Alias for field number 5

class pytcl.astronomical.StateVector(r, v)

Bases: NamedTuple

Cartesian state vector.

r

Position vector (km), shape (3,).

Type:

ndarray

v

Velocity vector (km/s), shape (3,).

Type:

ndarray

r: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 0

v: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 1

pytcl.astronomical.mean_to_eccentric_anomaly(M, e, tol=1e-12, max_iter=50)

Solve Kepler’s equation: M = E - e*sin(E).

Uses Newton-Raphson iteration to find eccentric anomaly E given mean anomaly M and eccentricity e.

Parameters:

• M (float) – Mean anomaly (radians).

• e (float) – Eccentricity (0 <= e < 1 for elliptic orbits).

• tol (float, optional) – Convergence tolerance. Default 1e-12.

• max_iter (int, optional) – Maximum iterations. Default 50.

Returns:

E – Eccentric anomaly (radians).

Return type:

float

Raises:

ValueError – If eccentricity is not in valid range or iteration fails.

Examples

>>> import numpy as np
>>> E = mean_to_eccentric_anomaly(np.pi/4, 0.5)
>>> print(f"E = {np.degrees(E):.4f} degrees")

pytcl.astronomical.mean_to_hyperbolic_anomaly(M, e, tol=1e-12, max_iter=50)

Solve hyperbolic Kepler’s equation: M = e*sinh(H) - H.

Uses Newton-Raphson iteration to find hyperbolic anomaly H given mean anomaly M and eccentricity e.

Parameters:

• M (float) – Mean anomaly (radians).

• e (float) – Eccentricity (e > 1 for hyperbolic orbits).

• tol (float, optional) – Convergence tolerance. Default 1e-12.

• max_iter (int, optional) – Maximum iterations. Default 50.

Returns:

H – Hyperbolic anomaly (radians).

Return type:

float

Examples

>>> H = mean_to_hyperbolic_anomaly(1.0, 1.5)
>>> abs(1.5 * np.sinh(H) - H - 1.0) < 1e-10
True

pytcl.astronomical.eccentric_to_true_anomaly(E, e)

Convert eccentric anomaly to true anomaly.

Parameters:

• E (float) – Eccentric anomaly (radians).

• e (float) – Eccentricity.

Returns:

nu – True anomaly (radians), in [0, 2*pi).

Return type:

float

Examples

>>> nu = eccentric_to_true_anomaly(np.pi/4, 0.5)
>>> 0 <= nu < 2 * np.pi
True

pytcl.astronomical.true_to_eccentric_anomaly(nu, e)

Convert true anomaly to eccentric anomaly.

Parameters:

• nu (float) – True anomaly (radians).

• e (float) – Eccentricity.

Returns:

E – Eccentric anomaly (radians), in [0, 2*pi).

Return type:

float

Examples

>>> E = true_to_eccentric_anomaly(np.pi/3, 0.5)
>>> 0 <= E < 2 * np.pi
True

pytcl.astronomical.hyperbolic_to_true_anomaly(H, e)

Convert hyperbolic anomaly to true anomaly.

Parameters:

• H (float) – Hyperbolic anomaly (radians).

• e (float) – Eccentricity (e > 1).

Returns:

nu – True anomaly (radians).

Return type:

float

Examples

>>> nu = hyperbolic_to_true_anomaly(0.5, 1.5)
>>> isinstance(nu, float)
True

pytcl.astronomical.true_to_hyperbolic_anomaly(nu, e)

Convert true anomaly to hyperbolic anomaly.

Parameters:

• nu (float) – True anomaly (radians).

• e (float) – Eccentricity (e > 1).

Returns:

H – Hyperbolic anomaly (radians).

Return type:

float

pytcl.astronomical.eccentric_to_mean_anomaly(E, e)

Convert eccentric anomaly to mean anomaly (Kepler’s equation).

Parameters:

• E (float) – Eccentric anomaly (radians).

• e (float) – Eccentricity.

Returns:

M – Mean anomaly (radians).

Return type:

float

Examples

>>> M = eccentric_to_mean_anomaly(np.pi/4, 0.5)
>>> 0 <= M < 2 * np.pi
True

pytcl.astronomical.mean_to_true_anomaly(M, e)

Convert mean anomaly to true anomaly.

Parameters:

• M (float) – Mean anomaly (radians).

• e (float) – Eccentricity.

Returns:

nu – True anomaly (radians).

Return type:

float

Examples

>>> nu = mean_to_true_anomaly(np.pi/4, 0.1)
>>> 0 <= nu < 2 * np.pi
True

pytcl.astronomical.true_to_mean_anomaly(nu, e)

Convert true anomaly to mean anomaly.

Parameters:

• nu (float) – True anomaly (radians).

• e (float) – Eccentricity.

Returns:

M – Mean anomaly (radians).

Return type:

float

pytcl.astronomical.orbital_elements_to_state(elements, mu=398600.4418)

Convert orbital elements to Cartesian state vector.

Parameters:

• elements (OrbitalElements) – Classical orbital elements.

• mu (float, optional) – Gravitational parameter (km^3/s^2). Default is Earth.

Returns:

state – Position and velocity vectors in inertial frame.

Return type:

StateVector

Examples

>>> elements = OrbitalElements(a=7000, e=0.01, i=0.5, raan=0, omega=0, nu=0)
>>> state = orbital_elements_to_state(elements)
>>> print(f"r = {state.r}")

pytcl.astronomical.state_to_orbital_elements(state, mu=398600.4418)

Convert Cartesian state vector to orbital elements.

Parameters:

• state (StateVector) – Position and velocity vectors.

• mu (float, optional) – Gravitational parameter (km^3/s^2). Default is Earth.

Returns:

elements – Classical orbital elements.

Return type:

OrbitalElements

Examples

>>> r = np.array([7000, 0, 0])
>>> v = np.array([0, 7.5, 0])
>>> state = StateVector(r=r, v=v)
>>> elements = state_to_orbital_elements(state)

pytcl.astronomical.kepler_propagate(elements, dt, mu=398600.4418)

Propagate orbital elements using Kepler’s equation.

This performs two-body propagation by advancing the mean anomaly and solving Kepler’s equation for the new true anomaly.

Parameters:

• elements (OrbitalElements) – Initial orbital elements.

• dt (float) – Time step (seconds).

• mu (float, optional) – Gravitational parameter (km^3/s^2).

Returns:

new_elements – Propagated orbital elements.

Return type:

OrbitalElements

Examples

>>> elements = OrbitalElements(a=7000, e=0.01, i=0.5, raan=0, omega=0, nu=0)
>>> new_elements = kepler_propagate(elements, 3600)  # 1 hour

pytcl.astronomical.kepler_propagate_state(state, dt, mu=398600.4418)

Propagate state vector using Kepler’s equation.

Parameters:

• state (StateVector) – Initial state vector.

• dt (float) – Time step (seconds).

• mu (float, optional) – Gravitational parameter (km^3/s^2).

Returns:

new_state – Propagated state vector.

Return type:

StateVector

Examples

>>> r = np.array([7000.0, 0.0, 0.0])
>>> v = np.array([0.0, 7.5, 0.0])
>>> state = StateVector(r=r, v=v)
>>> new_state = kepler_propagate_state(state, 3600)
>>> np.linalg.norm(new_state.r) > 0
True

pytcl.astronomical.orbital_period(a, mu=398600.4418)

Compute orbital period for elliptic orbit.

Parameters:

• a (float) – Semi-major axis (km).

• mu (float, optional) – Gravitational parameter (km^3/s^2).

Returns:

T – Orbital period (seconds).

Return type:

float

Examples

>>> T = orbital_period(7000)  # LEO satellite
>>> T / 60  # Convert to minutes
97.8...

pytcl.astronomical.mean_motion(a, mu=398600.4418)

Compute mean motion.

Parameters:

• a (float) – Semi-major axis (km).

• mu (float, optional) – Gravitational parameter (km^3/s^2).

Returns:

n – Mean motion (radians/second).

Return type:

float

Examples

>>> n = mean_motion(42164)  # GEO orbit
>>> revs_per_day = n * 86400 / (2 * np.pi)
>>> abs(revs_per_day - 1.0) < 0.01  # Approximately 1 rev/day
True

pytcl.astronomical.vis_viva(r, a, mu=398600.4418)

Compute orbital velocity using vis-viva equation.

Parameters:

• r (float) – Current orbital radius (km).

• a (float) – Semi-major axis (km).

• mu (float, optional) – Gravitational parameter (km^3/s^2).

Returns:

v – Orbital velocity (km/s).

Return type:

float

Examples

>>> v = vis_viva(7000, 7000)  # Circular orbit
>>> abs(v - circular_velocity(7000)) < 0.01
True

pytcl.astronomical.specific_angular_momentum(state)

Compute specific angular momentum vector.

Parameters:

state (StateVector) – State vector.

Returns:

h – Specific angular momentum vector (km^2/s).

Return type:

ndarray

Examples

>>> r = np.array([7000.0, 0.0, 0.0])
>>> v = np.array([0.0, 7.5, 0.0])
>>> state = StateVector(r=r, v=v)
>>> h = specific_angular_momentum(state)
>>> h[2]  # Angular momentum in z-direction
52500.0

pytcl.astronomical.specific_orbital_energy(state, mu=398600.4418)

Compute specific orbital energy.

Parameters:

• state (StateVector) – State vector.

• mu (float, optional) – Gravitational parameter (km^3/s^2).

Returns:

energy – Specific orbital energy (km^2/s^2). Negative for bound orbits, positive for escape trajectories.

Return type:

float

Examples

>>> r = np.array([7000.0, 0.0, 0.0])
>>> v = np.array([0.0, 7.5, 0.0])
>>> state = StateVector(r=r, v=v)
>>> energy = specific_orbital_energy(state)
>>> energy < 0  # Bound orbit
True

pytcl.astronomical.flight_path_angle(state)

Compute flight path angle.

Parameters:

state (StateVector) – State vector.

Returns:

gamma – Flight path angle (radians). Positive when climbing, negative when descending.

Return type:

float

Examples

>>> r = np.array([7000.0, 0.0, 0.0])
>>> v = np.array([0.0, 7.5, 0.0])  # Tangential velocity
>>> state = StateVector(r=r, v=v)
>>> gamma = flight_path_angle(state)
>>> abs(gamma) < 0.01  # Nearly zero for circular motion
True

pytcl.astronomical.periapsis_radius(a, e)

Compute periapsis radius.

Parameters:

• a (float) – Semi-major axis (km).

• e (float) – Eccentricity.

Returns:

r_p – Periapsis radius (km).

Return type:

float

Examples

>>> r_p = periapsis_radius(10000, 0.3)
>>> r_p
7000.0

pytcl.astronomical.apoapsis_radius(a, e)

Compute apoapsis radius.

Parameters:

• a (float) – Semi-major axis (km).

• e (float) – Eccentricity.

Returns:

r_a – Apoapsis radius (km). Infinite for parabolic/hyperbolic orbits.

Return type:

float

Examples

>>> r_a = apoapsis_radius(10000, 0.3)
>>> r_a
13000.0

pytcl.astronomical.time_since_periapsis(nu, a, e, mu=398600.4418)

Compute time since periapsis passage.

Parameters:

• nu (float) – True anomaly (radians).

• a (float) – Semi-major axis (km).

• e (float) – Eccentricity.

• mu (float, optional) – Gravitational parameter (km^3/s^2).

Returns:

t – Time since periapsis (seconds).

Return type:

float

Examples

>>> t = time_since_periapsis(np.pi, 7000, 0.1)  # At apoapsis
>>> T = orbital_period(7000)
>>> abs(t - T/2) < 1  # Approximately half the period
True

pytcl.astronomical.orbit_radius(nu, a, e)

Compute orbital radius at given true anomaly.

Parameters:

• nu (float) – True anomaly (radians).

• a (float) – Semi-major axis (km).

• e (float) – Eccentricity.

Returns:

r – Orbital radius (km).

Return type:

float

Examples

>>> r = orbit_radius(0, 10000, 0.3)  # At periapsis
>>> r
7000.0
>>> r = orbit_radius(np.pi, 10000, 0.3)  # At apoapsis
>>> r
13000.0

pytcl.astronomical.escape_velocity(r, mu=398600.4418)

Compute escape velocity at given radius.

Parameters:

• r (float) – Radial distance (km).

• mu (float, optional) – Gravitational parameter (km^3/s^2).

Returns:

v_esc – Escape velocity (km/s).

Return type:

float

Examples

>>> v_esc = escape_velocity(6378 + 400)  # At ISS altitude
>>> 10 < v_esc < 12  # About 11 km/s
True

pytcl.astronomical.circular_velocity(r, mu=398600.4418)

Compute circular orbital velocity at given radius.

Parameters:

• r (float) – Orbital radius (km).

• mu (float, optional) – Gravitational parameter (km^3/s^2).

Returns:

v_circ – Circular velocity (km/s).

Return type:

float

Examples

>>> v_circ = circular_velocity(6378 + 400)  # At ISS altitude
>>> 7 < v_circ < 8  # About 7.7 km/s
True

class pytcl.astronomical.LambertSolution(v1, v2, a, e, tof)

Bases: NamedTuple

Solution to Lambert’s problem.

v1

Velocity at first position (km/s), shape (3,).

Type:

ndarray

v2

Velocity at second position (km/s), shape (3,).

Type:

ndarray

a

Semi-major axis of transfer orbit (km).

Type:

float

e

Eccentricity of transfer orbit.

Type:

float

tof

Time of flight (seconds).

Type:

float

v1: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 0

v2: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 1

a: float

Alias for field number 2

e: float

Alias for field number 3

tof: float

Alias for field number 4

pytcl.astronomical.lambert_universal(r1, r2, tof, mu=398600.4418, prograde=True, low_path=True, max_iter=100, tol=1e-10)

Solve Lambert’s problem using universal variables.

Given two position vectors and time of flight, determine the transfer orbit connecting them.

Parameters:

• r1 (ndarray) – Initial position vector (km), shape (3,).

• r2 (ndarray) – Final position vector (km), shape (3,).

• tof (float) – Time of flight (seconds). Must be positive.

• mu (float, optional) – Gravitational parameter (km^3/s^2). Default is Earth.

• prograde (bool, optional) – If True, use prograde (counterclockwise) transfer. If False, use retrograde transfer. Default True.

• low_path (bool, optional) – If True, use low energy (short way) transfer. If False, use high energy (long way) transfer. Default True.

• max_iter (int, optional) – Maximum iterations. Default 100.

• tol (float, optional) – Convergence tolerance. Default 1e-10.

Returns:

solution – Solution containing velocities and orbital parameters.

Return type:

LambertSolution

Raises:

ValueError – If solution does not converge.

Examples

>>> r1 = np.array([5000, 10000, 2100])  # km
>>> r2 = np.array([-14600, 2500, 7000])  # km
>>> tof = 3600  # 1 hour
>>> sol = lambert_universal(r1, r2, tof)
>>> print(f"v1 = {sol.v1} km/s")

pytcl.astronomical.lambert_izzo(r1, r2, tof, mu=398600.4418, prograde=True, multi_rev=0, max_iter=100, tol=1e-10)

Solve Lambert’s problem using Izzo’s algorithm.

This is a more robust algorithm that handles multi-revolution transfers and edge cases better than the universal variable method.

Parameters:

• r1 (ndarray) – Initial position vector (km), shape (3,).

• r2 (ndarray) – Final position vector (km), shape (3,).

• tof (float) – Time of flight (seconds).

• mu (float, optional) – Gravitational parameter (km^3/s^2). Default is Earth.

• prograde (bool, optional) – If True, use prograde transfer. Default True.

• multi_rev (int, optional) – Number of complete revolutions. Default 0 (direct transfer).

• max_iter (int, optional) – Maximum iterations. Default 100.

• tol (float, optional) – Convergence tolerance. Default 1e-10.

Returns:

solution – Solution containing velocities and orbital parameters.

Return type:

LambertSolution

Notes

For multi-revolution transfers, there may be two solutions (low and high energy). This returns the low energy solution.

pytcl.astronomical.minimum_energy_transfer(r1, r2, mu=398600.4418, prograde=True)

Compute minimum energy transfer between two positions.

Parameters:

• r1 (ndarray) – Initial position vector (km).

• r2 (ndarray) – Final position vector (km).

• mu (float, optional) – Gravitational parameter (km^3/s^2).

• prograde (bool, optional) – If True, use prograde transfer.

Returns:

• tof_min (float) – Minimum energy time of flight (seconds).

• solution (LambertSolution) – Lambert solution at minimum energy.

Return type:

Tuple[float, LambertSolution]

pytcl.astronomical.hohmann_transfer(r1, r2, mu=398600.4418)

Compute Hohmann transfer between two circular orbits.

Parameters:

• r1 (float) – Initial orbit radius (km).

• r2 (float) – Final orbit radius (km).

• mu (float, optional) – Gravitational parameter (km^3/s^2).

Returns:

• dv1 (float) – Delta-v at first burn (km/s).

• dv2 (float) – Delta-v at second burn (km/s).

• tof (float) – Transfer time of flight (seconds).

Return type:

Tuple[float, float, float]

Examples

>>> dv1, dv2, tof = hohmann_transfer(6678, 42164)  # LEO to GEO
>>> print(f"Total dv = {dv1 + dv2:.3f} km/s")

pytcl.astronomical.bi_elliptic_transfer(r1, r2, r_intermediate, mu=398600.4418)

Compute bi-elliptic transfer between two circular orbits.

Parameters:

• r1 (float) – Initial orbit radius (km).

• r2 (float) – Final orbit radius (km).

• r_intermediate (float) – Intermediate apoapsis radius (km). Must be > max(r1, r2).

• mu (float, optional) – Gravitational parameter (km^3/s^2).

Returns:

• dv1 (float) – Delta-v at first burn (km/s).

• dv2 (float) – Delta-v at intermediate apoapsis (km/s).

• dv3 (float) – Delta-v at final circularization (km/s).

• tof (float) – Total transfer time (seconds).

Return type:

Tuple[float, float, float, float]

pytcl.astronomical.julian_centuries_j2000(jd)

Compute Julian centuries since J2000.0.

Parameters:

jd (float) – Julian date.

Returns:

T – Julian centuries since J2000.0.

Return type:

float

pytcl.astronomical.precession_angles_iau76(T)

Compute IAU 1976 precession angles.

Parameters:

T (float) – Julian centuries since J2000.0.

Returns:

• zeta (float) – Precession angle zeta (radians).

• theta (float) – Precession angle theta (radians).

• z (float) – Precession angle z (radians).

Return type:

Tuple[float, float, float]

pytcl.astronomical.precession_matrix_iau76(jd)

Compute IAU 1976 precession matrix from J2000 to date.

Parameters:

jd (float) – Julian date of epoch.

Returns:

P – Precession rotation matrix (3x3). Transforms from J2000 (GCRF) to mean of date.

Return type:

ndarray

Notes

Results are cached for repeated queries at the same epoch. Cache key is quantized to ~86ms precision.

pytcl.astronomical.nutation_angles_iau80(jd)

Compute IAU 1980 nutation angles (simplified).

Parameters:

jd (float) – Julian date.

Returns:

• dpsi (float) – Nutation in longitude (radians).

• deps (float) – Nutation in obliquity (radians).

Return type:

Tuple[float, float]

pytcl.astronomical.nutation_matrix(jd)

Compute nutation matrix.

Parameters:

jd (float) – Julian date.

Returns:

N – Nutation rotation matrix (3x3). Transforms from mean of date to true of date.

Return type:

ndarray

Notes

Results are cached for repeated queries at the same epoch. Cache key is quantized to ~86ms precision.

pytcl.astronomical.mean_obliquity_iau80(jd)

Compute mean obliquity of the ecliptic (IAU 1980).

Parameters:

jd (float) – Julian date.

Returns:

eps0 – Mean obliquity (radians).

Return type:

float

pytcl.astronomical.true_obliquity(jd_tt)

Compute true obliquity of the ecliptic.

Parameters:

jd_tt (float) – Julian date in TT.

Returns:

eps – True obliquity (radians).

Return type:

float

pytcl.astronomical.earth_rotation_angle(jd_ut1)

Compute Earth Rotation Angle (ERA).

Parameters:

jd_ut1 (float) – Julian date in UT1.

Returns:

theta – Earth rotation angle (radians), in [0, 2*pi).

Return type:

float

pytcl.astronomical.gmst_iau82(jd_ut1)

Compute Greenwich Mean Sidereal Time (GMST) using IAU 1982 model.

Parameters:

jd_ut1 (float) – Julian date in UT1.

Returns:

gmst – GMST (radians), in [0, 2*pi).

Return type:

float

pytcl.astronomical.gast_iau82(jd_ut1, jd_tt)

Compute Greenwich Apparent Sidereal Time (GAST).

Parameters:

• jd_ut1 (float) – Julian date in UT1.

• jd_tt (float) – Julian date in TT.

Returns:

gast – GAST (radians), in [0, 2*pi).

Return type:

float

pytcl.astronomical.sidereal_rotation_matrix(theta)

Compute Earth rotation matrix (sidereal rotation).

Parameters:

theta (float) – Sidereal angle (GMST or GAST) in radians.

Returns:

R – Rotation matrix (3x3). Transforms from true of date to pseudo Earth-fixed (PEF).

Return type:

ndarray

pytcl.astronomical.equation_of_equinoxes(jd_tt)

Compute equation of the equinoxes.

Parameters:

jd_tt (float) – Julian date in TT.

Returns:

eq_eq – Equation of the equinoxes (radians).

Return type:

float

pytcl.astronomical.polar_motion_matrix(xp, yp)

Compute polar motion rotation matrix.

Parameters:

• xp (float) – Polar motion x coordinate (radians).

• yp (float) – Polar motion y coordinate (radians).

Returns:

W – Polar motion matrix (3x3). Transforms from PEF to ITRF.

Return type:

ndarray

pytcl.astronomical.gcrf_to_itrf(r_gcrf, jd_ut1, jd_tt, xp=0.0, yp=0.0)

Transform position from GCRF (inertial) to ITRF (Earth-fixed).

Parameters:

• r_gcrf (ndarray) – Position in GCRF (km), shape (3,).

• jd_ut1 (float) – Julian date in UT1.

• jd_tt (float) – Julian date in TT.

• xp (float, optional) – Polar motion x (radians). Default 0.

• yp (float, optional) – Polar motion y (radians). Default 0.

Returns:

r_itrf – Position in ITRF (km), shape (3,).

Return type:

ndarray

Examples

>>> r_gcrf = np.array([5102.5096, 6123.01152, 6378.1363])
>>> r_itrf = gcrf_to_itrf(r_gcrf, 2453101.827406783, 2453101.828154745)

pytcl.astronomical.itrf_to_gcrf(r_itrf, jd_ut1, jd_tt, xp=0.0, yp=0.0)

Transform position from ITRF (Earth-fixed) to GCRF (inertial).

Parameters:

• r_itrf (ndarray) – Position in ITRF (km), shape (3,).

• jd_ut1 (float) – Julian date in UT1.

• jd_tt (float) – Julian date in TT.

• xp (float, optional) – Polar motion x (radians). Default 0.

• yp (float, optional) – Polar motion y (radians). Default 0.

Returns:

r_gcrf – Position in GCRF (km), shape (3,).

Return type:

ndarray

pytcl.astronomical.eci_to_ecef(r_eci, gmst)

Simple ECI to ECEF transformation (rotation only).

This is a simplified transformation that only accounts for Earth rotation, ignoring precession, nutation, and polar motion.

Parameters:

• r_eci (ndarray) – Position in ECI (km), shape (3,).

• gmst (float) – Greenwich Mean Sidereal Time (radians).

Returns:

r_ecef – Position in ECEF (km), shape (3,).

Return type:

ndarray

pytcl.astronomical.ecef_to_eci(r_ecef, gmst)

Simple ECEF to ECI transformation (rotation only).

Parameters:

• r_ecef (ndarray) – Position in ECEF (km), shape (3,).

• gmst (float) – Greenwich Mean Sidereal Time (radians).

Returns:

r_eci – Position in ECI (km), shape (3,).

Return type:

ndarray

pytcl.astronomical.ecliptic_to_equatorial(r_ecl, obliquity)

Transform from ecliptic to equatorial coordinates.

Parameters:

• r_ecl (ndarray) – Position in ecliptic frame (km), shape (3,).

• obliquity (float) – Obliquity of the ecliptic (radians).

Returns:

r_eq – Position in equatorial frame (km), shape (3,).

Return type:

ndarray

pytcl.astronomical.equatorial_to_ecliptic(r_eq, obliquity)

Transform from equatorial to ecliptic coordinates.

Parameters:

• r_eq (ndarray) – Position in equatorial frame (km), shape (3,).

• obliquity (float) – Obliquity of the ecliptic (radians).

Returns:

r_ecl – Position in ecliptic frame (km), shape (3,).

Return type:

ndarray

pytcl.astronomical.teme_to_pef(r_teme, jd_ut1)

Transform position from TEME to PEF (Pseudo Earth-Fixed).

TEME is the True Equator, Mean Equinox frame used by SGP4. This transformation applies only the GMST rotation.

Parameters:

• r_teme (ndarray) – Position in TEME frame (km), shape (3,).

• jd_ut1 (float) – Julian date in UT1.

Returns:

r_pef – Position in PEF frame (km), shape (3,).

Return type:

ndarray

pytcl.astronomical.pef_to_teme(r_pef, jd_ut1)

Transform position from PEF to TEME.

Parameters:

• r_pef (ndarray) – Position in PEF frame (km), shape (3,).

• jd_ut1 (float) – Julian date in UT1.

Returns:

r_teme – Position in TEME frame (km), shape (3,).

Return type:

ndarray

pytcl.astronomical.teme_to_itrf(r_teme, jd_ut1, xp=0.0, yp=0.0)

Transform position from TEME to ITRF (Earth-fixed).

TEME is the True Equator, Mean Equinox frame used by SGP4/SDP4. This is the frame in which TLE-propagated positions are expressed.

Parameters:

• r_teme (ndarray) – Position in TEME frame (km), shape (3,).

• jd_ut1 (float) – Julian date in UT1.

• xp (float, optional) – Polar motion x (radians). Default 0.

• yp (float, optional) – Polar motion y (radians). Default 0.

Returns:

r_itrf – Position in ITRF frame (km), shape (3,).

Return type:

ndarray

Notes

TEME is a quasi-inertial frame that uses the mean equinox instead of the true equinox. The transformation sequence is:

TEME -> PEF (via GMST rotation) -> ITRF (via polar motion)

Examples

>>> from pytcl.astronomical.sgp4 import sgp4_propagate
>>> from pytcl.astronomical.tle import parse_tle
>>> tle = parse_tle(line1, line2)
>>> state = sgp4_propagate(tle, 0.0)
>>> r_itrf = teme_to_itrf(state.r, jd_ut1)

pytcl.astronomical.itrf_to_teme(r_itrf, jd_ut1, xp=0.0, yp=0.0)

Transform position from ITRF to TEME.

Parameters:

• r_itrf (ndarray) – Position in ITRF frame (km), shape (3,).

• jd_ut1 (float) – Julian date in UT1.

• xp (float, optional) – Polar motion x (radians). Default 0.

• yp (float, optional) – Polar motion y (radians). Default 0.

Returns:

r_teme – Position in TEME frame (km), shape (3,).

Return type:

ndarray

pytcl.astronomical.teme_to_gcrf(r_teme, jd_tt)

Transform position from TEME to GCRF (inertial).

This transformation accounts for the difference between the mean and true equinox (equation of equinoxes) and then applies precession and nutation to go from TOD to GCRF.

Parameters:

• r_teme (ndarray) – Position in TEME frame (km), shape (3,).

• jd_tt (float) – Julian date in TT (Terrestrial Time).

Returns:

r_gcrf – Position in GCRF frame (km), shape (3,).

Return type:

ndarray

Notes

The transformation sequence is:

TEME -> TOD (via equation of equinoxes) TOD -> MOD (via nutation, inverse) MOD -> GCRF (via precession, inverse)

Examples

>>> state = sgp4_propagate(tle, 60.0)
>>> r_gcrf = teme_to_gcrf(state.r, jd_tt)

pytcl.astronomical.gcrf_to_teme(r_gcrf, jd_tt)

Transform position from GCRF to TEME.

Parameters:

• r_gcrf (ndarray) – Position in GCRF frame (km), shape (3,).

• jd_tt (float) – Julian date in TT.

Returns:

r_teme – Position in TEME frame (km), shape (3,).

Return type:

ndarray

pytcl.astronomical.teme_to_itrf_with_velocity(r_teme, v_teme, jd_ut1, xp=0.0, yp=0.0)

Transform position and velocity from TEME to ITRF.

This properly accounts for the velocity transformation including the Earth’s rotation rate.

Parameters:

• r_teme (ndarray) – Position in TEME frame (km), shape (3,).

• v_teme (ndarray) – Velocity in TEME frame (km/s), shape (3,).

• jd_ut1 (float) – Julian date in UT1.

• xp (float, optional) – Polar motion x (radians). Default 0.

• yp (float, optional) – Polar motion y (radians). Default 0.

Returns:

• r_itrf (ndarray) – Position in ITRF frame (km), shape (3,).

• v_itrf (ndarray) – Velocity in ITRF frame (km/s), shape (3,).

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

pytcl.astronomical.itrf_to_teme_with_velocity(r_itrf, v_itrf, jd_ut1, xp=0.0, yp=0.0)

Transform position and velocity from ITRF to TEME.

Parameters:

• r_itrf (ndarray) – Position in ITRF frame (km), shape (3,).

• v_itrf (ndarray) – Velocity in ITRF frame (km/s), shape (3,).

• jd_ut1 (float) – Julian date in UT1.

• xp (float, optional) – Polar motion x (radians). Default 0.

• yp (float, optional) – Polar motion y (radians). Default 0.

Returns:

• r_teme (ndarray) – Position in TEME frame (km), shape (3,).

• v_teme (ndarray) – Velocity in TEME frame (km/s), shape (3,).

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

pytcl.astronomical.gcrf_to_mod(r_gcrf, jd_tt)

Transform position from GCRF to MOD (Mean of Date).

MOD is the mean equator and mean equinox of date frame. This applies only the precession transformation.

Parameters:

• r_gcrf (ndarray) – Position in GCRF frame (km), shape (3,).

• jd_tt (float) – Julian date in TT (Terrestrial Time).

Returns:

r_mod – Position in MOD frame (km), shape (3,).

Return type:

ndarray

Notes

MOD is a legacy frame convention. For most modern applications, GCRF (J2000) is preferred. MOD was historically used in older software and publications.

The transformation is simply the precession matrix:

r_mod = P @ r_gcrf

See also

mod_to_gcrf

Inverse transformation.

gcrf_to_tod

Includes nutation for true of date.

pytcl.astronomical.mod_to_gcrf(r_mod, jd_tt)

Transform position from MOD (Mean of Date) to GCRF.

Parameters:

• r_mod (ndarray) – Position in MOD frame (km), shape (3,).

• jd_tt (float) – Julian date in TT.

Returns:

r_gcrf – Position in GCRF frame (km), shape (3,).

Return type:

ndarray

See also

gcrf_to_mod

Forward transformation.

pytcl.astronomical.gcrf_to_tod(r_gcrf, jd_tt)

Transform position from GCRF to TOD (True of Date).

TOD is the true equator and true equinox of date frame. This applies both precession and nutation transformations.

Parameters:

• r_gcrf (ndarray) – Position in GCRF frame (km), shape (3,).

• jd_tt (float) – Julian date in TT (Terrestrial Time).

Returns:

r_tod – Position in TOD frame (km), shape (3,).

Return type:

ndarray

Notes

TOD is a legacy frame convention. The transformation is:

r_mod = P @ r_gcrf r_tod = N @ r_mod

where P is the precession matrix and N is the nutation matrix.

See also

tod_to_gcrf

Inverse transformation.

gcrf_to_mod

Mean of date (without nutation).

pytcl.astronomical.tod_to_gcrf(r_tod, jd_tt)

Transform position from TOD (True of Date) to GCRF.

Parameters:

• r_tod (ndarray) – Position in TOD frame (km), shape (3,).

• jd_tt (float) – Julian date in TT.

Returns:

r_gcrf – Position in GCRF frame (km), shape (3,).

Return type:

ndarray

See also

gcrf_to_tod

Forward transformation.

pytcl.astronomical.mod_to_tod(r_mod, jd_tt)

Transform position from MOD (Mean of Date) to TOD (True of Date).

This applies only the nutation transformation.

Parameters:

• r_mod (ndarray) – Position in MOD frame (km), shape (3,).

• jd_tt (float) – Julian date in TT.

Returns:

r_tod – Position in TOD frame (km), shape (3,).

Return type:

ndarray

Notes

The transformation is simply the nutation matrix:

r_tod = N @ r_mod

See also

tod_to_mod

Inverse transformation.

pytcl.astronomical.tod_to_mod(r_tod, jd_tt)

Transform position from TOD (True of Date) to MOD (Mean of Date).

Parameters:

• r_tod (ndarray) – Position in TOD frame (km), shape (3,).

• jd_tt (float) – Julian date in TT.

Returns:

r_mod – Position in MOD frame (km), shape (3,).

Return type:

ndarray

See also

mod_to_tod

Forward transformation.

pytcl.astronomical.tod_to_itrf(r_tod, jd_ut1, jd_tt=None, xp=0.0, yp=0.0)

Transform position from TOD (True of Date) to ITRF.

Parameters:

• r_tod (ndarray) – Position in TOD frame (km), shape (3,).

• jd_ut1 (float) – Julian date in UT1.

• jd_tt (float, optional) – Julian date in TT. If not provided, assumed equal to jd_ut1.

• xp (float, optional) – Polar motion x (radians). Default 0.

• yp (float, optional) – Polar motion y (radians). Default 0.

Returns:

r_itrf – Position in ITRF frame (km), shape (3,).

Return type:

ndarray

Notes

The transformation applies the sidereal rotation (using GAST) and polar motion:

r_pef = R(GAST) @ r_tod r_itrf = W @ r_pef

See also

itrf_to_tod

Inverse transformation.

pytcl.astronomical.itrf_to_tod(r_itrf, jd_ut1, jd_tt=None, xp=0.0, yp=0.0)

Transform position from ITRF to TOD (True of Date).

Parameters:

• r_itrf (ndarray) – Position in ITRF frame (km), shape (3,).

• jd_ut1 (float) – Julian date in UT1.

• jd_tt (float, optional) – Julian date in TT. If not provided, assumed equal to jd_ut1.

• xp (float, optional) – Polar motion x (radians). Default 0.

• yp (float, optional) – Polar motion y (radians). Default 0.

Returns:

r_tod – Position in TOD frame (km), shape (3,).

Return type:

ndarray

See also

tod_to_itrf

Forward transformation.

class pytcl.astronomical.TLE(name, catalog_number, classification, int_designator, epoch_year, epoch_day, ndot, nddot, bstar, ephemeris_type, element_set_number, inclination, raan, eccentricity, arg_perigee, mean_anomaly, mean_motion, revolution_number, line1, line2)

Bases: NamedTuple

Two-Line Element set data structure.

Contains all orbital elements and metadata from a NORAD TLE.

name

Satellite name (from line 0, if present).

Type:

str

catalog_number

NORAD catalog number.

Type:

int

classification

Classification (‘U’ = unclassified).

Type:

str

int_designator

International designator (e.g., ‘98067A’).

Type:

str

epoch_year

Epoch year (4-digit).

Type:

int

epoch_day

Epoch day of year (fractional).

Type:

float

ndot

First derivative of mean motion (rev/day^2).

Type:

float

nddot

Second derivative of mean motion (rev/day^3).

Type:

float

bstar

BSTAR drag coefficient (1/Earth radii).

Type:

float

ephemeris_type

Ephemeris type (usually 0 for SGP4).

Type:

int

element_set_number

Element set number.

Type:

int

inclination

Inclination (radians).

Type:

float

raan

Right ascension of ascending node (radians).

Type:

float

eccentricity

Eccentricity (dimensionless).

Type:

float

arg_perigee

Argument of perigee (radians).

Type:

float

mean_anomaly

Mean anomaly (radians).

Type:

float

mean_motion

Mean motion (radians/minute).

Type:

float

revolution_number

Revolution number at epoch.

Type:

int

line1

Original TLE line 1.

Type:

str

line2

Original TLE line 2.

Type:

str

Notes

Angular quantities are stored in radians for consistency with the rest of the pytcl library. Mean motion is in radians/minute as required by SGP4/SDP4.

name: str

Alias for field number 0

catalog_number: int

Alias for field number 1

classification: str

Alias for field number 2

int_designator: str

Alias for field number 3

epoch_year: int

Alias for field number 4

epoch_day: float

Alias for field number 5

ndot: float

Alias for field number 6

nddot: float

Alias for field number 7

bstar: float

Alias for field number 8

ephemeris_type: int

Alias for field number 9

element_set_number: int

Alias for field number 10

inclination: float

Alias for field number 11

raan: float

Alias for field number 12

eccentricity: float

Alias for field number 13

arg_perigee: float

Alias for field number 14

mean_anomaly: float

Alias for field number 15

mean_motion: float

Alias for field number 16

revolution_number: int

Alias for field number 17

line1: str

Alias for field number 18

line2: str

Alias for field number 19

pytcl.astronomical.parse_tle(line1, line2, name='', verify_checksum=True)

Parse a Two-Line Element set.

Parameters:

• line1 (str) – First line of the TLE (69 characters).

• line2 (str) – Second line of the TLE (69 characters).

• name (str, optional) – Satellite name. Default empty.

• verify_checksum (bool, optional) – Whether to verify line checksums. Default True.

Returns:

tle – Parsed TLE data structure.

Return type:

TLE

Raises:

ValueError – If TLE format is invalid or checksum fails.

Examples

>>> line1 = "1 25544U 98067A   24001.50000000  .00016717  00000-0  10270-3 0  9993"
>>> line2 = "2 25544  51.6400 247.4627 0006703 130.5360 325.0288 15.49815350479001"
>>> tle = parse_tle(line1, line2, name="ISS (ZARYA)")
>>> print(f"Inclination: {np.degrees(tle.inclination):.4f} deg")

pytcl.astronomical.parse_tle_3line(lines)

Parse a three-line TLE (with satellite name).

Parameters:

lines (str) – Three-line TLE string (name, line1, line2).

Returns:

tle – Parsed TLE data structure.

Return type:

TLE

Examples

>>> tle_text = '''ISS (ZARYA)
... 1 25544U 98067A   24001.50000000  .00016717  00000-0  10270-3 0  9993
... 2 25544  51.6400 247.4627 0006703 130.5360 325.0288 15.49815350479001'''
>>> tle = parse_tle_3line(tle_text)
>>> print(tle.name)
ISS (ZARYA)

pytcl.astronomical.tle_epoch_to_jd(tle)

Convert TLE epoch to Julian date.

Parameters:

tle (TLE) – Parsed TLE.

Returns:

jd – Julian date of TLE epoch.

Return type:

float

Examples

>>> tle = parse_tle(line1, line2)
>>> jd = tle_epoch_to_jd(tle)

pytcl.astronomical.tle_epoch_to_datetime(tle)

Convert TLE epoch to Python datetime.

Parameters:

tle (TLE) – Parsed TLE.

Returns:

dt – TLE epoch as UTC datetime.

Return type:

datetime

pytcl.astronomical.format_tle(tle, include_name=True)

Format TLE data structure back to TLE string.

Parameters:

• tle (TLE) – TLE data structure.

• include_name (bool, optional) – Whether to include satellite name as line 0. Default True.

Returns:

Formatted TLE string.

Return type:

str

pytcl.astronomical.is_deep_space(tle)

Determine if TLE requires deep-space (SDP4) propagation.

Satellites with orbital period >= 225 minutes use SDP4 instead of SGP4 due to lunar-solar perturbations.

Parameters:

tle (TLE) – Parsed TLE.

Returns:

True if deep-space propagation (SDP4) is required.

Return type:

bool

pytcl.astronomical.semi_major_axis_from_mean_motion(n, mu=398600.4418)

Compute semi-major axis from mean motion.

Uses the relationship n = sqrt(mu / a^3) where n is in rad/s.

Parameters:

• n (float) – Mean motion (radians/minute).

• mu (float, optional) – Gravitational parameter (km^3/s^2). Default is Earth.

Returns:

a – Semi-major axis (km).

Return type:

float

pytcl.astronomical.orbital_period_from_tle(tle)

Compute orbital period from TLE mean motion.

Parameters:

tle (TLE) – Parsed TLE.

Returns:

period – Orbital period (seconds).

Return type:

float

class pytcl.astronomical.SGP4State(r, v, error)

Bases: NamedTuple

State vector from SGP4 propagation.

r

Position in TEME frame (km), shape (3,).

Type:

ndarray

v

Velocity in TEME frame (km/s), shape (3,).

Type:

ndarray

error

Error code (0 = success).

Type:

int

r: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 0

v: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 1

error: int

Alias for field number 2

class pytcl.astronomical.SGP4Satellite(tle)

Bases: object

SGP4 satellite propagator initialized from a TLE.

This class encapsulates the initialization and propagation logic for a satellite using the SGP4/SDP4 models.

Parameters:

tle (TLE) – Two-Line Element set.

tle

Original TLE data.

Type:

TLE

epoch_jd

Julian date of TLE epoch.

Type:

float

is_deep_space

True if SDP4 (deep-space) propagation is used.

Type:

bool

Examples

>>> tle = parse_tle(line1, line2, name="ISS")
>>> sat = SGP4Satellite(tle)
>>> state = sat.propagate(0.0)  # At epoch
>>> print(f"Position: {state.r} km")
>>> state = sat.propagate(60.0)  # 60 minutes later

__init__(tle)

Initialize SGP4 satellite from TLE.

propagate(tsince)

Propagate satellite to specified time.

Parameters:

tsince (float) – Time since epoch (minutes). Positive = after epoch.

Returns:

state – Position and velocity in TEME frame.

Return type:

SGP4State

Examples

>>> sat = SGP4Satellite(tle)
>>> state = sat.propagate(0.0)  # At TLE epoch
>>> state = sat.propagate(60.0)  # 60 minutes later
>>> state = sat.propagate(-30.0)  # 30 minutes before epoch

propagate_jd(jd)

Propagate satellite to specified Julian date.

Parameters:

jd (float) – Julian date.

Returns:

state – Position and velocity in TEME frame.

Return type:

SGP4State

pytcl.astronomical.sgp4_propagate(tle, tsince)

Propagate TLE using SGP4/SDP4 model.

Convenience function that creates an SGP4Satellite and propagates.

Parameters:

• tle (TLE) – Two-Line Element set.

• tsince (float) – Time since epoch (minutes).

Returns:

state – Position and velocity in TEME frame.

Return type:

SGP4State

Examples

>>> tle = parse_tle(line1, line2)
>>> state = sgp4_propagate(tle, 60.0)  # 60 minutes after epoch
>>> print(f"Position: {state.r} km")

pytcl.astronomical.sgp4_propagate_batch(tle, times)

Propagate TLE to multiple times.

Parameters:

• tle (TLE) – Two-Line Element set.

• times (ndarray) – Times since epoch (minutes), shape (n,).

Returns:

• positions (ndarray) – Positions in TEME frame (km), shape (n, 3).

• velocities (ndarray) – Velocities in TEME frame (km/s), shape (n, 3).

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Examples

>>> tle = parse_tle(line1, line2)
>>> times = np.linspace(0, 90, 100)  # 0 to 90 minutes
>>> r, v = sgp4_propagate_batch(tle, times)

class pytcl.astronomical.DEEphemeris(version='DE440')

Bases: object

High-precision JPL Development Ephemeris kernel wrapper.

This class manages access to JPL ephemeris files and provides methods for querying positions and velocities of celestial bodies.

Parameters:

version ({'DE405', 'DE430', 'DE432s', 'DE440'}, optional) – Ephemeris version to load. Default is ‘DE440’ (latest). - DE440: Latest JPL release (2020), covers 1550-2650 - DE432s: High-precision version (2013), covers 1350-3000 - DE430: Earlier release (2013), covers 1550-2650 - DE405: Older version (1998), compact, covers 1600-2200

version

Ephemeris version identifier

Type:

str

kernel

Loaded ephemeris kernel object

Type:

jplephem.SpiceKernel

_cache

Cache for frequently accessed positions

Type:

dict

Raises:

• DependencyError – If jplephem is not installed

• ValueError – If version is not recognized

Examples

>>> eph = DEEphemeris(version='DE440')
>>> r_sun, v_sun = eph.sun_position(2451545.0)

__init__(version='DE440')

Initialize ephemeris kernel.

Parameters:

version (str, optional) – Ephemeris version (default: ‘DE440’)

property kernel: object | None

Lazy-load ephemeris kernel on first access.

Note: This requires jplephem to be installed and the kernel file to be available locally or downloadable from the JPL servers.

sun_position(jd, frame='icrf')

Compute Sun position and velocity.

Parameters:

• jd (float) – Julian Day in Terrestrial Time (TT)

• frame ({'icrf', 'ecliptic'}, optional) – Coordinate frame (default: ‘icrf’). - ‘icrf’: International Celestial Reference Frame - ‘ecliptic’: Ecliptic coordinate system (J2000.0)

Returns:

• position (ndarray, shape (3,)) – Sun position in AU

• velocity (ndarray, shape (3,)) – Sun velocity in AU/day

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Notes

The Sun’s position is computed relative to the Solar System Barycenter (SSB) in the ICRF frame.

Examples

>>> eph = DEEphemeris()
>>> r, v = eph.sun_position(2451545.0)
>>> print(f"Distance: {np.linalg.norm(r):.6f} AU")

moon_position(jd, frame='icrf')

Compute Moon position and velocity.

Parameters:

• jd (float) – Julian Day in Terrestrial Time (TT)

• frame ({'icrf', 'ecliptic', 'earth_centered'}, optional) – Coordinate frame (default: ‘icrf’). - ‘icrf’: Moon position relative to Solar System Barycenter - ‘ecliptic’: Ecliptic coordinates - ‘earth_centered’: Position relative to Earth

Returns:

• position (ndarray, shape (3,)) – Moon position in AU (or relative to Earth for ‘earth_centered’)

• velocity (ndarray, shape (3,)) – Moon velocity in AU/day

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Notes

By default, returns Moon position relative to the Solar System Barycenter. Use frame=’earth_centered’ for geocentric coordinates.

Examples

>>> eph = DEEphemeris()
>>> r, v = eph.moon_position(2451545.0, frame='earth_centered')

planet_position(planet, jd, frame='icrf')

Compute planet position and velocity.

Parameters:

• planet (str) – Planet name: ‘mercury’, ‘venus’, ‘mars’, ‘jupiter’, ‘saturn’, ‘uranus’, ‘neptune’

• jd (float) – Julian Day in Terrestrial Time (TT)

• frame ({'icrf', 'ecliptic'}, optional) – Coordinate frame (default: ‘icrf’)

Returns:

• position (ndarray, shape (3,)) – Planet position in AU

• velocity (ndarray, shape (3,)) – Planet velocity in AU/day

Raises:

ValueError – If planet name is not recognized

Return type:

Tuple[ndarray[Any, Any], ndarray[Any, Any]]

Examples

>>> eph = DEEphemeris()
>>> r, v = eph.planet_position('mars', 2451545.0)

barycenter_position(body, jd)

Compute position of any body relative to Solar System Barycenter.

Parameters:

• body (str) – Body name (‘sun’, ‘moon’, ‘mercury’, …, ‘neptune’)

• jd (float) – Julian Day in Terrestrial Time (TT)

Returns:

• position (ndarray, shape (3,)) – Position in AU

• velocity (ndarray, shape (3,)) – Velocity in AU/day

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

clear_cache()

Clear internal position cache.

pytcl.astronomical.sun_position(jd, frame='icrf')

Convenience function: Compute Sun position and velocity.

Parameters:

• jd (float) – Julian Day in Terrestrial Time (TT)

• frame ({'icrf', 'ecliptic'}, optional) – Coordinate frame (default: ‘icrf’)

Returns:

• position (ndarray, shape (3,)) – Sun position in AU

• velocity (ndarray, shape (3,)) – Sun velocity in AU/day

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Examples

>>> r, v = sun_position(2451545.0)  # J2000.0
>>> np.linalg.norm(r)  # Distance from SSB
0.00...

See also

DEEphemeris.sun_position

Full ephemeris class with caching

pytcl.astronomical.moon_position(jd, frame='icrf')

Convenience function: Compute Moon position and velocity.

Parameters:

• jd (float) – Julian Day in Terrestrial Time (TT)

• frame ({'icrf', 'ecliptic', 'earth_centered'}, optional) – Coordinate frame (default: ‘icrf’)

Returns:

• position (ndarray, shape (3,)) – Moon position in AU

• velocity (ndarray, shape (3,)) – Moon velocity in AU/day

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Examples

>>> r, v = moon_position(2451545.0, 'earth_centered')
>>> np.linalg.norm(r) * 149597870.7  # Distance in km
402...

See also

DEEphemeris.moon_position

Full ephemeris class with caching

pytcl.astronomical.planet_position(planet, jd, frame='icrf')

Convenience function: Compute planet position and velocity.

Parameters:

• planet (str) – Planet name

• jd (float) – Julian Day in Terrestrial Time (TT)

• frame ({'icrf', 'ecliptic'}, optional) – Coordinate frame (default: ‘icrf’)

Returns:

• position (ndarray, shape (3,)) – Planet position in AU

• velocity (ndarray, shape (3,)) – Planet velocity in AU/day

Return type:

Tuple[ndarray[Any, Any], ndarray[Any, Any]]

See also

DEEphemeris.planet_position

Full ephemeris class with caching

pytcl.astronomical.barycenter_position(body, jd)

Convenience function: Position relative to Solar System Barycenter.

Parameters:

• body (str) – Body name (‘sun’, ‘moon’, ‘mercury’, …, ‘neptune’)

• jd (float) – Julian Day in Terrestrial Time (TT)

Returns:

• position (ndarray, shape (3,)) – Position in AU

• velocity (ndarray, shape (3,)) – Velocity in AU/day

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Examples

>>> r, v = barycenter_position('mars', 2451545.0)
>>> np.linalg.norm(r)  # Mars distance from SSB
1.4...

pytcl.astronomical.schwarzschild_radius(mass)

Compute Schwarzschild radius for a given mass.

The Schwarzschild radius is the radius at which an object becomes a black hole. It is given by r_s = 2GM/c^2.

Parameters:

mass (float) – Mass of the object (kg)

Returns:

Schwarzschild radius (m)

Return type:

float

Examples

>>> r_s_earth = schwarzschild_radius(5.972e24)  # Earth's mass
>>> print(f"Earth's Schwarzschild radius: {r_s_earth:.3e} m")
Earth's Schwarzschild radius: 8.870e-03 m

>>> r_s_sun = schwarzschild_radius(1.989e30)  # Sun's mass
>>> print(f"Sun's Schwarzschild radius: {r_s_sun:.3e} m")
Sun's Schwarzschild radius: 2.952e+03 m

pytcl.astronomical.gravitational_time_dilation(r, gm=398600441800000.0)

Compute gravitational time dilation factor sqrt(1 - 2GM/(rc^2)).

In general relativity, time passes slower in stronger gravitational fields. This function computes the metric coefficient g_00 for the Schwarzschild metric, which determines proper time relative to coordinate time.

Parameters:

• r (float) – Distance from the gravitational body (m)

• gm (float, optional) – Gravitational parameter GM of the body (m^3/s^2). Default is GM_EARTH.

Returns:

Time dilation factor in [0, 1]. A value less than 1 indicates time passes slower at radius r compared to infinity.

Return type:

float

Raises:

ValueError – If r is less than or equal to Schwarzschild radius

Examples

Compute time dilation at Earth’s surface (6371 km):

>>> r_earth = 6.371e6  # meters
>>> dilation = gravitational_time_dilation(r_earth)
>>> print(f"Time dilation at surface: {dilation:.15f}")
Time dilation at surface: 0.999999999300693

At GPS orbital altitude (~20,200 km):

>>> r_gps = 26.56e6  # meters
>>> dilation_gps = gravitational_time_dilation(r_gps)
>>> time_shift = (1 - dilation_gps) * 86400 * 1e9  # nanoseconds per day
>>> print(f"Time shift: {time_shift:.1f} ns/day")

pytcl.astronomical.proper_time_rate(v, r, gm=398600441800000.0)

Compute proper time rate accounting for both velocity and gravity.

The proper time rate combines special relativistic time dilation from velocity and general relativistic time dilation from the gravitational potential.

d(tau)/d(t) = sqrt(1 - v^2/c^2) * sqrt(1 - 2GM/(rc^2))

For small velocities and weak fields: 1 - v^2/(2c^2) - GM/(rc^2)

Parameters:

• v (float) – Velocity magnitude (m/s)

• r (float) – Distance from gravitational body (m)

• gm (float, optional) – Gravitational parameter GM (m^3/s^2). Default is GM_EARTH.

Returns:

Proper time rate. A value less than 1 indicates proper time passes slower than coordinate time.

Return type:

float

Examples

Proper time rate for a GPS satellite at ~3.87 km/s and 26.56 Mm altitude:

>>> v_gps = 3870.0  # m/s
>>> r_gps = 26.56e6  # m
>>> rate = proper_time_rate(v_gps, r_gps)
>>> print(f"Proper time rate: {rate:.15f}")
>>> time_shift = (1 - rate) * 86400  # seconds per day
>>> print(f"Daily time shift: {time_shift:.3f} s/day")

pytcl.astronomical.shapiro_delay(observer_pos, light_source_pos, gravitating_body_pos, gm=1.32712440018e+20)

Compute Shapiro time delay for light propagation through gravitational field.

The Shapiro delay is the additional propagation time experienced by light traveling through a gravitational field, compared to flat spacetime.

delay = (2GM/c^3) * ln((r_o + r_s + r_os) / (r_o + r_s - r_os))

where r_o is distance from body to observer, r_s is distance from body to source, and r_os is distance from observer to source.

Parameters:

• observer_pos (np.ndarray) – Position of observer (m), shape (3,)

• light_source_pos (np.ndarray) – Position of light source (m), shape (3,)

• gravitating_body_pos (np.ndarray) – Position of gravitating body (m), shape (3,)

• gm (float, optional) – Gravitational parameter GM (m^3/s^2). Default is GM_SUN.

Returns:

Shapiro delay (seconds)

Return type:

float

Examples

Earth-Sun-Spacecraft signal at superior conjunction (worst case):

>>> # Simplified geometry: Sun at origin, Earth at 1 AU, spacecraft beyond at distance
>>> sun_pos = np.array([0.0, 0.0, 0.0])
>>> earth_pos = np.array([1.496e11, 0.0, 0.0])  # 1 AU
>>> spacecraft_pos = np.array([1.496e11, 1.0e11, 0.0])  # Far from sun
>>> delay = shapiro_delay(earth_pos, spacecraft_pos, sun_pos, GM_SUN)
>>> print(f"Shapiro delay: {delay:.3e} seconds")
>>> print(f"Shapiro delay: {delay*1e6:.1f} microseconds")

pytcl.astronomical.schwarzschild_precession_per_orbit(a, e, gm=1.32712440018e+20)

Compute perihelion precession per orbit due to general relativity.

The advance of perihelion for an orbit around a central mass M is:

Δφ = (6π * GM) / (c^2 * a * (1 - e^2))

This effect is a key test of general relativity. For Mercury, the predicted precession is ~43 arcseconds per century.

Parameters:

• a (float) – Semi-major axis (m)

• e (float) – Eccentricity (dimensionless), must be in [0, 1)

• gm (float, optional) – Gravitational parameter GM (m^3/s^2). Default is GM_SUN.

Returns:

Perihelion precession per orbit (radians)

Return type:

float

Examples

Mercury’s perihelion precession (GR contribution):

>>> a_mercury = 0.38709927 * AU  # Semi-major axis in meters
>>> e_mercury = 0.20563593     # Eccentricity
>>> precession_rad = schwarzschild_precession_per_orbit(a_mercury, e_mercury, GM_SUN)
>>> precession_arcsec = precession_rad * 206265  # Convert to arcseconds
>>> orbital_period = 87.969  # days
>>> centuries = 36525 / orbital_period  # Orbits per century
>>> precession_per_century = precession_arcsec * centuries
>>> print(f"GR perihelion precession: {precession_per_century:.1f} arcsec/century")
GR perihelion precession: 42.98 arcsec/century

pytcl.astronomical.post_newtonian_acceleration(r_vec, v_vec, gm=398600441800000.0)

Compute post-Newtonian acceleration corrections (1PN order).

Extends Newtonian gravity with first-order post-Newtonian corrections.

a_PN = -GM/r^2 * u_r + a_1PN

where a_1PN includes velocity-dependent and metric perturbation terms.

Parameters:

• r_vec (np.ndarray) – Position vector (m), shape (3,)

• v_vec (np.ndarray) – Velocity vector (m/s), shape (3,)

• gm (float, optional) – Gravitational parameter GM (m^3/s^2). Default is GM_EARTH.

Returns:

Total acceleration including 1PN corrections (m/s^2), shape (3,)

Return type:

np.ndarray

Examples

Compare Newtonian and PN acceleration for LEO satellite:

>>> r = np.array([6.678e6, 0.0, 0.0])  # ~300 km altitude
>>> v = np.array([0.0, 7.7e3, 0.0])    # Circular orbit velocity
>>> a_total = post_newtonian_acceleration(r, v)
>>> a_newt = -GM_EARTH / np.linalg.norm(r)**3 * r
>>> correction_ratio = np.linalg.norm(a_total - a_newt) / np.linalg.norm(a_newt)
>>> print(f"PN correction: {correction_ratio*1e6:.1f} ppm")

pytcl.astronomical.geodetic_precession(a, e, inclination, gm=398600441800000.0)

Compute geodetic (de Sitter) precession rate of orbital plane.

The orbital plane of a satellite precesses due to frame-dragging effects and spacetime curvature. The geodetic precession rate is:

Ω_geodetic = -GM/(c^2 * a^3 * (1 - e^2)^2) * cos(i)

Parameters:

• a (float) – Semi-major axis (m)

• e (float) – Eccentricity (dimensionless)

• inclination (float) – Orbital inclination (radians)

• gm (float, optional) – Gravitational parameter (m^3/s^2). Default is GM_EARTH.

Returns:

Geodetic precession rate (radians per orbit)

Return type:

float

Examples

Geodetic precession for a typical Earth satellite:

>>> a = 6.678e6  # ~300 km altitude
>>> e = 0.0      # Circular
>>> i = np.radians(51.6)  # ISS-like inclination
>>> rate = geodetic_precession(a, e, i)
>>> print(f"Precession per orbit: {rate*206265:.3f} arcsec")

pytcl.astronomical.lense_thirring_precession(a, e, inclination, angular_momentum, gm=398600441800000.0)

Compute Lense-Thirring (frame-dragging) precession of orbital node.

A rotating central body drags the orbital plane of nearby objects. The nodal precession rate due to this effect is:

Ω_LT = (2GM * J_2 * a * ω) / (c^2 * p^2) * f(e, i)

where J_2 is the quadrupole moment, ω is angular velocity, and f depends on eccentricity and inclination.

Parameters:

• a (float) – Semi-major axis (m)

• e (float) – Eccentricity (dimensionless)

• inclination (float) – Orbital inclination (radians)

• angular_momentum (float) – Angular momentum of central body (kg·m^2/s)

• gm (float, optional) – Gravitational parameter (m^3/s^2). Default is GM_EARTH.

Returns:

Lense-Thirring precession rate (radians per orbit)

Return type:

float

Notes

This is a simplified version. For Earth, J_2 effects typically dominate classical nodal precession, while Lense-Thirring is a small correction (~1 arcsec per year for typical satellites).

Examples

Lense-Thirring effect for LAGEOS satellite:

>>> # LAGEOS parameters
>>> a = 12.27e6  # Semi-major axis
>>> e = 0.0045
>>> i = np.radians(109.9)
>>> L_earth = 7.05e33  # Earth's angular momentum
>>> rate = lense_thirring_precession(a, e, i, L_earth)
>>> print(f"LT precession per orbit: {rate*206265*1e3:.1f} milliarcsec")

pytcl.astronomical.relativistic_range_correction(distance, relative_velocity, gm=398600441800000.0)

Compute relativistic range correction for ranging measurements.

When measuring distance to a satellite or spacecraft using ranging (e.g., laser ranging), relativistic effects introduce corrections to the measured range.

The main contributions are: - Gravitational time dilation - Relativistic Doppler effect

Parameters:

• distance (float) – Distance to object (m)

• relative_velocity (float) – Radial velocity component (m/s, positive = receding)

• gm (float, optional) – Gravitational parameter (m^3/s^2). Default is GM_EARTH.

Returns:

Range correction (m)

Return type:

float

Examples

Range correction for lunar laser ranging:

>>> distance_to_moon = 3.84e8  # meters
>>> radial_velocity = 0.0  # Average over orbit
>>> correction = relativistic_range_correction(distance_to_moon, radial_velocity, GM_EARTH)
>>> print(f"Range correction: {correction:.1f} m")

class pytcl.astronomical.OrbitType(value)

Bases: Enum

Classification of orbit types based on eccentricity.

CIRCULAR = 0

ELLIPTICAL = 1

PARABOLIC = 2

HYPERBOLIC = 3

pytcl.astronomical.classify_orbit(e, tol=1e-09)

Classify orbit type based on eccentricity.

Parameters:

• e (float) – Eccentricity value.

• tol (float, optional) – Tolerance for parabolic classification (default 1e-9). Orbits with abs(e - 1) < tol are classified as parabolic.

Returns:

Classified orbit type.

Return type:

OrbitType

Raises:

ValueError – If eccentricity is negative or NaN.

pytcl.astronomical.mean_to_parabolic_anomaly(M, tol=1e-12, max_iter=100)

Solve parabolic Kepler’s equation: M = D + (1/3)*D^3.

For parabolic orbits (e=1), the “anomaly” is the parabolic anomaly D, related to true anomaly by: D = tan(nu/2).

The equation relates mean anomaly to parabolic anomaly: M = D + (1/3)*D^3

This is solved numerically using Newton-Raphson iteration.

Parameters:

• M (float) – Mean anomaly (radians).

• tol (float, optional) – Convergence tolerance (default 1e-12).

• max_iter (int, optional) – Maximum iterations (default 100).

Returns:

D – Parabolic anomaly (the parameter D such that tan(nu/2) = D).

Return type:

float

Notes

For parabolic orbits, mean anomaly relates to time as: M = sqrt(mu/rp^3) * t where rp is periapsis distance and mu is GM.

The solution D satisfies: D + (1/3)*D^3 = M

pytcl.astronomical.parabolic_anomaly_to_true_anomaly(D)

Convert parabolic anomaly to true anomaly.

For parabolic orbits, the parabolic anomaly D relates to true anomaly by: tan(nu/2) = D

Parameters:

D (float) – Parabolic anomaly (the parameter such that tan(nu/2) = D).

Returns:

nu – True anomaly (radians), in [-pi, pi].

Return type:

float

pytcl.astronomical.true_anomaly_to_parabolic_anomaly(nu)

Convert true anomaly to parabolic anomaly.

Parameters:

nu (float) – True anomaly (radians).

Returns:

D – Parabolic anomaly.

Return type:

float

pytcl.astronomical.mean_to_true_anomaly_parabolic(M, tol=1e-12)

Direct conversion from mean to true anomaly for parabolic orbits.

Parameters:

• M (float) – Mean anomaly (radians).

• tol (float, optional) – Convergence tolerance (default 1e-12).

Returns:

nu – True anomaly (radians).

Return type:

float

pytcl.astronomical.radius_parabolic(rp, nu)

Compute radius for parabolic orbit.

For a parabolic orbit with periapsis distance rp and true anomaly nu: r = 2*rp / (1 + cos(nu))

This formula is consistent with the general conic section equation with e=1: r = p/(1 + e*cos(nu)) where p = 2*rp (semi-latus rectum).

Parameters:

• rp (float) – Periapsis distance (km).

• nu (float) – True anomaly (radians).

Returns:

r – Orbital radius (km).

Return type:

float

Raises:

ValueError – If radius would be negative (nu near +pi for parabolic orbit).

pytcl.astronomical.velocity_parabolic(mu, rp, nu)

Compute velocity magnitude for parabolic orbit.

For a parabolic orbit (e=1), the specific orbital energy is zero, and the velocity relates to radius by: v = sqrt(2*mu/r)

Parameters:

• mu (float) – Standard gravitational parameter (km^3/s^2).

• rp (float) – Periapsis distance (km).

• nu (float) – True anomaly (radians).

Returns:

v – Velocity magnitude (km/s).

Return type:

float

pytcl.astronomical.hyperbolic_anomaly_to_true_anomaly(H, e)

Convert hyperbolic anomaly to true anomaly.

For hyperbolic orbits (e > 1), hyperbolic anomaly H relates to true anomaly by: tan(nu/2) = sqrt((e+1)/(e-1)) * tanh(H/2)

Parameters:

• H (float) – Hyperbolic anomaly (radians).

• e (float) – Eccentricity (e > 1 for hyperbolic).

Returns:

nu – True anomaly (radians).

Return type:

float

Raises:

ValueError – If eccentricity is not hyperbolic (e <= 1).

pytcl.astronomical.true_anomaly_to_hyperbolic_anomaly(nu, e)

Convert true anomaly to hyperbolic anomaly.

Parameters:

• nu (float) – True anomaly (radians).

• e (float) – Eccentricity (e > 1 for hyperbolic).

Returns:

H – Hyperbolic anomaly (radians).

Return type:

float

Raises:

ValueError – If eccentricity is not hyperbolic.

pytcl.astronomical.escape_velocity_at_radius(mu, r)

Compute escape velocity at a given radius.

Escape velocity is the minimum velocity needed to reach infinity with zero velocity, corresponding to a parabolic orbit: v_esc = sqrt(2*mu/r)

Parameters:

• mu (float) – Standard gravitational parameter (km^3/s^2).

• r (float) – Orbital radius (km).

Returns:

v_esc – Escape velocity (km/s).

Return type:

float

pytcl.astronomical.hyperbolic_excess_velocity(mu, a)

Compute hyperbolic excess velocity.

For a hyperbolic orbit with semi-major axis a (negative for hyperbolic), the excess velocity at infinity is: v_inf = sqrt(-mu/a)

Parameters:

• mu (float) – Standard gravitational parameter (km^3/s^2).

• a (float) – Semi-major axis (km). Must be negative for hyperbolic orbits.

Returns:

v_inf – Hyperbolic excess velocity (km/s).

Return type:

float

Raises:

ValueError – If semi-major axis is not negative.

pytcl.astronomical.semi_major_axis_from_energy(mu, specific_energy)

Compute semi-major axis from specific orbital energy.

The specific orbital energy relates to semi-major axis by: epsilon = -mu / (2*a)

Rearranging: a = -mu / (2*epsilon)

Parameters:

• mu (float) – Standard gravitational parameter (km^3/s^2).

• specific_energy (float) – Specific orbital energy (km^2/s^2).

Returns:

a – Semi-major axis (km). - a > 0 for elliptical orbits (epsilon < 0) - a < 0 for hyperbolic orbits (epsilon > 0) - a → ∞ for parabolic orbits (epsilon = 0)

Return type:

float

Raises:

ValueError – If specific energy is exactly zero (parabolic case).

pytcl.astronomical.hyperbolic_asymptote_angle(e)

Compute the asymptote angle for a hyperbolic orbit.

For a hyperbolic orbit with eccentricity e, the true anomaly asymptotically approaches ±nu_inf where: cos(nu_inf) = -1/e

The asymptote angle is nu_inf.

Parameters:

e (float) – Eccentricity (e > 1 for hyperbolic).

Returns:

nu_inf – Asymptote angle (radians), in (0, pi).

Return type:

float

Raises:

ValueError – If eccentricity is not hyperbolic.

pytcl.astronomical.hyperbolic_deflection_angle(e)

Compute the deflection angle for a hyperbolic orbit.

The deflection angle is the angle through which the velocity vector is deflected from its asymptotic direction: delta = pi - 2*nu_inf = pi - 2*arccos(-1/e)

Parameters:

e (float) – Eccentricity (e > 1 for hyperbolic).

Returns:

delta – Deflection angle (radians), in (0, pi).

Return type:

float

Raises:

ValueError – If eccentricity is not hyperbolic.

pytcl.astronomical.eccentricity_vector(r, v, mu)

Compute eccentricity vector from position and velocity.

The eccentricity vector e is defined as: e = (v^2/mu - 1/r) * r - (r·v/mu) * v

This works for all orbit types: elliptical, parabolic, and hyperbolic.

Parameters:

• r (ndarray) – Position vector (km), shape (3,).

• v (ndarray) – Velocity vector (km/s), shape (3,).

• mu (float) – Standard gravitational parameter (km^3/s^2).

Returns:

e – Eccentricity vector, shape (3,).

Return type:

ndarray

Orbital Mechanics

Kepler propagation, orbital elements, and state conversions.

Two-body orbital mechanics and Kepler’s equation.

This module provides functions for two-body orbit propagation, Kepler’s equation solvers, and orbital element conversions.

References

[1]

Vallado, D. A., “Fundamentals of Astrodynamics and Applications,” 4th ed., Microcosm Press, 2013.

[2]

Curtis, H. D., “Orbital Mechanics for Engineering Students,” 3rd ed., Butterworth-Heinemann, 2014.

class pytcl.astronomical.orbital_mechanics.OrbitalElements(a, e, i, raan, omega, nu)

Bases: NamedTuple

Classical (Keplerian) orbital elements.

a

Semi-major axis (km). Negative for hyperbolic orbits.

Type:

float

e

Eccentricity. 0=circular, 0<e<1=elliptic, e=1=parabolic, e>1=hyperbolic.

Type:

float

i

Inclination (radians), 0 to pi.

Type:

float

raan

Right ascension of ascending node (radians), 0 to 2*pi.

Type:

float

omega

Argument of periapsis (radians), 0 to 2*pi.

Type:

float

nu

True anomaly (radians), 0 to 2*pi.

Type:

float

a: float

Alias for field number 0

e: float

Alias for field number 1

i: float

Alias for field number 2

raan: float

Alias for field number 3

omega: float

Alias for field number 4

nu: float

Alias for field number 5

class pytcl.astronomical.orbital_mechanics.StateVector(r, v)

Bases: NamedTuple

Cartesian state vector.

r

Position vector (km), shape (3,).

Type:

ndarray

v

Velocity vector (km/s), shape (3,).

Type:

ndarray

r: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 0

v: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 1

pytcl.astronomical.orbital_mechanics.mean_to_eccentric_anomaly(M, e, tol=1e-12, max_iter=50)

Solve Kepler’s equation: M = E - e*sin(E).

Uses Newton-Raphson iteration to find eccentric anomaly E given mean anomaly M and eccentricity e.

Parameters:

• M (float) – Mean anomaly (radians).

• e (float) – Eccentricity (0 <= e < 1 for elliptic orbits).

• tol (float, optional) – Convergence tolerance. Default 1e-12.

• max_iter (int, optional) – Maximum iterations. Default 50.

Returns:

E – Eccentric anomaly (radians).

Return type:

float

Raises:

ValueError – If eccentricity is not in valid range or iteration fails.

Examples

>>> import numpy as np
>>> E = mean_to_eccentric_anomaly(np.pi/4, 0.5)
>>> print(f"E = {np.degrees(E):.4f} degrees")

pytcl.astronomical.orbital_mechanics.mean_to_hyperbolic_anomaly(M, e, tol=1e-12, max_iter=50)

Solve hyperbolic Kepler’s equation: M = e*sinh(H) - H.

Uses Newton-Raphson iteration to find hyperbolic anomaly H given mean anomaly M and eccentricity e.

Parameters:

• M (float) – Mean anomaly (radians).

• e (float) – Eccentricity (e > 1 for hyperbolic orbits).

• tol (float, optional) – Convergence tolerance. Default 1e-12.

• max_iter (int, optional) – Maximum iterations. Default 50.

Returns:

H – Hyperbolic anomaly (radians).

Return type:

float

Examples

>>> H = mean_to_hyperbolic_anomaly(1.0, 1.5)
>>> abs(1.5 * np.sinh(H) - H - 1.0) < 1e-10
True

pytcl.astronomical.orbital_mechanics.eccentric_to_true_anomaly(E, e)

Convert eccentric anomaly to true anomaly.

Parameters:

• E (float) – Eccentric anomaly (radians).

• e (float) – Eccentricity.

Returns:

nu – True anomaly (radians), in [0, 2*pi).

Return type:

float

Examples

>>> nu = eccentric_to_true_anomaly(np.pi/4, 0.5)
>>> 0 <= nu < 2 * np.pi
True

pytcl.astronomical.orbital_mechanics.true_to_eccentric_anomaly(nu, e)

Convert true anomaly to eccentric anomaly.

Parameters:

• nu (float) – True anomaly (radians).

• e (float) – Eccentricity.

Returns:

E – Eccentric anomaly (radians), in [0, 2*pi).

Return type:

float

Examples

>>> E = true_to_eccentric_anomaly(np.pi/3, 0.5)
>>> 0 <= E < 2 * np.pi
True

pytcl.astronomical.orbital_mechanics.hyperbolic_to_true_anomaly(H, e)

Convert hyperbolic anomaly to true anomaly.

Parameters:

• H (float) – Hyperbolic anomaly (radians).

• e (float) – Eccentricity (e > 1).

Returns:

nu – True anomaly (radians).

Return type:

float

Examples

>>> nu = hyperbolic_to_true_anomaly(0.5, 1.5)
>>> isinstance(nu, float)
True

pytcl.astronomical.orbital_mechanics.true_to_hyperbolic_anomaly(nu, e)

Convert true anomaly to hyperbolic anomaly.

Parameters:

• nu (float) – True anomaly (radians).

• e (float) – Eccentricity (e > 1).

Returns:

H – Hyperbolic anomaly (radians).

Return type:

float

pytcl.astronomical.orbital_mechanics.eccentric_to_mean_anomaly(E, e)

Convert eccentric anomaly to mean anomaly (Kepler’s equation).

Parameters:

• E (float) – Eccentric anomaly (radians).

• e (float) – Eccentricity.

Returns:

M – Mean anomaly (radians).

Return type:

float

Examples

>>> M = eccentric_to_mean_anomaly(np.pi/4, 0.5)
>>> 0 <= M < 2 * np.pi
True

pytcl.astronomical.orbital_mechanics.mean_to_true_anomaly(M, e)

Convert mean anomaly to true anomaly.

Parameters:

• M (float) – Mean anomaly (radians).

• e (float) – Eccentricity.

Returns:

nu – True anomaly (radians).

Return type:

float

Examples

>>> nu = mean_to_true_anomaly(np.pi/4, 0.1)
>>> 0 <= nu < 2 * np.pi
True

pytcl.astronomical.orbital_mechanics.true_to_mean_anomaly(nu, e)

Convert true anomaly to mean anomaly.

Parameters:

• nu (float) – True anomaly (radians).

• e (float) – Eccentricity.

Returns:

M – Mean anomaly (radians).

Return type:

float

pytcl.astronomical.orbital_mechanics.orbital_elements_to_state(elements, mu=398600.4418)

Convert orbital elements to Cartesian state vector.

Parameters:

• elements (OrbitalElements) – Classical orbital elements.

• mu (float, optional) – Gravitational parameter (km^3/s^2). Default is Earth.

Returns:

state – Position and velocity vectors in inertial frame.

Return type:

StateVector

Examples

>>> elements = OrbitalElements(a=7000, e=0.01, i=0.5, raan=0, omega=0, nu=0)
>>> state = orbital_elements_to_state(elements)
>>> print(f"r = {state.r}")

pytcl.astronomical.orbital_mechanics.state_to_orbital_elements(state, mu=398600.4418)

Convert Cartesian state vector to orbital elements.

Parameters:

• state (StateVector) – Position and velocity vectors.

• mu (float, optional) – Gravitational parameter (km^3/s^2). Default is Earth.

Returns:

elements – Classical orbital elements.

Return type:

OrbitalElements

Examples

>>> r = np.array([7000, 0, 0])
>>> v = np.array([0, 7.5, 0])
>>> state = StateVector(r=r, v=v)
>>> elements = state_to_orbital_elements(state)

pytcl.astronomical.orbital_mechanics.kepler_propagate(elements, dt, mu=398600.4418)

Propagate orbital elements using Kepler’s equation.

This performs two-body propagation by advancing the mean anomaly and solving Kepler’s equation for the new true anomaly.

Parameters:

• elements (OrbitalElements) – Initial orbital elements.

• dt (float) – Time step (seconds).

• mu (float, optional) – Gravitational parameter (km^3/s^2).

Returns:

new_elements – Propagated orbital elements.

Return type:

OrbitalElements

Examples

>>> elements = OrbitalElements(a=7000, e=0.01, i=0.5, raan=0, omega=0, nu=0)
>>> new_elements = kepler_propagate(elements, 3600)  # 1 hour

pytcl.astronomical.orbital_mechanics.kepler_propagate_state(state, dt, mu=398600.4418)

Propagate state vector using Kepler’s equation.

Parameters:

• state (StateVector) – Initial state vector.

• dt (float) – Time step (seconds).

• mu (float, optional) – Gravitational parameter (km^3/s^2).

Returns:

new_state – Propagated state vector.

Return type:

StateVector

Examples

>>> r = np.array([7000.0, 0.0, 0.0])
>>> v = np.array([0.0, 7.5, 0.0])
>>> state = StateVector(r=r, v=v)
>>> new_state = kepler_propagate_state(state, 3600)
>>> np.linalg.norm(new_state.r) > 0
True

pytcl.astronomical.orbital_mechanics.orbital_period(a, mu=398600.4418)

Compute orbital period for elliptic orbit.

Parameters:

• a (float) – Semi-major axis (km).

• mu (float, optional) – Gravitational parameter (km^3/s^2).

Returns:

T – Orbital period (seconds).

Return type:

float

Examples

>>> T = orbital_period(7000)  # LEO satellite
>>> T / 60  # Convert to minutes
97.8...

pytcl.astronomical.orbital_mechanics.mean_motion(a, mu=398600.4418)

Compute mean motion.

Parameters:

• a (float) – Semi-major axis (km).

• mu (float, optional) – Gravitational parameter (km^3/s^2).

Returns:

n – Mean motion (radians/second).

Return type:

float

Examples

>>> n = mean_motion(42164)  # GEO orbit
>>> revs_per_day = n * 86400 / (2 * np.pi)
>>> abs(revs_per_day - 1.0) < 0.01  # Approximately 1 rev/day
True

pytcl.astronomical.orbital_mechanics.vis_viva(r, a, mu=398600.4418)

Compute orbital velocity using vis-viva equation.

Parameters:

• r (float) – Current orbital radius (km).

• a (float) – Semi-major axis (km).

• mu (float, optional) – Gravitational parameter (km^3/s^2).

Returns:

v – Orbital velocity (km/s).

Return type:

float

Examples

>>> v = vis_viva(7000, 7000)  # Circular orbit
>>> abs(v - circular_velocity(7000)) < 0.01
True

pytcl.astronomical.orbital_mechanics.specific_angular_momentum(state)

Compute specific angular momentum vector.

Parameters:

state (StateVector) – State vector.

Returns:

h – Specific angular momentum vector (km^2/s).

Return type:

ndarray

Examples

>>> r = np.array([7000.0, 0.0, 0.0])
>>> v = np.array([0.0, 7.5, 0.0])
>>> state = StateVector(r=r, v=v)
>>> h = specific_angular_momentum(state)
>>> h[2]  # Angular momentum in z-direction
52500.0

pytcl.astronomical.orbital_mechanics.specific_orbital_energy(state, mu=398600.4418)

Compute specific orbital energy.

Parameters:

• state (StateVector) – State vector.

• mu (float, optional) – Gravitational parameter (km^3/s^2).

Returns:

energy – Specific orbital energy (km^2/s^2). Negative for bound orbits, positive for escape trajectories.

Return type:

float

Examples

>>> r = np.array([7000.0, 0.0, 0.0])
>>> v = np.array([0.0, 7.5, 0.0])
>>> state = StateVector(r=r, v=v)
>>> energy = specific_orbital_energy(state)
>>> energy < 0  # Bound orbit
True

pytcl.astronomical.orbital_mechanics.flight_path_angle(state)

Compute flight path angle.

Parameters:

state (StateVector) – State vector.

Returns:

gamma – Flight path angle (radians). Positive when climbing, negative when descending.

Return type:

float

Examples

>>> r = np.array([7000.0, 0.0, 0.0])
>>> v = np.array([0.0, 7.5, 0.0])  # Tangential velocity
>>> state = StateVector(r=r, v=v)
>>> gamma = flight_path_angle(state)
>>> abs(gamma) < 0.01  # Nearly zero for circular motion
True

pytcl.astronomical.orbital_mechanics.periapsis_radius(a, e)

Compute periapsis radius.

Parameters:

• a (float) – Semi-major axis (km).

• e (float) – Eccentricity.

Returns:

r_p – Periapsis radius (km).

Return type:

float

Examples

>>> r_p = periapsis_radius(10000, 0.3)
>>> r_p
7000.0

pytcl.astronomical.orbital_mechanics.apoapsis_radius(a, e)

Compute apoapsis radius.

Parameters:

• a (float) – Semi-major axis (km).

• e (float) – Eccentricity.

Returns:

r_a – Apoapsis radius (km). Infinite for parabolic/hyperbolic orbits.

Return type:

float

Examples

>>> r_a = apoapsis_radius(10000, 0.3)
>>> r_a
13000.0

pytcl.astronomical.orbital_mechanics.time_since_periapsis(nu, a, e, mu=398600.4418)

Compute time since periapsis passage.

Parameters:

• nu (float) – True anomaly (radians).

• a (float) – Semi-major axis (km).

• e (float) – Eccentricity.

• mu (float, optional) – Gravitational parameter (km^3/s^2).

Returns:

t – Time since periapsis (seconds).

Return type:

float

Examples

>>> t = time_since_periapsis(np.pi, 7000, 0.1)  # At apoapsis
>>> T = orbital_period(7000)
>>> abs(t - T/2) < 1  # Approximately half the period
True

pytcl.astronomical.orbital_mechanics.orbit_radius(nu, a, e)

Compute orbital radius at given true anomaly.

Parameters:

• nu (float) – True anomaly (radians).

• a (float) – Semi-major axis (km).

• e (float) – Eccentricity.

Returns:

r – Orbital radius (km).

Return type:

float

Examples

>>> r = orbit_radius(0, 10000, 0.3)  # At periapsis
>>> r
7000.0
>>> r = orbit_radius(np.pi, 10000, 0.3)  # At apoapsis
>>> r
13000.0

pytcl.astronomical.orbital_mechanics.escape_velocity(r, mu=398600.4418)

Compute escape velocity at given radius.

Parameters:

• r (float) – Radial distance (km).

• mu (float, optional) – Gravitational parameter (km^3/s^2).

Returns:

v_esc – Escape velocity (km/s).

Return type:

float

Examples

>>> v_esc = escape_velocity(6378 + 400)  # At ISS altitude
>>> 10 < v_esc < 12  # About 11 km/s
True

pytcl.astronomical.orbital_mechanics.circular_velocity(r, mu=398600.4418)

Compute circular orbital velocity at given radius.

Parameters:

• r (float) – Orbital radius (km).

• mu (float, optional) – Gravitational parameter (km^3/s^2).

Returns:

v_circ – Circular velocity (km/s).

Return type:

float

Examples

>>> v_circ = circular_velocity(6378 + 400)  # At ISS altitude
>>> 7 < v_circ < 8  # About 7.7 km/s
True

Lambert Problem

Lambert problem solvers for orbit determination.

Lambert’s problem solver.

Lambert’s problem determines the orbit connecting two position vectors given the time of flight. This is fundamental for orbital transfer calculations, trajectory design, and orbit determination.

References

[1]

Vallado, D. A., “Fundamentals of Astrodynamics and Applications,” 4th ed., Microcosm Press, 2013.

[2]

Izzo, D., “Revisiting Lambert’s problem,” Celestial Mechanics and Dynamical Astronomy, 2015.

[3]

Gooding, R. H., “A procedure for the solution of Lambert’s orbital boundary-value problem,” Celestial Mechanics, 1990.

class pytcl.astronomical.lambert.LambertSolution(v1, v2, a, e, tof)

Bases: NamedTuple

Solution to Lambert’s problem.

v1

Velocity at first position (km/s), shape (3,).

Type:

ndarray

v2

Velocity at second position (km/s), shape (3,).

Type:

ndarray

a

Semi-major axis of transfer orbit (km).

Type:

float

e

Eccentricity of transfer orbit.

Type:

float

tof

Time of flight (seconds).

Type:

float

v1: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 0

v2: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 1

a: float

Alias for field number 2

e: float

Alias for field number 3

tof: float

Alias for field number 4

pytcl.astronomical.lambert.lambert_universal(r1, r2, tof, mu=398600.4418, prograde=True, low_path=True, max_iter=100, tol=1e-10)

Solve Lambert’s problem using universal variables.

Given two position vectors and time of flight, determine the transfer orbit connecting them.

Parameters:

• r1 (ndarray) – Initial position vector (km), shape (3,).

• r2 (ndarray) – Final position vector (km), shape (3,).

• tof (float) – Time of flight (seconds). Must be positive.

• mu (float, optional) – Gravitational parameter (km^3/s^2). Default is Earth.

• prograde (bool, optional) – If True, use prograde (counterclockwise) transfer. If False, use retrograde transfer. Default True.

• low_path (bool, optional) – If True, use low energy (short way) transfer. If False, use high energy (long way) transfer. Default True.

• max_iter (int, optional) – Maximum iterations. Default 100.

• tol (float, optional) – Convergence tolerance. Default 1e-10.

Returns:

solution – Solution containing velocities and orbital parameters.

Return type:

LambertSolution

Raises:

ValueError – If solution does not converge.

Examples

>>> r1 = np.array([5000, 10000, 2100])  # km
>>> r2 = np.array([-14600, 2500, 7000])  # km
>>> tof = 3600  # 1 hour
>>> sol = lambert_universal(r1, r2, tof)
>>> print(f"v1 = {sol.v1} km/s")

pytcl.astronomical.lambert.lambert_izzo(r1, r2, tof, mu=398600.4418, prograde=True, multi_rev=0, max_iter=100, tol=1e-10)

Solve Lambert’s problem using Izzo’s algorithm.

This is a more robust algorithm that handles multi-revolution transfers and edge cases better than the universal variable method.

Parameters:

• r1 (ndarray) – Initial position vector (km), shape (3,).

• r2 (ndarray) – Final position vector (km), shape (3,).

• tof (float) – Time of flight (seconds).

• mu (float, optional) – Gravitational parameter (km^3/s^2). Default is Earth.

• prograde (bool, optional) – If True, use prograde transfer. Default True.

• multi_rev (int, optional) – Number of complete revolutions. Default 0 (direct transfer).

• max_iter (int, optional) – Maximum iterations. Default 100.

• tol (float, optional) – Convergence tolerance. Default 1e-10.

Returns:

solution – Solution containing velocities and orbital parameters.

Return type:

LambertSolution

Notes

For multi-revolution transfers, there may be two solutions (low and high energy). This returns the low energy solution.

pytcl.astronomical.lambert.minimum_energy_transfer(r1, r2, mu=398600.4418, prograde=True)

Compute minimum energy transfer between two positions.

Parameters:

• r1 (ndarray) – Initial position vector (km).

• r2 (ndarray) – Final position vector (km).

• mu (float, optional) – Gravitational parameter (km^3/s^2).

• prograde (bool, optional) – If True, use prograde transfer.

Returns:

• tof_min (float) – Minimum energy time of flight (seconds).

• solution (LambertSolution) – Lambert solution at minimum energy.

Return type:

Tuple[float, LambertSolution]

pytcl.astronomical.lambert.hohmann_transfer(r1, r2, mu=398600.4418)

Compute Hohmann transfer between two circular orbits.

Parameters:

• r1 (float) – Initial orbit radius (km).

• r2 (float) – Final orbit radius (km).

• mu (float, optional) – Gravitational parameter (km^3/s^2).

Returns:

• dv1 (float) – Delta-v at first burn (km/s).

• dv2 (float) – Delta-v at second burn (km/s).

• tof (float) – Transfer time of flight (seconds).

Return type:

Tuple[float, float, float]

Examples

>>> dv1, dv2, tof = hohmann_transfer(6678, 42164)  # LEO to GEO
>>> print(f"Total dv = {dv1 + dv2:.3f} km/s")

pytcl.astronomical.lambert.bi_elliptic_transfer(r1, r2, r_intermediate, mu=398600.4418)

Compute bi-elliptic transfer between two circular orbits.

Parameters:

• r1 (float) – Initial orbit radius (km).

• r2 (float) – Final orbit radius (km).

• r_intermediate (float) – Intermediate apoapsis radius (km). Must be > max(r1, r2).

• mu (float, optional) – Gravitational parameter (km^3/s^2).

Returns:

• dv1 (float) – Delta-v at first burn (km/s).

• dv2 (float) – Delta-v at intermediate apoapsis (km/s).

• dv3 (float) – Delta-v at final circularization (km/s).

• tof (float) – Total transfer time (seconds).

Return type:

Tuple[float, float, float, float]

Reference Frames

Celestial and terrestrial reference frame transformations.

Reference frame transformations.

This module provides transformations between inertial and Earth-fixed reference frames, including precession, nutation, and Earth rotation.

Frames: - GCRF: Geocentric Celestial Reference Frame (inertial) - ITRF: International Terrestrial Reference Frame (Earth-fixed) - J2000: Mean equator and equinox of J2000.0 - TOD: True of Date (true equator and equinox) - MOD: Mean of Date (mean equator and equinox) - PEF: Pseudo Earth-Fixed (rotated by Earth rotation only)

References

[1]

Vallado, D. A., “Fundamentals of Astrodynamics and Applications,” 4th ed., Microcosm Press, 2013.

[2]

IERS Conventions (2010), IERS Technical Note No. 36.

[3]

Capitaine et al., “Expressions for IAU 2000 precession quantities,” A&A, 2003.

pytcl.astronomical.reference_frames.julian_centuries_j2000(jd)

Compute Julian centuries since J2000.0.

Parameters:

jd (float) – Julian date.

Returns:

T – Julian centuries since J2000.0.

Return type:

float

pytcl.astronomical.reference_frames.precession_angles_iau76(T)

Compute IAU 1976 precession angles.

Parameters:

T (float) – Julian centuries since J2000.0.

Returns:

• zeta (float) – Precession angle zeta (radians).

• theta (float) – Precession angle theta (radians).

• z (float) – Precession angle z (radians).

Return type:

Tuple[float, float, float]

pytcl.astronomical.reference_frames.precession_matrix_iau76(jd)

Compute IAU 1976 precession matrix from J2000 to date.

Parameters:

jd (float) – Julian date of epoch.

Returns:

P – Precession rotation matrix (3x3). Transforms from J2000 (GCRF) to mean of date.

Return type:

ndarray

Notes

Results are cached for repeated queries at the same epoch. Cache key is quantized to ~86ms precision.

pytcl.astronomical.reference_frames.nutation_angles_iau80(jd)

Compute IAU 1980 nutation angles (simplified).

Parameters:

jd (float) – Julian date.

Returns:

• dpsi (float) – Nutation in longitude (radians).

• deps (float) – Nutation in obliquity (radians).

Return type:

Tuple[float, float]

pytcl.astronomical.reference_frames.nutation_matrix(jd)

Compute nutation matrix.

Parameters:

jd (float) – Julian date.

Returns:

N – Nutation rotation matrix (3x3). Transforms from mean of date to true of date.

Return type:

ndarray

Notes

Results are cached for repeated queries at the same epoch. Cache key is quantized to ~86ms precision.

pytcl.astronomical.reference_frames.mean_obliquity_iau80(jd)

Compute mean obliquity of the ecliptic (IAU 1980).

Parameters:

jd (float) – Julian date.

Returns:

eps0 – Mean obliquity (radians).

Return type:

float

pytcl.astronomical.reference_frames.true_obliquity(jd_tt)

Compute true obliquity of the ecliptic.

Parameters:

jd_tt (float) – Julian date in TT.

Returns:

eps – True obliquity (radians).

Return type:

float

pytcl.astronomical.reference_frames.earth_rotation_angle(jd_ut1)

Compute Earth Rotation Angle (ERA).

Parameters:

jd_ut1 (float) – Julian date in UT1.

Returns:

theta – Earth rotation angle (radians), in [0, 2*pi).

Return type:

float

pytcl.astronomical.reference_frames.gmst_iau82(jd_ut1)

Compute Greenwich Mean Sidereal Time (GMST) using IAU 1982 model.

Parameters:

jd_ut1 (float) – Julian date in UT1.

Returns:

gmst – GMST (radians), in [0, 2*pi).

Return type:

float

pytcl.astronomical.reference_frames.gast_iau82(jd_ut1, jd_tt)

Compute Greenwich Apparent Sidereal Time (GAST).

Parameters:

• jd_ut1 (float) – Julian date in UT1.

• jd_tt (float) – Julian date in TT.

Returns:

gast – GAST (radians), in [0, 2*pi).

Return type:

float

pytcl.astronomical.reference_frames.sidereal_rotation_matrix(theta)

Compute Earth rotation matrix (sidereal rotation).

Parameters:

theta (float) – Sidereal angle (GMST or GAST) in radians.

Returns:

R – Rotation matrix (3x3). Transforms from true of date to pseudo Earth-fixed (PEF).

Return type:

ndarray

pytcl.astronomical.reference_frames.equation_of_equinoxes(jd_tt)

Compute equation of the equinoxes.

Parameters:

jd_tt (float) – Julian date in TT.

Returns:

eq_eq – Equation of the equinoxes (radians).

Return type:

float

pytcl.astronomical.reference_frames.polar_motion_matrix(xp, yp)

Compute polar motion rotation matrix.

Parameters:

• xp (float) – Polar motion x coordinate (radians).

• yp (float) – Polar motion y coordinate (radians).

Returns:

W – Polar motion matrix (3x3). Transforms from PEF to ITRF.

Return type:

ndarray

pytcl.astronomical.reference_frames.gcrf_to_itrf(r_gcrf, jd_ut1, jd_tt, xp=0.0, yp=0.0)

Transform position from GCRF (inertial) to ITRF (Earth-fixed).

Parameters:

• r_gcrf (ndarray) – Position in GCRF (km), shape (3,).

• jd_ut1 (float) – Julian date in UT1.

• jd_tt (float) – Julian date in TT.

• xp (float, optional) – Polar motion x (radians). Default 0.

• yp (float, optional) – Polar motion y (radians). Default 0.

Returns:

r_itrf – Position in ITRF (km), shape (3,).

Return type:

ndarray

Examples

>>> r_gcrf = np.array([5102.5096, 6123.01152, 6378.1363])
>>> r_itrf = gcrf_to_itrf(r_gcrf, 2453101.827406783, 2453101.828154745)

pytcl.astronomical.reference_frames.itrf_to_gcrf(r_itrf, jd_ut1, jd_tt, xp=0.0, yp=0.0)

Transform position from ITRF (Earth-fixed) to GCRF (inertial).

Parameters:

• r_itrf (ndarray) – Position in ITRF (km), shape (3,).

• jd_ut1 (float) – Julian date in UT1.

• jd_tt (float) – Julian date in TT.

• xp (float, optional) – Polar motion x (radians). Default 0.

• yp (float, optional) – Polar motion y (radians). Default 0.

Returns:

r_gcrf – Position in GCRF (km), shape (3,).

Return type:

ndarray

pytcl.astronomical.reference_frames.gcrf_to_pef(r_gcrf, jd_ut1, jd_tt)

Transform position from GCRF (inertial) to PEF (Earth-fixed, rotation only).

PEF (Pseudo-Earth Fixed) is an intermediate reference frame between GCRF and ITRF. It includes precession, nutation, and Earth rotation, but excludes polar motion.

Parameters:

• r_gcrf (ndarray) – Position in GCRF (km), shape (3,).

• jd_ut1 (float) – Julian date in UT1.

• jd_tt (float) – Julian date in TT.

Returns:

r_pef – Position in PEF (km), shape (3,).

Return type:

ndarray

Notes

The transformation chain is: GCRF -> MOD -> TOD -> PEF - Precession: GCRF -> MOD - Nutation: MOD -> TOD - Sidereal rotation: TOD -> PEF

See also

pef_to_gcrf

Inverse transformation

gcrf_to_itrf

Includes polar motion

References

[1]

Vallado et al., “Fundamentals of Astrodynamics and Applications”, 4th ed.

pytcl.astronomical.reference_frames.pef_to_gcrf(r_pef, jd_ut1, jd_tt)

Transform position from PEF (Earth-fixed, rotation only) to GCRF (inertial).

Inverse of gcrf_to_pef.

Parameters:

• r_pef (ndarray) – Position in PEF (km), shape (3,).

• jd_ut1 (float) – Julian date in UT1.

• jd_tt (float) – Julian date in TT.

Returns:

r_gcrf – Position in GCRF (km), shape (3,).

Return type:

ndarray

See also

gcrf_to_pef

Forward transformation

pytcl.astronomical.reference_frames.eci_to_ecef(r_eci, gmst)

Simple ECI to ECEF transformation (rotation only).

This is a simplified transformation that only accounts for Earth rotation, ignoring precession, nutation, and polar motion.

Parameters:

• r_eci (ndarray) – Position in ECI (km), shape (3,).

• gmst (float) – Greenwich Mean Sidereal Time (radians).

Returns:

r_ecef – Position in ECEF (km), shape (3,).

Return type:

ndarray

pytcl.astronomical.reference_frames.ecef_to_eci(r_ecef, gmst)

Simple ECEF to ECI transformation (rotation only).

Parameters:

• r_ecef (ndarray) – Position in ECEF (km), shape (3,).

• gmst (float) – Greenwich Mean Sidereal Time (radians).

Returns:

r_eci – Position in ECI (km), shape (3,).

Return type:

ndarray

pytcl.astronomical.reference_frames.ecliptic_to_equatorial(r_ecl, obliquity)

Transform from ecliptic to equatorial coordinates.

Parameters:

• r_ecl (ndarray) – Position in ecliptic frame (km), shape (3,).

• obliquity (float) – Obliquity of the ecliptic (radians).

Returns:

r_eq – Position in equatorial frame (km), shape (3,).

Return type:

ndarray

pytcl.astronomical.reference_frames.equatorial_to_ecliptic(r_eq, obliquity)

Transform from equatorial to ecliptic coordinates.

Parameters:

• r_eq (ndarray) – Position in equatorial frame (km), shape (3,).

• obliquity (float) – Obliquity of the ecliptic (radians).

Returns:

r_ecl – Position in ecliptic frame (km), shape (3,).

Return type:

ndarray

pytcl.astronomical.reference_frames.teme_to_pef(r_teme, jd_ut1)

Transform position from TEME to PEF (Pseudo Earth-Fixed).

TEME is the True Equator, Mean Equinox frame used by SGP4. This transformation applies only the GMST rotation.

Parameters:

• r_teme (ndarray) – Position in TEME frame (km), shape (3,).

• jd_ut1 (float) – Julian date in UT1.

Returns:

r_pef – Position in PEF frame (km), shape (3,).

Return type:

ndarray

pytcl.astronomical.reference_frames.pef_to_teme(r_pef, jd_ut1)

Transform position from PEF to TEME.

Parameters:

• r_pef (ndarray) – Position in PEF frame (km), shape (3,).

• jd_ut1 (float) – Julian date in UT1.

Returns:

r_teme – Position in TEME frame (km), shape (3,).

Return type:

ndarray

pytcl.astronomical.reference_frames.teme_to_itrf(r_teme, jd_ut1, xp=0.0, yp=0.0)

Transform position from TEME to ITRF (Earth-fixed).

TEME is the True Equator, Mean Equinox frame used by SGP4/SDP4. This is the frame in which TLE-propagated positions are expressed.

Parameters:

• r_teme (ndarray) – Position in TEME frame (km), shape (3,).

• jd_ut1 (float) – Julian date in UT1.

• xp (float, optional) – Polar motion x (radians). Default 0.

• yp (float, optional) – Polar motion y (radians). Default 0.

Returns:

r_itrf – Position in ITRF frame (km), shape (3,).

Return type:

ndarray

Notes

TEME is a quasi-inertial frame that uses the mean equinox instead of the true equinox. The transformation sequence is:

TEME -> PEF (via GMST rotation) -> ITRF (via polar motion)

Examples

>>> from pytcl.astronomical.sgp4 import sgp4_propagate
>>> from pytcl.astronomical.tle import parse_tle
>>> tle = parse_tle(line1, line2)
>>> state = sgp4_propagate(tle, 0.0)
>>> r_itrf = teme_to_itrf(state.r, jd_ut1)

pytcl.astronomical.reference_frames.itrf_to_teme(r_itrf, jd_ut1, xp=0.0, yp=0.0)

Transform position from ITRF to TEME.

Parameters:

• r_itrf (ndarray) – Position in ITRF frame (km), shape (3,).

• jd_ut1 (float) – Julian date in UT1.

• xp (float, optional) – Polar motion x (radians). Default 0.

• yp (float, optional) – Polar motion y (radians). Default 0.

Returns:

r_teme – Position in TEME frame (km), shape (3,).

Return type:

ndarray

pytcl.astronomical.reference_frames.teme_to_gcrf(r_teme, jd_tt)

Transform position from TEME to GCRF (inertial).

This transformation accounts for the difference between the mean and true equinox (equation of equinoxes) and then applies precession and nutation to go from TOD to GCRF.

Parameters:

• r_teme (ndarray) – Position in TEME frame (km), shape (3,).

• jd_tt (float) – Julian date in TT (Terrestrial Time).

Returns:

r_gcrf – Position in GCRF frame (km), shape (3,).

Return type:

ndarray

Notes

The transformation sequence is:

TEME -> TOD (via equation of equinoxes) TOD -> MOD (via nutation, inverse) MOD -> GCRF (via precession, inverse)

Examples

>>> state = sgp4_propagate(tle, 60.0)
>>> r_gcrf = teme_to_gcrf(state.r, jd_tt)

pytcl.astronomical.reference_frames.gcrf_to_teme(r_gcrf, jd_tt)

Transform position from GCRF to TEME.

Parameters:

• r_gcrf (ndarray) – Position in GCRF frame (km), shape (3,).

• jd_tt (float) – Julian date in TT.

Returns:

r_teme – Position in TEME frame (km), shape (3,).

Return type:

ndarray

pytcl.astronomical.reference_frames.teme_to_itrf_with_velocity(r_teme, v_teme, jd_ut1, xp=0.0, yp=0.0)

Transform position and velocity from TEME to ITRF.

This properly accounts for the velocity transformation including the Earth’s rotation rate.

Parameters:

• r_teme (ndarray) – Position in TEME frame (km), shape (3,).

• v_teme (ndarray) – Velocity in TEME frame (km/s), shape (3,).

• jd_ut1 (float) – Julian date in UT1.

• xp (float, optional) – Polar motion x (radians). Default 0.

• yp (float, optional) – Polar motion y (radians). Default 0.

Returns:

• r_itrf (ndarray) – Position in ITRF frame (km), shape (3,).

• v_itrf (ndarray) – Velocity in ITRF frame (km/s), shape (3,).

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

pytcl.astronomical.reference_frames.itrf_to_teme_with_velocity(r_itrf, v_itrf, jd_ut1, xp=0.0, yp=0.0)

Transform position and velocity from ITRF to TEME.

Parameters:

• r_itrf (ndarray) – Position in ITRF frame (km), shape (3,).

• v_itrf (ndarray) – Velocity in ITRF frame (km/s), shape (3,).

• jd_ut1 (float) – Julian date in UT1.

• xp (float, optional) – Polar motion x (radians). Default 0.

• yp (float, optional) – Polar motion y (radians). Default 0.

Returns:

• r_teme (ndarray) – Position in TEME frame (km), shape (3,).

• v_teme (ndarray) – Velocity in TEME frame (km/s), shape (3,).

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

pytcl.astronomical.reference_frames.gcrf_to_mod(r_gcrf, jd_tt)

Transform position from GCRF to MOD (Mean of Date).

MOD is the mean equator and mean equinox of date frame. This applies only the precession transformation.

Parameters:

• r_gcrf (ndarray) – Position in GCRF frame (km), shape (3,).

• jd_tt (float) – Julian date in TT (Terrestrial Time).

Returns:

r_mod – Position in MOD frame (km), shape (3,).

Return type:

ndarray

Notes

MOD is a legacy frame convention. For most modern applications, GCRF (J2000) is preferred. MOD was historically used in older software and publications.

The transformation is simply the precession matrix:

r_mod = P @ r_gcrf

See also

mod_to_gcrf

Inverse transformation.

gcrf_to_tod

Includes nutation for true of date.

pytcl.astronomical.reference_frames.mod_to_gcrf(r_mod, jd_tt)

Transform position from MOD (Mean of Date) to GCRF.

Parameters:

• r_mod (ndarray) – Position in MOD frame (km), shape (3,).

• jd_tt (float) – Julian date in TT.

Returns:

r_gcrf – Position in GCRF frame (km), shape (3,).

Return type:

ndarray

See also

gcrf_to_mod

Forward transformation.

pytcl.astronomical.reference_frames.gcrf_to_tod(r_gcrf, jd_tt)

Transform position from GCRF to TOD (True of Date).

TOD is the true equator and true equinox of date frame. This applies both precession and nutation transformations.

Parameters:

• r_gcrf (ndarray) – Position in GCRF frame (km), shape (3,).

• jd_tt (float) – Julian date in TT (Terrestrial Time).

Returns:

r_tod – Position in TOD frame (km), shape (3,).

Return type:

ndarray

Notes

TOD is a legacy frame convention. The transformation is:

r_mod = P @ r_gcrf r_tod = N @ r_mod

where P is the precession matrix and N is the nutation matrix.

See also

tod_to_gcrf

Inverse transformation.

gcrf_to_mod

Mean of date (without nutation).

pytcl.astronomical.reference_frames.tod_to_gcrf(r_tod, jd_tt)

Transform position from TOD (True of Date) to GCRF.

Parameters:

• r_tod (ndarray) – Position in TOD frame (km), shape (3,).

• jd_tt (float) – Julian date in TT.

Returns:

r_gcrf – Position in GCRF frame (km), shape (3,).

Return type:

ndarray

See also

gcrf_to_tod

Forward transformation.

pytcl.astronomical.reference_frames.mod_to_tod(r_mod, jd_tt)

Transform position from MOD (Mean of Date) to TOD (True of Date).

This applies only the nutation transformation.

Parameters:

• r_mod (ndarray) – Position in MOD frame (km), shape (3,).

• jd_tt (float) – Julian date in TT.

Returns:

r_tod – Position in TOD frame (km), shape (3,).

Return type:

ndarray

Notes

The transformation is simply the nutation matrix:

r_tod = N @ r_mod

See also

tod_to_mod

Inverse transformation.

pytcl.astronomical.reference_frames.tod_to_mod(r_tod, jd_tt)

Transform position from TOD (True of Date) to MOD (Mean of Date).

Parameters:

• r_tod (ndarray) – Position in TOD frame (km), shape (3,).

• jd_tt (float) – Julian date in TT.

Returns:

r_mod – Position in MOD frame (km), shape (3,).

Return type:

ndarray

See also

mod_to_tod

Forward transformation.

pytcl.astronomical.reference_frames.tod_to_itrf(r_tod, jd_ut1, jd_tt=None, xp=0.0, yp=0.0)

Transform position from TOD (True of Date) to ITRF.

Parameters:

• r_tod (ndarray) – Position in TOD frame (km), shape (3,).

• jd_ut1 (float) – Julian date in UT1.

• jd_tt (float, optional) – Julian date in TT. If not provided, assumed equal to jd_ut1.

• xp (float, optional) – Polar motion x (radians). Default 0.

• yp (float, optional) – Polar motion y (radians). Default 0.

Returns:

r_itrf – Position in ITRF frame (km), shape (3,).

Return type:

ndarray

Notes

The transformation applies the sidereal rotation (using GAST) and polar motion:

r_pef = R(GAST) @ r_tod r_itrf = W @ r_pef

See also

itrf_to_tod

Inverse transformation.

pytcl.astronomical.reference_frames.itrf_to_tod(r_itrf, jd_ut1, jd_tt=None, xp=0.0, yp=0.0)

Transform position from ITRF to TOD (True of Date).

Parameters:

• r_itrf (ndarray) – Position in ITRF frame (km), shape (3,).

• jd_ut1 (float) – Julian date in UT1.

• jd_tt (float, optional) – Julian date in TT. If not provided, assumed equal to jd_ut1.

• xp (float, optional) – Polar motion x (radians). Default 0.

• yp (float, optional) – Polar motion y (radians). Default 0.

Returns:

r_tod – Position in TOD frame (km), shape (3,).

Return type:

ndarray

See also

tod_to_itrf

Forward transformation.

pytcl.astronomical.reference_frames.clear_transformation_cache()

Clear cached transformation matrices.

Call this function to clear all cached precession and nutation matrices. Useful when memory is constrained or after processing a batch of observations at different epochs.

pytcl.astronomical.reference_frames.get_cache_info()

Get cache statistics for transformation matrices.

Returns:

Dictionary with ‘precession’ and ‘nutation’ keys, each containing CacheInfo namedtuple with hits, misses, maxsize, currsize.

Return type:

dict

Time Systems

Astronomical time systems and conversions.

Time system conversions for astronomical and tracking applications.

This module provides conversions between various time systems: - UTC (Coordinated Universal Time) - TAI (International Atomic Time) - TT (Terrestrial Time) - GPS (GPS Time) - Julian Date (JD) and Modified Julian Date (MJD) - Unix time - Sidereal time (GMST, GAST)

pytcl.astronomical.time_systems.cal_to_jd(year, month, day, hour=0, minute=0, second=0.0)

Convert calendar date to Julian Date.

Parameters:

• year (int) – Year (negative for BCE).

• month (int) – Month (1-12).

• day (int) – Day of month (1-31).

• hour (int, optional) – Hour (0-23).

• minute (int, optional) – Minute (0-59).

• second (float, optional) – Second (0-59.999…).

Returns:

Julian Date.

Return type:

float

Examples

>>> cal_to_jd(2000, 1, 1, 12, 0, 0)  # J2000.0
2451545.0
>>> cal_to_jd(1970, 1, 1, 0, 0, 0)  # Unix epoch
2440587.5

Notes

Uses the algorithm from Meeus, “Astronomical Algorithms”, 2nd ed. Valid for dates after October 15, 1582 (Gregorian calendar).

pytcl.astronomical.time_systems.jd_to_cal(jd)

Convert Julian Date to calendar date.

Parameters:

jd (float) – Julian Date.

Returns:

(year, month, day, hour, minute, second).

Return type:

tuple

Examples

>>> jd_to_cal(2451545.0)  # J2000.0
(2000, 1, 1, 12, 0, 0.0)

pytcl.astronomical.time_systems.mjd_to_jd(mjd)

Convert Modified Julian Date to Julian Date.

Parameters:

mjd (float) – Modified Julian Date.

Returns:

Julian Date.

Return type:

float

Examples

>>> mjd = 44239.0  # 1980-01-01
>>> jd = mjd_to_jd(mjd)
>>> jd
2444239.5

Notes

MJD = JD - 2400000.5

pytcl.astronomical.time_systems.jd_to_mjd(jd)

Convert Julian Date to Modified Julian Date.

Parameters:

jd (float) – Julian Date.

Returns:

Modified Julian Date.

Return type:

float

Examples

>>> jd = 2444239.5  # 1980-01-01
>>> mjd = jd_to_mjd(jd)
>>> mjd
44239.0

pytcl.astronomical.time_systems.utc_to_tai(year, month, day, hour=0, minute=0, second=0.0)

Convert UTC to TAI (Julian Date).

Parameters:

• year (int) – Calendar date.

• month (int) – Calendar date.

• day (int) – Calendar date.

• hour (int, optional) – Time of day.

• minute (int, optional) – Time of day.

• second (float, optional) – Seconds (including fractional).

Returns:

TAI as Julian Date.

Return type:

float

Notes

TAI = UTC + leap_seconds

pytcl.astronomical.time_systems.tai_to_utc(jd_tai)

Convert TAI (Julian Date) to UTC.

Parameters:

jd_tai (float) – TAI as Julian Date.

Returns:

• jd_utc (float) – UTC as Julian Date.

• leap_seconds (int) – Number of leap seconds applied.

Return type:

Tuple[float, int]

Notes

This is an approximate conversion that may have small errors near leap second boundaries.

pytcl.astronomical.time_systems.tai_to_tt(jd_tai)

Convert TAI (Julian Date) to TT (Terrestrial Time).

Parameters:

jd_tai (float) – TAI as Julian Date.

Returns:

TT as Julian Date.

Return type:

float

Notes

TT = TAI + 32.184 seconds (constant offset)

pytcl.astronomical.time_systems.tt_to_tai(jd_tt)

Convert TT (Terrestrial Time) to TAI (Julian Date).

Parameters:

jd_tt (float) – TT as Julian Date.

Returns:

TAI as Julian Date.

Return type:

float

pytcl.astronomical.time_systems.utc_to_tt(year, month, day, hour=0, minute=0, second=0.0)

Convert UTC to TT (Terrestrial Time).

Parameters:

• year (int) – Calendar date.

• month (int) – Calendar date.

• day (int) – Calendar date.

• hour (int, optional) – Time of day.

• minute (int, optional) – Time of day.

• second (float, optional) – Seconds.

Returns:

TT as Julian Date.

Return type:

float

Notes

TT = UTC + leap_seconds + 32.184

pytcl.astronomical.time_systems.tt_to_utc(jd_tt)

Convert TT to UTC.

Parameters:

jd_tt (float) – TT as Julian Date.

Returns:

• jd_utc (float) – UTC as Julian Date.

• leap_seconds (int) – Number of leap seconds applied.

Return type:

Tuple[float, int]

pytcl.astronomical.time_systems.tai_to_gps(jd_tai)

Convert TAI to GPS time.

Parameters:

jd_tai (float) – TAI as Julian Date.

Returns:

GPS time as Julian Date.

Return type:

float

Notes

GPS time = TAI - 19 seconds (offset at GPS epoch) GPS time does not include leap seconds after 1980.

pytcl.astronomical.time_systems.gps_to_tai(jd_gps)

Convert GPS time to TAI.

Parameters:

jd_gps (float) – GPS time as Julian Date.

Returns:

TAI as Julian Date.

Return type:

float

pytcl.astronomical.time_systems.utc_to_gps(year, month, day, hour=0, minute=0, second=0.0)

Convert UTC to GPS time.

Parameters:

• year (int) – Calendar date.

• month (int) – Calendar date.

• day (int) – Calendar date.

• hour (int, optional) – Time of day.

• minute (int, optional) – Time of day.

• second (float, optional) – Seconds.

Returns:

GPS time as Julian Date.

Return type:

float

pytcl.astronomical.time_systems.gps_to_utc(jd_gps)

Convert GPS time to UTC.

Parameters:

jd_gps (float) – GPS time as Julian Date.

Returns:

• jd_utc (float) – UTC as Julian Date.

• leap_seconds (int) – Number of leap seconds applied.

Return type:

Tuple[float, int]

pytcl.astronomical.time_systems.unix_to_jd(unix_time)

Convert Unix timestamp to Julian Date.

Parameters:

unix_time (float) – Unix timestamp (seconds since 1970-01-01 00:00:00 UTC).

Returns:

Julian Date.

Return type:

float

Examples

>>> unix_to_jd(0.0)  # Unix epoch
2440587.5

pytcl.astronomical.time_systems.jd_to_unix(jd)

Convert Julian Date to Unix timestamp.

Parameters:

jd (float) – Julian Date.

Returns:

Unix timestamp.

Return type:

float

Examples

>>> jd = 2440587.5  # 1970-01-01 00:00:00 UTC
>>> unix_to_jd = jd_to_unix(jd)
>>> unix_to_jd
0.0

pytcl.astronomical.time_systems.gps_week_seconds(jd_gps)

Convert GPS Julian Date to GPS week and seconds of week.

Parameters:

jd_gps (float) – GPS time as Julian Date.

Returns:

• week (int) – GPS week number.

• seconds (float) – Seconds into the week.

Return type:

Tuple[int, float]

Examples

>>> gps_week_seconds(JD_GPS_EPOCH)
(0, 0.0)

pytcl.astronomical.time_systems.gps_week_to_utc(week, seconds)

Convert GPS week and seconds to UTC.

Parameters:

• week (int) – GPS week number.

• seconds (float) – Seconds into the week.

Returns:

• jd_utc (float) – UTC as Julian Date.

• leap_seconds (int) – Number of leap seconds applied.

Return type:

Tuple[float, int]

pytcl.astronomical.time_systems.gmst(jd_ut1)

Compute Greenwich Mean Sidereal Time.

Parameters:

jd_ut1 (float) – UT1 as Julian Date.

Returns:

GMST in radians.

Return type:

float

Notes

Uses the IAU 1982 expression for GMST.

References

[1]

Explanatory Supplement to the Astronomical Almanac, 3rd ed.

pytcl.astronomical.time_systems.gast(jd_ut1, dpsi=0.0, eps=0.0)

Compute Greenwich Apparent Sidereal Time.

Parameters:

• jd_ut1 (float) – UT1 as Julian Date.

• dpsi (float, optional) – Nutation in longitude (radians). Default 0.

• eps (float, optional) – Mean obliquity of ecliptic (radians). Default uses approximate value.

Returns:

GAST in radians.

Return type:

float

Notes

GAST = GMST + equation of equinoxes The equation of equinoxes = dpsi * cos(eps)

For high precision applications, nutation parameters should be computed from the IAU 2006/2000A precession-nutation model.

pytcl.astronomical.time_systems.get_leap_seconds(year, month, day)

Get the number of leap seconds (TAI-UTC offset) for a given date.

Parameters:

• year (int) – Year.

• month (int) – Month (1-12).

• day (int) – Day of month.

Returns:

TAI-UTC offset in seconds.

Return type:

int

Examples

>>> get_leap_seconds(2020, 1, 1)
37
>>> get_leap_seconds(1980, 1, 6)
19

class pytcl.astronomical.time_systems.LeapSecondTable

Bases: object

Table of leap seconds for UTC-TAI conversion.

The leap second table contains the dates when leap seconds were added and the cumulative TAI-UTC offset.

entries

List of (year, month, day, tai_utc_offset) tuples.

Type:

list of tuple

Notes

This table must be updated when new leap seconds are announced. Last update: 2017-01-01 (37 seconds).

As of 2024, no new leap seconds have been added since 2017.

__init__()

get_offset(year, month, day)

Get TAI-UTC offset for a given date.

Parameters:

• year (int) – Year.

• month (int) – Month (1-12).

• day (int) – Day of month.

Returns:

TAI-UTC offset in seconds.

Return type:

int

JPL Ephemerides

High-precision celestial body positions using JPL Development Ephemeris files.

JPL Ephemerides for High-Precision Celestial Mechanics

This module provides access to JPL Development Ephemeris (DE) files for computing high-precision positions and velocities of celestial bodies (Sun, Moon, planets).

The module leverages the jplephem library, which provides optimized Fortran-based interpolation of ephemeris kernels. Multiple DE versions are supported (DE405, DE430, DE432s, DE440).

Constants

AU_PER_KMfloat

Astronomical Unit in kilometers (1 AU = 149597870.7 km)

KM_PER_DAY_TO_AU_PER_DAYfloat

Conversion factor for velocity from km/day to AU/day

Examples

>>> from pytcl.astronomical.ephemerides import DEEphemeris
>>> from datetime import datetime
>>>
>>> # Load ephemeris (auto-downloads if needed)
>>> eph = DEEphemeris(version='DE440')
>>>
>>> # Query Sun position (AU)
>>> jd = 2451545.0  # J2000.0
>>> r_sun, v_sun = eph.sun_position(jd)
>>> print(f"Sun distance: {np.linalg.norm(r_sun):.6f} AU")
Sun distance: 0.983327 AU
>>>
>>> # Query Moon position
>>> r_moon, v_moon = eph.moon_position(jd)

Notes

• Ephemeris files are auto-downloaded to ~/.jplephem/ on first use

• Time input is Julian Day (JD) in Terrestrial Time (TT) scale

• Positions returned in AU, velocities in AU/day in ICRF frame

• For highest precision, use DE440 (latest release) or DE432s (2013)

References

[1]

Standish, E. M. (1995). “Report of the IAU WGAS Sub-group on Numerical Standards”. In Highlights of Astronomy (Vol. 10).

[2]

Folkner, W. M., Williams, J. G., Boggs, D. H., Park, R. S., & Kuchynka, P. (2014). “The Planetary and Lunar Ephemeris DE430 and DE431”. Interplanetary Network Progress Report, 42(196), 1-81.

class pytcl.astronomical.ephemerides.DEEphemeris(version='DE440')

Bases: object

High-precision JPL Development Ephemeris kernel wrapper.

This class manages access to JPL ephemeris files and provides methods for querying positions and velocities of celestial bodies.

Parameters:

version ({'DE405', 'DE430', 'DE432s', 'DE440'}, optional) – Ephemeris version to load. Default is ‘DE440’ (latest). - DE440: Latest JPL release (2020), covers 1550-2650 - DE432s: High-precision version (2013), covers 1350-3000 - DE430: Earlier release (2013), covers 1550-2650 - DE405: Older version (1998), compact, covers 1600-2200

version

Ephemeris version identifier

Type:

str

kernel

Loaded ephemeris kernel object

Type:

jplephem.SpiceKernel

_cache

Cache for frequently accessed positions

Type:

dict

Raises:

• DependencyError – If jplephem is not installed

• ValueError – If version is not recognized

Examples

>>> eph = DEEphemeris(version='DE440')
>>> r_sun, v_sun = eph.sun_position(2451545.0)

__init__(version='DE440')

Initialize ephemeris kernel.

Parameters:

version (str, optional) – Ephemeris version (default: ‘DE440’)

property kernel: object | None

Lazy-load ephemeris kernel on first access.

Note: This requires jplephem to be installed and the kernel file to be available locally or downloadable from the JPL servers.

sun_position(jd, frame='icrf')

Compute Sun position and velocity.

Parameters:

• jd (float) – Julian Day in Terrestrial Time (TT)

• frame ({'icrf', 'ecliptic'}, optional) – Coordinate frame (default: ‘icrf’). - ‘icrf’: International Celestial Reference Frame - ‘ecliptic’: Ecliptic coordinate system (J2000.0)

Returns:

• position (ndarray, shape (3,)) – Sun position in AU

• velocity (ndarray, shape (3,)) – Sun velocity in AU/day

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Notes

The Sun’s position is computed relative to the Solar System Barycenter (SSB) in the ICRF frame.

Examples

>>> eph = DEEphemeris()
>>> r, v = eph.sun_position(2451545.0)
>>> print(f"Distance: {np.linalg.norm(r):.6f} AU")

moon_position(jd, frame='icrf')

Compute Moon position and velocity.

Parameters:

• jd (float) – Julian Day in Terrestrial Time (TT)

• frame ({'icrf', 'ecliptic', 'earth_centered'}, optional) – Coordinate frame (default: ‘icrf’). - ‘icrf’: Moon position relative to Solar System Barycenter - ‘ecliptic’: Ecliptic coordinates - ‘earth_centered’: Position relative to Earth

Returns:

• position (ndarray, shape (3,)) – Moon position in AU (or relative to Earth for ‘earth_centered’)

• velocity (ndarray, shape (3,)) – Moon velocity in AU/day

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Notes

By default, returns Moon position relative to the Solar System Barycenter. Use frame=’earth_centered’ for geocentric coordinates.

Examples

>>> eph = DEEphemeris()
>>> r, v = eph.moon_position(2451545.0, frame='earth_centered')

planet_position(planet, jd, frame='icrf')

Compute planet position and velocity.

Parameters:

• planet (str) – Planet name: ‘mercury’, ‘venus’, ‘mars’, ‘jupiter’, ‘saturn’, ‘uranus’, ‘neptune’

• jd (float) – Julian Day in Terrestrial Time (TT)

• frame ({'icrf', 'ecliptic'}, optional) – Coordinate frame (default: ‘icrf’)

Returns:

• position (ndarray, shape (3,)) – Planet position in AU

• velocity (ndarray, shape (3,)) – Planet velocity in AU/day

Raises:

ValueError – If planet name is not recognized

Return type:

Tuple[ndarray[Any, Any], ndarray[Any, Any]]

Examples

>>> eph = DEEphemeris()
>>> r, v = eph.planet_position('mars', 2451545.0)

barycenter_position(body, jd)

Compute position of any body relative to Solar System Barycenter.

Parameters:

• body (str) – Body name (‘sun’, ‘moon’, ‘mercury’, …, ‘neptune’)

• jd (float) – Julian Day in Terrestrial Time (TT)

Returns:

• position (ndarray, shape (3,)) – Position in AU

• velocity (ndarray, shape (3,)) – Velocity in AU/day

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

clear_cache()

Clear internal position cache.

pytcl.astronomical.ephemerides.sun_position(jd, frame='icrf')

Convenience function: Compute Sun position and velocity.

Parameters:

• jd (float) – Julian Day in Terrestrial Time (TT)

• frame ({'icrf', 'ecliptic'}, optional) – Coordinate frame (default: ‘icrf’)

Returns:

• position (ndarray, shape (3,)) – Sun position in AU

• velocity (ndarray, shape (3,)) – Sun velocity in AU/day

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Examples

>>> r, v = sun_position(2451545.0)  # J2000.0
>>> np.linalg.norm(r)  # Distance from SSB
0.00...

See also

DEEphemeris.sun_position

Full ephemeris class with caching

pytcl.astronomical.ephemerides.moon_position(jd, frame='icrf')

Convenience function: Compute Moon position and velocity.

Parameters:

• jd (float) – Julian Day in Terrestrial Time (TT)

• frame ({'icrf', 'ecliptic', 'earth_centered'}, optional) – Coordinate frame (default: ‘icrf’)

Returns:

• position (ndarray, shape (3,)) – Moon position in AU

• velocity (ndarray, shape (3,)) – Moon velocity in AU/day

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Examples

>>> r, v = moon_position(2451545.0, 'earth_centered')
>>> np.linalg.norm(r) * 149597870.7  # Distance in km
402...

See also

DEEphemeris.moon_position

Full ephemeris class with caching

pytcl.astronomical.ephemerides.planet_position(planet, jd, frame='icrf')

Convenience function: Compute planet position and velocity.

Parameters:

• planet (str) – Planet name

• jd (float) – Julian Day in Terrestrial Time (TT)

• frame ({'icrf', 'ecliptic'}, optional) – Coordinate frame (default: ‘icrf’)

Returns:

• position (ndarray, shape (3,)) – Planet position in AU

• velocity (ndarray, shape (3,)) – Planet velocity in AU/day

Return type:

Tuple[ndarray[Any, Any], ndarray[Any, Any]]

See also

DEEphemeris.planet_position

Full ephemeris class with caching

pytcl.astronomical.ephemerides.barycenter_position(body, jd)

Convenience function: Position relative to Solar System Barycenter.

Parameters:

• body (str) – Body name (‘sun’, ‘moon’, ‘mercury’, …, ‘neptune’)

• jd (float) – Julian Day in Terrestrial Time (TT)

Returns:

• position (ndarray, shape (3,)) – Position in AU

• velocity (ndarray, shape (3,)) – Velocity in AU/day

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Examples

>>> r, v = barycenter_position('mars', 2451545.0)
>>> np.linalg.norm(r)  # Mars distance from SSB
1.4...

Relativistic Corrections

Relativistic effects in orbital mechanics and space-time geometry.

Relativistic corrections for precision astronomy and satellite positioning.

This module provides utilities for computing relativistic effects in orbital mechanics, including gravitational time dilation, Shapiro delay, and coordinate transformations in the Schwarzschild metric. These effects are critical for high-precision applications such as GPS, pulsar timing, and celestial mechanics.

Key Physical Constants:

• Schwarzschild radius: r_s = 2GM/c^2

• Gravitational parameter for Earth: μ = GM = 3.986004418e14 m^3/s^2

• Speed of light: c = 299792458 m/s

• Gravitational constant: G = 6.67430e-11 m^3/(kg·s^2)

References:

• Soffel et al. (2003): The IAU 2000 Resolutions for Astrometry, Celestial Mechanics, and Reference Frames

• Will, C. M. (2014): The Confrontation between General Relativity and Experiment

• Ries et al. (1992): Preliminary Analysis of LAGEOS II Observations for the Determination of Relativistic Effects

pytcl.astronomical.relativity.schwarzschild_radius(mass)

Compute Schwarzschild radius for a given mass.

The Schwarzschild radius is the radius at which an object becomes a black hole. It is given by r_s = 2GM/c^2.

Parameters:

mass (float) – Mass of the object (kg)

Returns:

Schwarzschild radius (m)

Return type:

float

Examples

>>> r_s_earth = schwarzschild_radius(5.972e24)  # Earth's mass
>>> print(f"Earth's Schwarzschild radius: {r_s_earth:.3e} m")
Earth's Schwarzschild radius: 8.870e-03 m

>>> r_s_sun = schwarzschild_radius(1.989e30)  # Sun's mass
>>> print(f"Sun's Schwarzschild radius: {r_s_sun:.3e} m")
Sun's Schwarzschild radius: 2.952e+03 m

pytcl.astronomical.relativity.gravitational_time_dilation(r, gm=398600441800000.0)

Compute gravitational time dilation factor sqrt(1 - 2GM/(rc^2)).

In general relativity, time passes slower in stronger gravitational fields. This function computes the metric coefficient g_00 for the Schwarzschild metric, which determines proper time relative to coordinate time.

Parameters:

• r (float) – Distance from the gravitational body (m)

• gm (float, optional) – Gravitational parameter GM of the body (m^3/s^2). Default is GM_EARTH.

Returns:

Time dilation factor in [0, 1]. A value less than 1 indicates time passes slower at radius r compared to infinity.

Return type:

float

Raises:

ValueError – If r is less than or equal to Schwarzschild radius

Examples

Compute time dilation at Earth’s surface (6371 km):

>>> r_earth = 6.371e6  # meters
>>> dilation = gravitational_time_dilation(r_earth)
>>> print(f"Time dilation at surface: {dilation:.15f}")
Time dilation at surface: 0.999999999300693

At GPS orbital altitude (~20,200 km):

>>> r_gps = 26.56e6  # meters
>>> dilation_gps = gravitational_time_dilation(r_gps)
>>> time_shift = (1 - dilation_gps) * 86400 * 1e9  # nanoseconds per day
>>> print(f"Time shift: {time_shift:.1f} ns/day")

pytcl.astronomical.relativity.proper_time_rate(v, r, gm=398600441800000.0)

Compute proper time rate accounting for both velocity and gravity.

The proper time rate combines special relativistic time dilation from velocity and general relativistic time dilation from the gravitational potential.

d(tau)/d(t) = sqrt(1 - v^2/c^2) * sqrt(1 - 2GM/(rc^2))

For small velocities and weak fields: 1 - v^2/(2c^2) - GM/(rc^2)

Parameters:

• v (float) – Velocity magnitude (m/s)

• r (float) – Distance from gravitational body (m)

• gm (float, optional) – Gravitational parameter GM (m^3/s^2). Default is GM_EARTH.

Returns:

Proper time rate. A value less than 1 indicates proper time passes slower than coordinate time.

Return type:

float

Examples

Proper time rate for a GPS satellite at ~3.87 km/s and 26.56 Mm altitude:

>>> v_gps = 3870.0  # m/s
>>> r_gps = 26.56e6  # m
>>> rate = proper_time_rate(v_gps, r_gps)
>>> print(f"Proper time rate: {rate:.15f}")
>>> time_shift = (1 - rate) * 86400  # seconds per day
>>> print(f"Daily time shift: {time_shift:.3f} s/day")

pytcl.astronomical.relativity.shapiro_delay(observer_pos, light_source_pos, gravitating_body_pos, gm=1.32712440018e+20)

Compute Shapiro time delay for light propagation through gravitational field.

The Shapiro delay is the additional propagation time experienced by light traveling through a gravitational field, compared to flat spacetime.

delay = (2GM/c^3) * ln((r_o + r_s + r_os) / (r_o + r_s - r_os))

where r_o is distance from body to observer, r_s is distance from body to source, and r_os is distance from observer to source.

Parameters:

• observer_pos (np.ndarray) – Position of observer (m), shape (3,)

• light_source_pos (np.ndarray) – Position of light source (m), shape (3,)

• gravitating_body_pos (np.ndarray) – Position of gravitating body (m), shape (3,)

• gm (float, optional) – Gravitational parameter GM (m^3/s^2). Default is GM_SUN.

Returns:

Shapiro delay (seconds)

Return type:

float

Examples

Earth-Sun-Spacecraft signal at superior conjunction (worst case):

>>> # Simplified geometry: Sun at origin, Earth at 1 AU, spacecraft beyond at distance
>>> sun_pos = np.array([0.0, 0.0, 0.0])
>>> earth_pos = np.array([1.496e11, 0.0, 0.0])  # 1 AU
>>> spacecraft_pos = np.array([1.496e11, 1.0e11, 0.0])  # Far from sun
>>> delay = shapiro_delay(earth_pos, spacecraft_pos, sun_pos, GM_SUN)
>>> print(f"Shapiro delay: {delay:.3e} seconds")
>>> print(f"Shapiro delay: {delay*1e6:.1f} microseconds")

pytcl.astronomical.relativity.schwarzschild_precession_per_orbit(a, e, gm=1.32712440018e+20)

Compute perihelion precession per orbit due to general relativity.

The advance of perihelion for an orbit around a central mass M is:

Δφ = (6π * GM) / (c^2 * a * (1 - e^2))

This effect is a key test of general relativity. For Mercury, the predicted precession is ~43 arcseconds per century.

Parameters:

• a (float) – Semi-major axis (m)

• e (float) – Eccentricity (dimensionless), must be in [0, 1)

• gm (float, optional) – Gravitational parameter GM (m^3/s^2). Default is GM_SUN.

Returns:

Perihelion precession per orbit (radians)

Return type:

float

Examples

Mercury’s perihelion precession (GR contribution):

>>> a_mercury = 0.38709927 * AU  # Semi-major axis in meters
>>> e_mercury = 0.20563593     # Eccentricity
>>> precession_rad = schwarzschild_precession_per_orbit(a_mercury, e_mercury, GM_SUN)
>>> precession_arcsec = precession_rad * 206265  # Convert to arcseconds
>>> orbital_period = 87.969  # days
>>> centuries = 36525 / orbital_period  # Orbits per century
>>> precession_per_century = precession_arcsec * centuries
>>> print(f"GR perihelion precession: {precession_per_century:.1f} arcsec/century")
GR perihelion precession: 42.98 arcsec/century

pytcl.astronomical.relativity.post_newtonian_acceleration(r_vec, v_vec, gm=398600441800000.0)

Compute post-Newtonian acceleration corrections (1PN order).

Extends Newtonian gravity with first-order post-Newtonian corrections.

a_PN = -GM/r^2 * u_r + a_1PN

where a_1PN includes velocity-dependent and metric perturbation terms.

Parameters:

• r_vec (np.ndarray) – Position vector (m), shape (3,)

• v_vec (np.ndarray) – Velocity vector (m/s), shape (3,)

• gm (float, optional) – Gravitational parameter GM (m^3/s^2). Default is GM_EARTH.

Returns:

Total acceleration including 1PN corrections (m/s^2), shape (3,)

Return type:

np.ndarray

Examples

Compare Newtonian and PN acceleration for LEO satellite:

>>> r = np.array([6.678e6, 0.0, 0.0])  # ~300 km altitude
>>> v = np.array([0.0, 7.7e3, 0.0])    # Circular orbit velocity
>>> a_total = post_newtonian_acceleration(r, v)
>>> a_newt = -GM_EARTH / np.linalg.norm(r)**3 * r
>>> correction_ratio = np.linalg.norm(a_total - a_newt) / np.linalg.norm(a_newt)
>>> print(f"PN correction: {correction_ratio*1e6:.1f} ppm")

pytcl.astronomical.relativity.geodetic_precession(a, e, inclination, gm=398600441800000.0)

Compute geodetic (de Sitter) precession rate of orbital plane.

The orbital plane of a satellite precesses due to frame-dragging effects and spacetime curvature. The geodetic precession rate is:

Ω_geodetic = -GM/(c^2 * a^3 * (1 - e^2)^2) * cos(i)

Parameters:

• a (float) – Semi-major axis (m)

• e (float) – Eccentricity (dimensionless)

• inclination (float) – Orbital inclination (radians)

• gm (float, optional) – Gravitational parameter (m^3/s^2). Default is GM_EARTH.

Returns:

Geodetic precession rate (radians per orbit)

Return type:

float

Examples

Geodetic precession for a typical Earth satellite:

>>> a = 6.678e6  # ~300 km altitude
>>> e = 0.0      # Circular
>>> i = np.radians(51.6)  # ISS-like inclination
>>> rate = geodetic_precession(a, e, i)
>>> print(f"Precession per orbit: {rate*206265:.3f} arcsec")

pytcl.astronomical.relativity.lense_thirring_precession(a, e, inclination, angular_momentum, gm=398600441800000.0)

Compute Lense-Thirring (frame-dragging) precession of orbital node.

A rotating central body drags the orbital plane of nearby objects. The nodal precession rate due to this effect is:

Ω_LT = (2GM * J_2 * a * ω) / (c^2 * p^2) * f(e, i)

where J_2 is the quadrupole moment, ω is angular velocity, and f depends on eccentricity and inclination.

Parameters:

• a (float) – Semi-major axis (m)

• e (float) – Eccentricity (dimensionless)

• inclination (float) – Orbital inclination (radians)

• angular_momentum (float) – Angular momentum of central body (kg·m^2/s)

• gm (float, optional) – Gravitational parameter (m^3/s^2). Default is GM_EARTH.

Returns:

Lense-Thirring precession rate (radians per orbit)

Return type:

float

Notes

This is a simplified version. For Earth, J_2 effects typically dominate classical nodal precession, while Lense-Thirring is a small correction (~1 arcsec per year for typical satellites).

Examples

Lense-Thirring effect for LAGEOS satellite:

>>> # LAGEOS parameters
>>> a = 12.27e6  # Semi-major axis
>>> e = 0.0045
>>> i = np.radians(109.9)
>>> L_earth = 7.05e33  # Earth's angular momentum
>>> rate = lense_thirring_precession(a, e, i, L_earth)
>>> print(f"LT precession per orbit: {rate*206265*1e3:.1f} milliarcsec")

pytcl.astronomical.relativity.relativistic_range_correction(distance, relative_velocity, gm=398600441800000.0)

Compute relativistic range correction for ranging measurements.

When measuring distance to a satellite or spacecraft using ranging (e.g., laser ranging), relativistic effects introduce corrections to the measured range.

The main contributions are: - Gravitational time dilation - Relativistic Doppler effect

Parameters:

• distance (float) – Distance to object (m)

• relative_velocity (float) – Radial velocity component (m/s, positive = receding)

• gm (float, optional) – Gravitational parameter (m^3/s^2). Default is GM_EARTH.

Returns:

Range correction (m)

Return type:

float

Examples

Range correction for lunar laser ranging:

>>> distance_to_moon = 3.84e8  # meters
>>> radial_velocity = 0.0  # Average over orbit
>>> correction = relativistic_range_correction(distance_to_moon, radial_velocity, GM_EARTH)
>>> print(f"Range correction: {correction:.1f} m")