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
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
Vallado, D. A., “Fundamentals of Astrodynamics and Applications,” 4th ed., Microcosm Press, 2013.
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)[source]
Bases:
NamedTupleClassical (Keplerian) orbital elements.
- Variables:
a (float) – Semi-major axis (km). Negative for hyperbolic orbits.
e (float) – Eccentricity. 0=circular, 0<e<1=elliptic, e=1=parabolic, e>1=hyperbolic.
i (float) – Inclination (radians), 0 to pi.
raan (float) – Right ascension of ascending node (radians), 0 to 2*pi.
omega (float) – Argument of periapsis (radians), 0 to 2*pi.
nu (float) – True anomaly (radians), 0 to 2*pi.
- class pytcl.astronomical.orbital_mechanics.StateVector(r, v)[source]
Bases:
NamedTupleCartesian state vector.
- Variables:
r (ndarray) – Position vector (km), shape (3,).
v (ndarray) – Velocity vector (km/s), shape (3,).
- pytcl.astronomical.orbital_mechanics.mean_to_eccentric_anomaly(M, e, tol=1e-12, max_iter=50)[source]
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:
- Returns:
E – Eccentric anomaly (radians).
- Return type:
- 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") E = 72.2903 degrees
- pytcl.astronomical.orbital_mechanics.mean_to_hyperbolic_anomaly(M, e, tol=1e-12, max_iter=50)[source]
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:
- Returns:
H – Hyperbolic anomaly (radians).
- Return type:
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)[source]
Convert eccentric anomaly to true anomaly.
- Parameters:
- Returns:
nu – True anomaly (radians), in [0, 2*pi).
- Return type:
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)[source]
Convert true anomaly to eccentric anomaly.
- Parameters:
- Returns:
E – Eccentric anomaly (radians), in [0, 2*pi).
- Return type:
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)[source]
Convert hyperbolic anomaly to true anomaly.
- Parameters:
- Returns:
nu – True anomaly (radians).
- Return type:
Examples
>>> nu = hyperbolic_to_true_anomaly(0.5, 1.5) >>> isinstance(nu, float) True
- pytcl.astronomical.orbital_mechanics.true_to_hyperbolic_anomaly(nu, e)[source]
Convert true anomaly to hyperbolic anomaly.
- pytcl.astronomical.orbital_mechanics.eccentric_to_mean_anomaly(E, e)[source]
Convert eccentric anomaly to mean anomaly (Kepler’s equation).
- Parameters:
- Returns:
M – Mean anomaly (radians).
- Return type:
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)[source]
Convert mean anomaly to true anomaly.
- Parameters:
- Returns:
nu – True anomaly (radians).
- Return type:
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)[source]
Convert true anomaly to mean anomaly.
- pytcl.astronomical.orbital_mechanics.orbital_elements_to_state(elements, mu=398600.4418)[source]
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:
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}") # Perigee radius a*(1-e) = 6930 km along x r = [6930. 0. 0.]
- pytcl.astronomical.orbital_mechanics.state_to_orbital_elements(state, mu=398600.4418)[source]
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:
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)[source]
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:
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)[source]
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:
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)[source]
Compute orbital period for elliptic orbit.
- Parameters:
- Returns:
T – Orbital period (seconds).
- Return type:
Examples
>>> T = orbital_period(7000) # LEO satellite >>> T / 60 # Convert to minutes 97.8...
- pytcl.astronomical.orbital_mechanics.mean_motion(a, mu=398600.4418)[source]
Compute mean motion.
- Parameters:
- Returns:
n – Mean motion (radians/second).
- Return type:
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)[source]
Compute orbital velocity using vis-viva equation.
- Parameters:
- Returns:
v – Orbital velocity (km/s).
- Return type:
Examples
>>> v = vis_viva(7000, 7000) # Circular orbit >>> abs(v - circular_velocity(7000)) < 0.01 True
- pytcl.astronomical.orbital_mechanics.specific_angular_momentum(state)[source]
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)[source]
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:
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)[source]
Compute flight path angle.
- Parameters:
state (StateVector) – State vector.
- Returns:
gamma – Flight path angle (radians). Positive when climbing, negative when descending.
- Return type:
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)[source]
Compute periapsis radius.
- Parameters:
- Returns:
r_p – Periapsis radius (km).
- Return type:
Examples
>>> r_p = periapsis_radius(10000, 0.3) >>> r_p 7000.0
- pytcl.astronomical.orbital_mechanics.apoapsis_radius(a, e)[source]
Compute apoapsis radius.
- Parameters:
- Returns:
r_a – Apoapsis radius (km). Infinite for parabolic/hyperbolic orbits.
- Return type:
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)[source]
Compute time since periapsis passage.
- Parameters:
- Returns:
t – Time since periapsis (seconds).
- Return type:
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)[source]
Compute orbital radius at given true anomaly.
- Parameters:
- Returns:
r – Orbital radius (km).
- Return type:
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)[source]
Compute escape velocity at given radius.
- Parameters:
- Returns:
v_esc – Escape velocity (km/s).
- Return type:
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)[source]
Compute circular orbital velocity at given radius.
- Parameters:
- Returns:
v_circ – Circular velocity (km/s).
- Return type:
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
Vallado, D. A., “Fundamentals of Astrodynamics and Applications,” 4th ed., Microcosm Press, 2013.
Izzo, D., “Revisiting Lambert’s problem,” Celestial Mechanics and Dynamical Astronomy, 2015.
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)[source]
Bases:
NamedTupleSolution to Lambert’s problem.
- Variables:
- pytcl.astronomical.lambert.lambert_universal(r1, r2, tof, mu=398600.4418, prograde=True, low_path=True, max_iter=100, tol=1e-10)[source]
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:
- 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) >>> bool(np.allclose(sol.v1, [-5.9925, 1.9254, 3.2456], atol=1e-3)) # Curtis Ex. 5.2 True
- pytcl.astronomical.lambert.lambert_izzo(r1, r2, tof, mu=398600.4418, prograde=True, multi_rev=0, max_iter=100, tol=1e-10)[source]
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:
Notes
For multi-revolution transfers, there may be two solutions (low and high energy). This returns the low energy solution.
- Raises:
ValueError – If the iteration does not converge (e.g. the requested time of flight is below the minimum for the requested number of revolutions).
- pytcl.astronomical.lambert.minimum_energy_transfer(r1, r2, mu=398600.4418, prograde=True)[source]
Compute minimum energy transfer between two positions.
- Parameters:
- Returns:
tof_min (float) – Minimum energy time of flight (seconds).
solution (LambertSolution) – Lambert solution at minimum energy.
- Return type:
- pytcl.astronomical.lambert.hohmann_transfer(r1, r2, mu=398600.4418)[source]
Compute Hohmann transfer between two circular orbits.
- Parameters:
- 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:
Examples
>>> dv1, dv2, tof = hohmann_transfer(6678, 42164) # LEO to GEO >>> print(f"Total dv = {dv1 + dv2:.3f} km/s") Total dv = 3.893 km/s
- pytcl.astronomical.lambert.bi_elliptic_transfer(r1, r2, r_intermediate, mu=398600.4418)[source]
Compute bi-elliptic transfer between two circular orbits.
- Parameters:
- 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:
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
Vallado, D. A., “Fundamentals of Astrodynamics and Applications,” 4th ed., Microcosm Press, 2013.
IERS Conventions (2010), IERS Technical Note No. 36.
Capitaine et al., “Expressions for IAU 2000 precession quantities,” A&A, 2003.
- pytcl.astronomical.reference_frames.julian_centuries_j2000(jd)[source]
Compute Julian centuries since J2000.0.
- pytcl.astronomical.reference_frames.precession_angles_iau76(T)[source]
Compute IAU 1976 precession angles.
- pytcl.astronomical.reference_frames.precession_matrix_iau76(jd)[source]
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)[source]
Compute IAU 1980 nutation angles (simplified).
- pytcl.astronomical.reference_frames.nutation_matrix(jd)[source]
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)[source]
Compute mean obliquity of the ecliptic (IAU 1980).
- pytcl.astronomical.reference_frames.true_obliquity(jd_tt)[source]
Compute true obliquity of the ecliptic.
- pytcl.astronomical.reference_frames.earth_rotation_angle(jd_ut1)[source]
Compute Earth Rotation Angle (ERA).
- pytcl.astronomical.reference_frames.gmst_iau82(jd_ut1)[source]
Compute Greenwich Mean Sidereal Time (GMST) using IAU 1982 model.
- pytcl.astronomical.reference_frames.gast_iau82(jd_ut1, jd_tt)[source]
Compute Greenwich Apparent Sidereal Time (GAST).
- pytcl.astronomical.reference_frames.sidereal_rotation_matrix(theta)[source]
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)[source]
Compute equation of the equinoxes.
- pytcl.astronomical.reference_frames.polar_motion_matrix(xp, yp)[source]
Compute polar motion rotation matrix.
- pytcl.astronomical.reference_frames.gcrf_to_itrf(r_gcrf, jd_ut1, jd_tt, xp=0.0, yp=0.0)[source]
Transform position from GCRF (inertial) to ITRF (Earth-fixed).
- Parameters:
- 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)[source]
Transform position from ITRF (Earth-fixed) to GCRF (inertial).
- Parameters:
- 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)[source]
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:
- 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_gcrfInverse transformation
gcrf_to_itrfIncludes polar motion
References
Vallado et al., “Fundamentals of Astrodynamics and Applications”, 4th ed.
- pytcl.astronomical.reference_frames.pef_to_gcrf(r_pef, jd_ut1, jd_tt)[source]
Transform position from PEF (Earth-fixed, rotation only) to GCRF (inertial).
Inverse of gcrf_to_pef.
- Parameters:
- Returns:
r_gcrf – Position in GCRF (km), shape (3,).
- Return type:
ndarray
See also
gcrf_to_pefForward transformation
- pytcl.astronomical.reference_frames.eci_to_ecef(r_eci, gmst)[source]
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)[source]
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)[source]
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)[source]
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)[source]
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)[source]
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)[source]
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:
- 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
>>> r_teme = np.array([7000.0, 0.0, 0.0]) # km >>> r_itrf = teme_to_itrf(r_teme, 2460311.0) >>> bool(np.isclose(np.linalg.norm(r_itrf), 7000.0)) # Rotation preserves norm True
Typical use with SGP4 output (requires a parsed TLE):
>>> 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)[source]
Transform position from ITRF to TEME.
- Parameters:
- Returns:
r_teme – Position in TEME frame (km), shape (3,).
- Return type:
ndarray
- pytcl.astronomical.reference_frames.teme_to_gcrf(r_teme, jd_tt)[source]
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
>>> r_teme = np.array([7000.0, 0.0, 0.0]) # km >>> r_gcrf = teme_to_gcrf(r_teme, 2460311.0) >>> bool(np.isclose(np.linalg.norm(r_gcrf), 7000.0)) # Rotation preserves norm True
Typical use with SGP4 output (requires a parsed TLE):
>>> 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)[source]
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)[source]
Transform position and velocity from TEME to ITRF.
This properly accounts for the velocity transformation including the Earth’s rotation rate.
- pytcl.astronomical.reference_frames.itrf_to_teme_with_velocity(r_itrf, v_itrf, jd_ut1, xp=0.0, yp=0.0)[source]
Transform position and velocity from ITRF to TEME.
- pytcl.astronomical.reference_frames.gcrf_to_mod(r_gcrf, jd_tt)[source]
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_gcrfInverse transformation.
gcrf_to_todIncludes nutation for true of date.
- pytcl.astronomical.reference_frames.mod_to_gcrf(r_mod, jd_tt)[source]
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_modForward transformation.
- pytcl.astronomical.reference_frames.gcrf_to_tod(r_gcrf, jd_tt)[source]
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_gcrfInverse transformation.
gcrf_to_modMean of date (without nutation).
- pytcl.astronomical.reference_frames.tod_to_gcrf(r_tod, jd_tt)[source]
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_todForward transformation.
- pytcl.astronomical.reference_frames.mod_to_tod(r_mod, jd_tt)[source]
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_modInverse transformation.
- pytcl.astronomical.reference_frames.tod_to_mod(r_tod, jd_tt)[source]
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_todForward transformation.
- pytcl.astronomical.reference_frames.tod_to_itrf(r_tod, jd_ut1, jd_tt=None, xp=0.0, yp=0.0)[source]
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_todInverse transformation.
- pytcl.astronomical.reference_frames.itrf_to_tod(r_itrf, jd_ut1, jd_tt=None, xp=0.0, yp=0.0)[source]
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_itrfForward transformation.
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)[source]
Convert calendar date to Julian Date.
- Parameters:
- Returns:
Julian Date.
- Return type:
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)[source]
Convert Julian Date to calendar date.
- Parameters:
jd (float) – Julian Date.
- Returns:
(year, month, day, hour, minute, second).
- Return type:
Examples
>>> jd_to_cal(2451545.0) # J2000.0 (2000, 1, 1, 12, 0, 0.0)
- pytcl.astronomical.time_systems.mjd_to_jd(mjd)[source]
Convert Modified Julian Date to Julian Date.
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)[source]
Convert Julian Date to Modified Julian Date.
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)[source]
Convert UTC to TAI (Julian Date).
- Parameters:
- Returns:
TAI as Julian Date.
- Return type:
Notes
TAI = UTC + leap_seconds
- pytcl.astronomical.time_systems.tai_to_utc(jd_tai)[source]
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:
Notes
This is an approximate conversion that may have small errors near leap second boundaries.
Notes
The leap-second lookup is by table, so the conversion is exact except within one second of an insertion. Instants inside a leap second itself – 23:59:60 on an insertion date – have no distinct representation here and are attributed to the following second (gh-25). Sub-second work spanning a leap-second boundary needs a library that models UTC as a discontinuous scale, such as astropy.
- pytcl.astronomical.time_systems.tai_to_tt(jd_tai)[source]
Convert TAI (Julian Date) to TT (Terrestrial Time).
Notes
TT = TAI + 32.184 seconds (constant offset)
- pytcl.astronomical.time_systems.tt_to_tai(jd_tt)[source]
Convert TT (Terrestrial Time) to TAI (Julian Date).
- pytcl.astronomical.time_systems.utc_to_tt(year, month, day, hour=0, minute=0, second=0.0)[source]
Convert UTC to TT (Terrestrial Time).
- Parameters:
- Returns:
TT as Julian Date.
- Return type:
Notes
TT = UTC + leap_seconds + 32.184
- pytcl.astronomical.time_systems.tai_to_gps(jd_tai)[source]
Convert TAI to GPS time.
- Parameters:
jd_tai (float) – TAI as Julian Date.
- Returns:
GPS time as Julian Date.
- Return type:
Notes
GPS time = TAI - 19 seconds (offset at GPS epoch) GPS time does not include leap seconds after 1980.
- pytcl.astronomical.time_systems.utc_to_gps(year, month, day, hour=0, minute=0, second=0.0)[source]
Convert UTC to GPS time.
- pytcl.astronomical.time_systems.unix_to_jd(unix_time)[source]
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:
Examples
>>> unix_to_jd(0.0) # Unix epoch 2440587.5
- pytcl.astronomical.time_systems.jd_to_unix(jd)[source]
Convert Julian Date to Unix timestamp.
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)[source]
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:
Examples
>>> gps_week_seconds(JD_GPS_EPOCH) (0, 0.0)
- pytcl.astronomical.time_systems.gps_week_to_utc(week, seconds)[source]
Convert GPS week and seconds to UTC.
- pytcl.astronomical.time_systems.gmst(jd_ut1)[source]
Compute Greenwich Mean Sidereal Time.
Notes
Uses the IAU 1982 expression for GMST.
References
Explanatory Supplement to the Astronomical Almanac, 3rd ed.
- pytcl.astronomical.time_systems.gast(jd_ut1, dpsi=0.0, eps=0.0)[source]
Compute Greenwich Apparent Sidereal Time.
- Parameters:
- Returns:
GAST in radians.
- Return type:
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.
Notes
With the default ``dpsi=0`` this returns GMST, not GAST. The equation of the equinoxes is
dpsi * cos(eps), so leaving both at zero removes it entirely and the two are identical to the last bit (gh-25). Supply the nutation in longitude and the obliquity – from ~pytcl.astronomical.reference_frames.nutation_angles_iau80 – to get apparent sidereal time. The difference is up to about 1.1 arcseconds, which is 0.02 arcseconds of longitude error at the equator per milliarcsecond of neglected nutation.
- pytcl.astronomical.time_systems.get_leap_seconds(year, month, day)[source]
Get the number of leap seconds (TAI-UTC offset) for a given date.
- Parameters:
- Returns:
TAI-UTC offset in seconds.
- Return type:
Examples
>>> get_leap_seconds(2020, 1, 1) 37 >>> get_leap_seconds(1980, 1, 6) 19
- class pytcl.astronomical.time_systems.LeapSecondTable[source]
Bases:
objectTable 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.
Notes
This table must be updated when new leap seconds are announced. Last update: 2017-01-01 (37 seconds).
No leap second has been added since 2017; the table is current as of the last library release.
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
>>>
>>> # Load ephemeris (auto-downloads if needed)
>>> eph = DEEphemeris(version='DE440')
>>>
>>> # Query Sun position relative to the Solar System Barycenter (AU).
>>> # The Sun orbits the SSB within ~1-2 solar radii, so this is small.
>>> jd = 2451545.0 # J2000.0
>>> r_sun, v_sun = eph.sun_position(jd)
>>> print(f"Sun distance from SSB: {np.linalg.norm(r_sun):.6f} AU")
Sun distance from SSB: 0.007668 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, relative to the Solar System Barycenter (SSB) unless stated otherwise
For highest precision, use DE440 (latest release) or DE430 (2013); DE432s trades span for size (1950-2050 only, not higher precision)
References
Standish, E. M. (1995). “Report of the IAU WGAS Sub-group on Numerical Standards”. In Highlights of Astronomy (Vol. 10).
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')[source]
Bases:
objectHigh-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: Short-span version (2013) – the “s” denotes the reduced time range, not extra precision; covers 1950-2050 per NAIF’s generic_kernels/spk/planets/aareadme.txt (full-span DE432 covers 1550-2650)
DE430: Earlier release (2013), covers 1550-2650
DE405: Older version (1998), compact, covers 1600-2200
- Variables:
- 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')[source]
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')[source]
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 from SSB: {np.linalg.norm(r):.6f} AU") Distance from SSB: 0.007668 AU
- moon_position(jd, frame='icrf')[source]
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')[source]
Compute planet position and velocity.
- Parameters:
- 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:
Examples
>>> eph = DEEphemeris() >>> r, v = eph.planet_position('mars', 2451545.0)
- barycenter_position(body, jd)[source]
Compute position of any body relative to Solar System Barycenter.
- pytcl.astronomical.ephemerides.sun_position(jd, frame='icrf')[source]
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_positionFull ephemeris class with caching
- pytcl.astronomical.ephemerides.moon_position(jd, frame='icrf')[source]
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_positionFull ephemeris class with caching
- pytcl.astronomical.ephemerides.planet_position(planet, jd, frame='icrf')[source]
Convenience function: Compute planet position and velocity.
- Parameters:
- Returns:
position (ndarray, shape (3,)) – Planet position in AU
velocity (ndarray, shape (3,)) – Planet velocity in AU/day
- Return type:
See also
DEEphemeris.planet_positionFull ephemeris class with caching
- pytcl.astronomical.ephemerides.barycenter_position(body, jd)[source]
Convenience function: Position relative to Solar System Barycenter.
- Parameters:
- 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)[source]
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:
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.954e+03 m
- pytcl.astronomical.relativity.gravitational_time_dilation(r, gm=398600441800000.0)[source]
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:
- Returns:
Time dilation factor in [0, 1]. A value less than 1 indicates time passes slower at radius r compared to infinity.
- Return type:
- 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.999999999303873
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") Time shift: 14427.2 ns/day
- pytcl.astronomical.relativity.proper_time_rate(v, r, gm=398600441800000.0)[source]
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:
- Returns:
Proper time rate. A value less than 1 indicates proper time passes slower than coordinate time.
- Return type:
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}") Proper time rate: 0.999999999749698 >>> time_shift = (1 - rate) * 86400 * 1e6 # microseconds per day >>> print(f"Daily time shift: {time_shift:.1f} us/day") Daily time shift: 21.6 us/day
- pytcl.astronomical.relativity.shapiro_delay(observer_pos, light_source_pos, gravitating_body_pos, gm=1.32712440018e+20)[source]
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:
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") Shapiro delay: 6.173e-06 seconds >>> print(f"Shapiro delay: {delay*1e6:.1f} microseconds") Shapiro delay: 6.2 microseconds
- pytcl.astronomical.relativity.schwarzschild_precession_per_orbit(a, e, gm=1.32712440018e+20)[source]
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:
- Returns:
Perihelion precession per orbit (radians)
- Return type:
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:.2f} arcsec/century") GR perihelion precession: 42.98 arcsec/century
- pytcl.astronomical.relativity.post_newtonian_acceleration(r_vec, v_vec, gm=398600441800000.0)[source]
Compute total acceleration including first post-Newtonian corrections.
Extends Newtonian gravity with the first post-Newtonian (Schwarzschild) correction in harmonic coordinates:
a_1PN = (GM / (c^2 * r^3)) * [(4*GM/r - v.v) * r_vec + 4*(r_vec . v_vec) * v_vec]
- 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 a circular 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 >>> ratio = np.linalg.norm(a_total - a_newt) / np.linalg.norm(a_newt) >>> print(f"PN correction ratio: {ratio:.3e}") PN correction ratio: 1.997e-09
- pytcl.astronomical.relativity.geodetic_precession(a, e, inclination, gm=398600441800000.0)[source]
Compute geodetic (de Sitter) precession per orbit.
A gyroscope (or the orbital angular momentum vector of a satellite) transported around a central mass precesses due to spacetime curvature. The geodetic precession accumulated over one orbit is:
Δω = 3π * GM / (c^2 * a * (1 - e^2))
The result is positive (prograde) and its magnitude is independent of the orbital inclination.
- Parameters:
- Returns:
Geodetic precession per orbit (radians), positive (prograde).
- Return type:
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:.3e} rad") Precession per orbit: 6.259e-09 rad
De Sitter precession of the Earth-Moon gyroscope orbiting the Sun:
>>> rate = geodetic_precession(AU, 0.0167, 0.0, GM_SUN) >>> per_century = rate * 100 * 206264.806 # ~100 orbits per century >>> print(f"De Sitter precession: {per_century:.2f} arcsec/century") De Sitter precession: 1.92 arcsec/century
- pytcl.astronomical.relativity.lense_thirring_precession(a, e, inclination, angular_momentum, gm=398600441800000.0)[source]
Compute Lense-Thirring (frame-dragging) nodal precession rate.
A rotating central body drags the orbital plane of nearby objects. The precession rate of the ascending node due to this effect is:
Ω̇_LT = 2*G*J / (c^2 * a^3 * (1 - e^2)^(3/2))
where J is the spin angular momentum of the central body. The rate is independent of the orbital inclination. Multiply by the orbital period to obtain the precession per orbit.
- Parameters:
a (float) – Semi-major axis (m)
e (float) – Eccentricity (dimensionless)
inclination (float) – Orbital inclination (radians). Retained for API compatibility; the nodal Lense-Thirring rate does not depend on it.
angular_momentum (float) – Spin angular momentum J of central body (kg·m^2/s). For Earth, J ≈ 5.86e33 kg·m^2/s.
gm (float, optional) – Gravitational parameter (m^3/s^2). Default is GM_EARTH. Retained for API compatibility; the rate depends on G*J, not GM.
- Returns:
Lense-Thirring nodal precession rate (radians per second), positive (prograde).
- Return type:
Notes
For Earth satellites, classical J_2 nodal precession dominates by many orders of magnitude; Lense-Thirring is a small correction (~31 mas/year for LAGEOS), first measured with the LAGEOS satellites.
Examples
Lense-Thirring effect for the LAGEOS satellite:
>>> a = 12.27e6 # Semi-major axis (m) >>> e = 0.0045 >>> i = np.radians(109.9) >>> J_earth = 5.86e33 # Earth's spin angular momentum (kg m^2/s) >>> rate = lense_thirring_precession(a, e, i, J_earth) >>> mas_per_year = rate * 86400 * 365.25 * 206264.806 * 1e3 >>> print(f"LT nodal precession: {mas_per_year:.1f} mas/year") LT nodal precession: 30.7 mas/year
- pytcl.astronomical.relativity.relativistic_range_correction(r1, r2, rho, gm=398600441800000.0)[source]
Compute relativistic (Shapiro) range correction for ranging measurements.
When measuring the distance between two points through a gravitational field (e.g., satellite laser ranging or GNSS pseudoranging), spacetime curvature delays the signal. The equivalent one-way range correction is:
Δρ = (2*GM/c^2) * ln((r1 + r2 + ρ) / (r1 + r2 - ρ))
where r1 and r2 are the distances of the endpoints from the center of the gravitating body and ρ is the geometric distance between them.
- Parameters:
- Returns:
One-way range correction (m), always positive (the measured range is longer than the geometric range).
- Return type:
Examples
Shapiro range correction for GPS-to-ground ranging (satellite overhead):
>>> r_ground = 6.371e6 # Ground station radius (m) >>> r_gps = 26.561e6 # GPS orbit radius (m) >>> rho = r_gps - r_ground >>> correction = relativistic_range_correction(r_ground, r_gps, rho) >>> print(f"Range correction: {correction*1e3:.2f} mm") Range correction: 12.66 mm
Sgp4
SGP4/SDP4 Satellite Propagation Models.
This module implements the Simplified General Perturbations model (SGP4) and its deep-space extension (SDP4) for propagating satellite orbits from Two-Line Element (TLE) sets.
SGP4 models the effects of:
Atmospheric drag (via the B* term)
J2, J3, J4 gravitational harmonics
Secular and periodic variations
SDP4 additionally models (for orbital periods >= 225 min, computed from the recovered mean motion):
Lunar gravitational perturbations (secular and periodic)
Solar gravitational perturbations (secular and periodic)
Geopotential resonance for 12-hour and 24-hour orbits, integrated with the Euler-Maclaurin scheme of the reference algorithm
These are implemented as the four standard deep-space routines DSCOM, DSINIT, DSPACE and DPPER.
The output is in the TEME (True Equator, Mean Equinox) reference frame, which is a quasi-inertial frame used by NORAD.
Limitations
Only the WGS-72 gravity model is supported; the reference implementation also offers WGS-72 (old) and WGS-84.
Only the “improved” operation mode (
opsmode = 'i') is implemented, which is the recommended mode; the legacy AFSPC mode differs in the sidereal-time formula and in two angle-wrapping conventions.
References
Hoots, F. R. and Roehrich, R. L., “Spacetrack Report No. 3: Models for Propagation of NORAD Element Sets,” 1980.
Vallado, D. A., Crawford, P., Hujsak, R., and Kelso, T.S., “Revisiting Spacetrack Report #3,” AIAA 2006-6753.
Vallado, D. A., “Fundamentals of Astrodynamics and Applications,” 4th ed., Microcosm Press, 2013.
- class pytcl.astronomical.sgp4.SGP4State(r, v, error)[source]
Bases:
NamedTupleState vector from SGP4 propagation.
- Variables:
r (ndarray) – Position in TEME frame (km), shape (3,).
v (ndarray) – Velocity in TEME frame (km/s), shape (3,).
error (int) – Error code (0 = success).
- class pytcl.astronomical.sgp4.SGP4Satellite(tle)[source]
Bases:
objectSGP4 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.
- Variables:
Examples
>>> from pytcl.astronomical.tle import parse_tle >>> line1 = "1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9997" >>> line2 = "2 25544 51.6400 247.4627 0006703 130.5360 325.0288 15.49815350479003" >>> tle = parse_tle(line1, line2, name="ISS") >>> sat = SGP4Satellite(tle) >>> state = sat.propagate(0.0) # At epoch >>> bool(6700 < np.linalg.norm(state.r) < 6900) # ISS orbital radius (km) True >>> state = sat.propagate(60.0) # 60 minutes later
- propagate(tsince)[source]
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:
Examples
>>> from pytcl.astronomical.tle import parse_tle >>> tle = parse_tle( ... "1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9997", ... "2 25544 51.6400 247.4627 0006703 130.5360 325.0288 15.49815350479003") >>> 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
- pytcl.astronomical.sgp4.sgp4_propagate(tle, tsince)[source]
Propagate TLE using SGP4/SDP4 model.
Convenience function that creates an SGP4Satellite and propagates.
- Parameters:
- Returns:
state – Position and velocity in TEME frame.
- Return type:
Examples
>>> from pytcl.astronomical.tle import parse_tle >>> tle = parse_tle( ... "1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9997", ... "2 25544 51.6400 247.4627 0006703 130.5360 325.0288 15.49815350479003") >>> state = sgp4_propagate(tle, 60.0) # 60 minutes after epoch >>> bool(6700 < np.linalg.norm(state.r) < 6900) # ISS orbital radius (km) True
- pytcl.astronomical.sgp4.sgp4_propagate_batch(tle, times)[source]
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
>>> from pytcl.astronomical.tle import parse_tle >>> tle = parse_tle( ... "1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9997", ... "2 25544 51.6400 247.4627 0006703 130.5360 325.0288 15.49815350479003") >>> times = np.linspace(0, 90, 100) # 0 to 90 minutes >>> r, v = sgp4_propagate_batch(tle, times) >>> r.shape (100, 3)
- pytcl.astronomical.sgp4.unkozai_mean_motion(no, inclination, eccentricity)[source]
Recover the un-Kozai’d (Brouwer) mean motion from a TLE mean motion.
TLE mean motions are Kozai mean motions; SGP4 works with the Brouwer mean motion obtained by removing the secular J2 contribution. The recovered value also decides whether the deep-space model applies.
- Parameters:
- Returns:
Un-Kozai’d mean motion (rad/min).
- Return type:
Examples
>>> n = unkozai_mean_motion(0.0676, 0.9013, 0.0006703) >>> bool(n < 0.0676) True
Special Orbits
Special orbit cases: parabolic and advanced hyperbolic orbits.
This module extends orbital_mechanics.py with handling for edge cases: - Parabolic orbits (e = 1, unbounded trajectory with zero energy) - Advanced hyperbolic orbit calculations (escape trajectories) - Unified orbit type detection and handling
References
Vallado, D. A., “Fundamentals of Astrodynamics and Applications,” 4th ed., Microcosm Press, 2013.
Curtis, H. D., “Orbital Mechanics for Engineering Students,” 3rd ed., Butterworth-Heinemann, 2014.
Battin, R. H., “An Introduction to the Mathematics and Methods of Astrodynamics,” 2nd ed., AIAA, 1999.
- class pytcl.astronomical.special_orbits.OrbitType(value)[source]
Bases:
EnumClassification of orbit types based on eccentricity.
- CIRCULAR = 0
- ELLIPTICAL = 1
- PARABOLIC = 2
- HYPERBOLIC = 3
- class pytcl.astronomical.special_orbits.ParabolicElements(rp, i, raan, omega, nu)[source]
Bases:
NamedTupleParabolic (escape) orbit elements.
For a parabolic orbit (e=1), the semi-major axis is infinite. Instead, we use the periapsis distance and orientation parameters.
- Variables:
rp (float) – Periapsis distance (km). Also called pericenter or closest approach.
i (float) – Inclination (radians), 0 to pi.
raan (float) – Right ascension of ascending node (radians), 0 to 2*pi.
omega (float) – Argument of periapsis (radians), 0 to 2*pi.
nu (float) – True anomaly (radians), typically in [-pi, pi] for parabolic orbits.
- pytcl.astronomical.special_orbits.classify_orbit(e, tol=1e-09)[source]
Classify orbit type based on eccentricity.
- Parameters:
- Returns:
Classified orbit type.
- Return type:
- Raises:
ValueError – If eccentricity is negative or NaN.
- pytcl.astronomical.special_orbits.mean_to_parabolic_anomaly(M, tol=1e-12, max_iter=100)[source]
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:
- Returns:
D – Parabolic anomaly (the parameter D such that tan(nu/2) = D).
- Return type:
Notes
For parabolic orbits, mean anomaly relates to time as:
M = sqrt(mu / (2 * rp^3)) * t
with
rpthe periapsis distance andmuthe gravitational parameter. This used to be writtensqrt(mu/rp^3), missing the factor of 2 under the root – a sqrt(2) error in the time-to-anomaly conversion (gh-25). The solver below is unaffected: it takesMas given and is self-consistent either way, so only a caller converting from time was misled.The solution D satisfies: D + (1/3)*D^3 = M
- pytcl.astronomical.special_orbits.parabolic_anomaly_to_true_anomaly(D)[source]
Convert parabolic anomaly to true anomaly.
For parabolic orbits, the parabolic anomaly D relates to true anomaly by: tan(nu/2) = D
- pytcl.astronomical.special_orbits.true_anomaly_to_parabolic_anomaly(nu)[source]
Convert true anomaly to parabolic anomaly.
- pytcl.astronomical.special_orbits.mean_to_true_anomaly_parabolic(M, tol=1e-12)[source]
Direct conversion from mean to true anomaly for parabolic orbits.
- pytcl.astronomical.special_orbits.radius_parabolic(rp, nu)[source]
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:
- Returns:
r – Orbital radius (km).
- Return type:
- Raises:
ValueError – If radius would be negative (nu near +pi for parabolic orbit).
- pytcl.astronomical.special_orbits.velocity_parabolic(mu, rp, nu)[source]
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)
- pytcl.astronomical.special_orbits.hyperbolic_anomaly_to_true_anomaly(H, e)[source]
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:
- Returns:
nu – True anomaly (radians).
- Return type:
- Raises:
ValueError – If eccentricity is not hyperbolic (e <= 1).
- pytcl.astronomical.special_orbits.true_anomaly_to_hyperbolic_anomaly(nu, e)[source]
Convert true anomaly to hyperbolic anomaly.
- Parameters:
- Returns:
H – Hyperbolic anomaly (radians).
- Return type:
- Raises:
ValueError – If eccentricity is not hyperbolic.
- pytcl.astronomical.special_orbits.escape_velocity_at_radius(mu, r)[source]
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)
- pytcl.astronomical.special_orbits.hyperbolic_excess_velocity(mu, a)[source]
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:
- Returns:
v_inf – Hyperbolic excess velocity (km/s).
- Return type:
- Raises:
ValueError – If semi-major axis is not negative.
- pytcl.astronomical.special_orbits.hyperbolic_asymptote_angle(e)[source]
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:
- Raises:
ValueError – If eccentricity is not hyperbolic.
- pytcl.astronomical.special_orbits.hyperbolic_deflection_angle(e)[source]
Compute the deflection angle for a hyperbolic orbit.
The deflection (turn) angle is the angle through which the velocity vector is rotated between the incoming and outgoing asymptotes: delta = 2*arcsin(1/e), equivalently 2*nu_inf - pi.
- Parameters:
e (float) – Eccentricity (e > 1 for hyperbolic).
- Returns:
delta – Deflection angle (radians), in (0, pi).
- Return type:
- Raises:
ValueError – If eccentricity is not hyperbolic.
- pytcl.astronomical.special_orbits.semi_major_axis_from_energy(mu, specific_energy)[source]
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:
- 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:
- Raises:
ValueError – If
|specific_energy| < 1e-15(the near-parabolic case; not only exactly zero).
- pytcl.astronomical.special_orbits.eccentricity_vector(r, v, mu)[source]
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
Tle
Two-Line Element (TLE) set parsing and data structures.
This module provides functions for parsing NORAD Two-Line Element sets, which are the standard format for distributing satellite orbital data.
TLE Format
A TLE consists of two 69-character lines containing orbital elements in a specialized format used by NORAD/USSPACECOM.
- Line 1 format:
1 NNNNNC NNNNNAAA NNNNN.NNNNNNNN +.NNNNNNNN +NNNNN-N +NNNNN-N N NNNNN
Column Description 01 Line number (1) 03-07 Satellite catalog number 08 Classification (U=unclassified) 10-11 International designator (last two digits of launch year) 12-14 International designator (launch number of the year) 15-17 International designator (piece of the launch) 19-20 Epoch year (last two digits) 21-32 Epoch (day of the year and fractional portion of the day) 34-43 First derivative of mean motion (ballistic coefficient) 45-52 Second derivative of mean motion (decimal point assumed) 54-61 BSTAR drag term (decimal point assumed) 63 Ephemeris type 65-68 Element set number 69 Checksum (modulo 10)
- Line 2 format:
2 NNNNN NNN.NNNN NNN.NNNN NNNNNNN NNN.NNNN NNN.NNNN NN.NNNNNNNNNNNNNN
Column Description 01 Line number (2) 03-07 Satellite catalog number 09-16 Inclination (degrees) 18-25 Right ascension of ascending node (degrees) 27-33 Eccentricity (decimal point assumed) 35-42 Argument of perigee (degrees) 44-51 Mean anomaly (degrees) 53-63 Mean motion (revolutions per day) 64-68 Revolution number at epoch 69 Checksum (modulo 10)
References
Vallado, D. A., “Fundamentals of Astrodynamics and Applications,” 4th ed., Microcosm Press, 2013, Appendix C.
Hoots, F. R. and Roehrich, R. L., “Spacetrack Report No. 3: Models for Propagation of NORAD Element Sets,” 1980.
CelesTrak, https://celestrak.org/NORAD/documentation/tle-fmt.php
- class pytcl.astronomical.tle.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)[source]
Bases:
NamedTupleTwo-Line Element set data structure.
Contains all orbital elements and metadata from a NORAD TLE.
- Variables:
name (str) – Satellite name (from line 0, if present).
catalog_number (int) – NORAD catalog number.
classification (str) – Classification (‘U’ = unclassified).
int_designator (str) – International designator (e.g., ‘98067A’).
epoch_year (int) – Epoch year (4-digit).
epoch_day (float) – Epoch day of year (fractional).
ndot (float) – First derivative of mean motion (rev/day^2).
nddot (float) – Second derivative of mean motion (rev/day^3).
bstar (float) – BSTAR drag coefficient (1/Earth radii).
ephemeris_type (int) – Ephemeris type (usually 0 for SGP4).
element_set_number (int) – Element set number.
inclination (float) – Inclination (radians).
raan (float) – Right ascension of ascending node (radians).
eccentricity (float) – Eccentricity (dimensionless).
arg_perigee (float) – Argument of perigee (radians).
mean_anomaly (float) – Mean anomaly (radians).
mean_motion (float) – Mean motion (radians/minute).
revolution_number (int) – Revolution number at epoch.
line1 (str) – Original TLE line 1.
line2 (str) – Original TLE line 2.
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.
- pytcl.astronomical.tle.parse_tle(line1, line2, name='', verify_checksum=True)[source]
Parse a Two-Line Element set.
- Parameters:
- Returns:
tle – Parsed TLE data structure.
- Return type:
- Raises:
ValueError – If TLE format is invalid or checksum fails.
Examples
>>> line1 = "1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9997" >>> line2 = "2 25544 51.6400 247.4627 0006703 130.5360 325.0288 15.49815350479003" >>> tle = parse_tle(line1, line2, name="ISS (ZARYA)") >>> print(f"Inclination: {np.degrees(tle.inclination):.4f} deg") Inclination: 51.6400 deg
- pytcl.astronomical.tle.parse_tle_3line(lines)[source]
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:
Examples
>>> tle_text = '''ISS (ZARYA) ... 1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9997 ... 2 25544 51.6400 247.4627 0006703 130.5360 325.0288 15.49815350479003''' >>> tle = parse_tle_3line(tle_text) >>> print(tle.name) ISS (ZARYA)
- pytcl.astronomical.tle.tle_epoch_to_jd(tle)[source]
Convert TLE epoch to Julian date.
Examples
>>> line1 = "1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9997" >>> line2 = "2 25544 51.6400 247.4627 0006703 130.5360 325.0288 15.49815350479003" >>> tle = parse_tle(line1, line2) >>> jd = tle_epoch_to_jd(tle) >>> print(f"{jd:.1f}") # 2024 January 1.5 2460311.0
- pytcl.astronomical.tle.tle_epoch_to_datetime(tle)[source]
Convert TLE epoch to Python datetime.
- Parameters:
tle (TLE) – Parsed TLE.
- Returns:
dt – TLE epoch as UTC datetime.
- Return type:
datetime
- pytcl.astronomical.tle.format_tle(tle, include_name=True)[source]
Format TLE data structure back to TLE string.
- pytcl.astronomical.tle.is_deep_space(tle)[source]
Determine if TLE requires deep-space (SDP4) propagation.
Satellites with orbital period >= 225 minutes use SDP4 instead of SGP4 due to lunar-solar perturbations. The period is formed from the recovered (un-Kozai’d) mean motion, as in the reference algorithm – using the raw TLE mean motion picks the wrong model for element sets within ~0.1 min of the 225-minute boundary.