Kepler Orbit Propagation

Kepler Orbit Propagation (space-systems/orbit-mechanics/kepler-orbit-propagation)

Use when you must propagate a spacecraft orbit in time from its classical orbital elements: mean motion from Kepler's third law, the Kepler equation M = E - e sin E solved by Newton iteration, the eccentric-to-true anomaly conversion, the radius at any anomaly, the time since periapsis, and the full inertial position and velocity after an elapsed time. This leaf implements the standard two-body propagation model in pure Python, stdlib only. It pairs with space-systems/orbit-mechanics/keplerian-elements, which extracts elements from a state; this leaf is the forward direction, elements to state in time. It does NOT extract elements from a state vector, add perturbation drift, size impulsive maneuvers or solve two-position targeting; those are sibling leaves.

Domain quick reference

• Inputs: semimajor axis a (km), eccentricity e in [0, 1), inclination i, RAAN, argument of periapsis argp and initial true anomaly nu0 (rad), elapsed time dt (s), mu in km^3/s^2 (default 398600.4418, Earth).
• Mean motion: n = sqrt(mu / a^3); period T = 2 pi / n.
• Initial mean anomaly from nu0: E0 via the inverse half-angle map, then M0 = E0 - e sin E0. Propagation: M = M0 + n dt.
• Kepler equation: M = E - e sin E, solved for E by Newton iteration E -= (E - e sin E - M) / (1 - e cos E) from E = M + e sin M to a 1e-12 residual (cap 100 iterations).
• Eccentric anomaly to true anomaly: nu = 2 atan2(sqrt(1 + e) sin(E / 2), sqrt(1 - e) cos(E / 2)), branch-safe, folded to (-pi, pi]; the inverse map uses the mirrored half-angle form.
• Radius: r = a (1 - e cos E), identical to the conic form r = a (1 - e^2) / (1 + e cos nu).
• Inertial state: perifocal position r (cos nu, sin nu) and velocity sqrt(mu / p) (-sin nu, e + cos nu) with p = a (1 - e^2), rotated by M = R3(-RAAN) R1(-i) R3(-argp) (classical rotation matrix forms).
• Time since periapsis: t = (E - e sin E) / n for the anomaly nu, folded into [0, T).
• Non-physical inputs raise ValueError: a <= 0, mu <= 0, e outside [0, 1), dt < 0.

Workflow

1. Validate the element state and mu; fix a > 0, e in [0, 1), dt >= 0.
2. Get the mean motion and period with mean_motion and orbital_period.
3. Map the starting true anomaly nu0 back to E0 with eccentric_anomaly_from_true, then form M = M0 + n dt.
4. Solve the Kepler equation with kepler_solve(M, e) for E.
5. Convert E to nu with true_anomaly_from_eccentric and read the radius from r = a (1 - e cos E), cross-checked by radius_at_anomaly with the conic form.
6. Call propagate_kepler for the packed result: mean, eccentric and true anomaly, radius, position_km, velocity_kms and period_s.
7. For event timing, get the time since periapsis of any anomaly with time_since_periapsis; confirm the Kepler identities with the contract test scripts/test_kepler_orbit_propagation.py.

Worked example

Reference orbit a = 12000 km, e = 0.35, i = 30 deg, RAAN = 45 deg, argp = 20 deg, starting at periapsis (nu0 = 0), dt = 3600 s, Earth mu. Module outputs:

• n = 4.80283e-4 rad/s; T = 13082.262 s (3.634 h).
• M = 1.729018 rad; E = 2.041030 rad; nu = 2.336674 rad (133.8815 deg).
• r = 13902.9969 km; speed = 4.911570 km/s, consistent with the specific orbital energy relation v^2 = mu (2 / r - 1 / a).
• Inertial state: r = (-12575.051, -5079.004, 3060.248) km, v = (-0.292236, -4.579770, -1.750378) km/s.
• After one period (dt = T): nu returns to 0.000000 (mod 2 pi) and r = 7800.000000 km = a (1 - e); M = E = 2 pi.
• time_since_periapsis of the propagated nu recovers 3600.000000 s (inverse time of flight).
• Radius identity: a (1 - e cos E) equals a (1 - e^2) / (1 + e cos nu) to 1.8e-12 km. Recovering a and e from the propagated r, v state (energy and angular momentum identities) gives a to 3.0e-16 relative and e to 1.1e-16 absolute.

Pitfalls

• Starting the Newton solve at the wrong guess: the iteration launches from E = M + e sin M and caps at 100 iterations; feeding an initial anomaly far from the true branch still converges, but the high-eccentricity cases (e = 0.99 in the contract test) are where a poor guess costs sweep count and accuracy.
• Confusing mean and eccentric anomaly: M is the time-linear quantity (M = M0 + n dt) and E is the transcendental one solved from the Kepler equation; mixing them in the radius formula r = a (1 - e cos E) corrupts the state.
• Using the atan form of the anomaly map: the leaf converts with the branch-safe nu = 2 atan2(sqrt(1+e) sin(E/2), sqrt(1-e) cos(E/2)) folded to (-pi, pi]; the naive arctan of the tangent loses the quadrant and the round trip breaks.
• Treating dt as elapsed from nu0 without the mean-anomaly wrap: propagation adds n dt to M0, so a dt of many periods folds correctly only because M stays the running quantity; the final anomaly is then folded, not the time.
• Forgetting that a = 12000 km is a radius, not an altitude: every element and output uses the orbital radius scale (the worked example returns r = 13902.997 km), so altitude inputs must be converted before the call.
• Assuming extraction works here too: this leaf is the forward map, elements to state in time; recovering elements from a state vector belongs to keplerian-elements, so do not look for an inverse element solver in this module.

Verification

• Confirm kepler_solve(1.729018..., 0.35) returns E = 2.041030 rad with a residual below 1e-12 and that E = M exactly for e = 0.
• Confirm the radius from r = a (1 - e cos E) matches the conic form to 1e-9 relative at the worked state.
• Confirm a full-period propagation returns the initial anomaly and radius a (1 - e) to 1e-9, and that dt = 0 leaves the state unchanged.
• Confirm time_since_periapsis inverts the propagation for any dt inside one period (to 1e-6 and better).
• Confirm the propagated r, v lie in the orbital plane: the angular momentum direction equals (sin RAAN sin i, -cos RAAN sin i, cos i) to 1e-9 and |h| = sqrt(mu a (1 - e^2)).
• Confirm ValueError rejection of a <= 0, mu <= 0, e < 0 or e >= 1, dt < 0 and malformed vectors.
• Run the contract test offline: python3 scripts/test_kepler_orbit_propagation.py (35 tests, deterministic).

Related leaves

• space-systems/orbit-mechanics/keplerian-elements: the inverse map, a position and velocity state back to classical elements; this leaf is elements to state in time.
• gnc-autonomy/space/orbit-dynamics: the two-body mission context that consumes propagated states.
• space-systems/orbit-mechanics/lambert-transfer: two-position targeting, the boundary case this leaf does not cover.

Behavior contract (gate 3)

Run the deterministic contract test (stdlib unittest, offline):

python3 scripts/test_kepler_orbit_propagation.py

The test covers the worked-example anchors (mean motion, period, M, E, nu, radius, speed), the Newton solver residual across anomaly and eccentricity grids including e = 0.99, exact e = 0 behavior, the anomaly maps and their round trips, the radius identity, the one-period return, dt = 0 invariance, the inverse time of flight, the perifocal rotation against the closed-form periapsis direction, the angular momentum direction and magnitude, orbit-plane membership, determinism, and ValueError rejection of every non-physical input class.

Compliance

• Standards referenced, not reproduced: ECSS series text is copyright ESA and freely downloadable (standards-map.yaml); the two-body propagation relations above are common astrodynamics methodology, summary-only.
• compliance: STANDARDS-REF, gated: false.