Native astrodynamics in Rust

pykep-rust is an independent native Rust port of the numerical C++ library in pykep version 3 (kep3). The reusable pykep-core crate contains the astrodynamics implementation without a C or C++ runtime dependency. The optional pykep-rust Python distribution exposes the same implementation through PyO3.

This book connects the narrative material for both interfaces:

• examples and quick starts provide runnable Rust and Python entry points;
• numerical conventions define units, epochs, frames, array layouts, tolerances, and error behavior;
• the dynamics, ephemeris, propagation, and low-thrust guides explain algorithm-specific contracts;
• validation and stabilization evidence record parity, independent checks, coverage, performance, and release limitations.

Exact Rust types, methods, and compiled examples are in the generated pykep-core API reference. Python users should start with the Python API contract and migration matrix.

This project is not an official ESA release. Its pinned upstream source and MPL-2.0 adaptation policy are recorded in the GitHub repository.

Documentation

The documentation describes the implemented native Rust core and its thin Python interface:

• AI problem-solving context: selecting physical models, APIs, units, frames, solver settings, and validation for a user problem;
• examples.md: Rust/Python quick starts and runnable examples;
• conventions.md: units, epochs, frames, array layouts, and numerical behavior;
• dynamics.md, pontryagin.md, zero-order-hold.md, ephemerides.md, low-thrust-legs.md, and zoh-leg.md: module guides;
• python-api.md: Python units, arrays, defaults, errors, ownership, GIL, and typing contract;
• python-migration.md: upstream C++/Python names, deliberate differences, deferrals, and unsupported ecosystem modules;
• decisions/: architecture and dependency decisions;
• performance.md: benchmark methodology and results;
• stabilization.md: matched distributions, profiling, regression limits, Miri/fuzz/Valgrind evidence, and release blockers;
• status.md, source-map.md, and validation.md: limitations, provenance, and validation evidence;
• development.md: contributor commands and repository policy.

Documentation is updated only for implemented behavior. Explicitly deferred or unavailable functionality is labelled as such in the migration and status documents.

Examples and quick starts

The examples are deliberately small and deterministic. Runtime statements are orientation for a release build, not performance guarantees; use the maintained benchmark harnesses for measurements.

Rust quick start

Add the native crate and propagate a normalized circular orbit:

use pykep_core::astro::propagation::propagate_lagrangian;

fn main() -> pykep_core::Result<()> {
    let initial = [1.0, 0.0, 0.0, 0.0, 1.0, 0.0];
    let final_state =
        propagate_lagrangian(&initial, core::f64::consts::FRAC_PI_2, 1.0)?;
    assert!((final_state[1] - 1.0).abs() < 1e-12);
    Ok(())
}

From this workspace, run cargo run --release -p pykep-examples --bin propagation. The default feature set includes the embedded VSOP2013 coefficients; use pykep-core with default-features = false when they are not needed.

Python quick start

Install a wheel in a virtual environment and use the same Rust core:

python -m venv .venv
.venv/bin/python -m pip install target/wheels/pykep_rust-*.whl
.venv/bin/python python/examples/elements_propagation.py

import numpy as np
import pykep_rust as pk

states = np.tile([1.0, 0.0, 0.0, 0.0, 1.0, 0.0], (4096, 1))
times = np.linspace(0.0, 1.0, len(states))
result = pk.propagate_lagrangian_batch(states, times, 1.0, workers=8)
assert result.shape == (4096, 6)

The Python package ships type information and accepts NumPy batch arrays only at float64; see python-api.md and batch-processing.md.

Runnable matrix

Every Rust binary is compiled by the workspace test and clippy gates. Every Python script is executed by python/tests/test_examples.py, including from the clean-wheel CI matrix.

Documentation map

• conventions.md defines units, epochs, frames, layout, and accuracy interpretation.
• Dynamics, Pontryagin, zero-order-hold, ephemeris, Sims–Flanagan, and generic ZOH-leg behavior have dedicated module guides in this directory.
• python-migration.md maps the C++/Python upstream surface and records deliberate differences and unavailable ecosystem areas.
• The decisions/ records explain fixed-size types, errors, VSOP data, and adaptive-integration dependencies.
• performance.md gives benchmark methodology and results.
• status.md, source-map.md, and validation.md state implemented scope, limitations, and evidence without implying support for deferred modules.

Numerical conventions

Unless an API states otherwise:

• scalar computations use IEEE-754 binary64 (f64);
• position is in metres and velocity is in metres per second;
• gravitational parameters are in cubic metres per square second;
• durations are in seconds, while Julian-date conversions operate in days;
• angles are in radians;
• scalar ephemeris epochs use MJD2000;
• Cartesian state ordering is [x, y, z, vx, vy, vz];
• classical element ordering is [a, e, i, Ω, ω, ν];
• modified equinoctial ordering is [p, f, g, h, k, L];
• matrices use row-major nested arrays in Rust and C-contiguous row-major arrays at the Python boundary;
• public functions reject NaN and positive or negative infinity;
• deterministic batch APIs preserve input order; workers=0 uses Rayon’s shared pool, workers=1 is serial, and workers=N uses a cached pool with exactly N workers.

Some upstream functions propagate non-finite values or use them as invalid domain sentinels. The Rust API reports explicit errors instead. Algorithm guides document any narrower valid domain and all intentional deviations.

Classical conversion reports circular and equatorial states as singular because their node/periapsis angles are undefined. Modified equinoctial elements cover those states: choose the prograde convention except at inclination π, and the retrograde convention except at inclination zero. Element Jacobians are 6 × 6, with output components as rows and input components as columns.

Two-body propagation accepts any caller-consistent position, velocity, time, and gravitational-parameter units. Negative durations propagate backward. Time grids are relative to their first entry. State-transition matrices use ∂state_final/∂state_initial, with output components as rows.

Lambert solutions are ordered deterministically: the zero-revolution solution first, then left and right solutions for each increasing revolution count. clockwise = false selects prograde motion as viewed from positive z; collinear endpoints are rejected because that automatic direction convention is undefined.

Ephemeris scalar epochs use MJD2000 days. Returned states use SI units for built-in Solar System providers and caller-consistent units for constructed Keplerian providers. Provider metadata uses Option rather than upstream negative sentinels, and hyperbolic periods are None.

Ephemerides

JPL low-precision planets

JplLowPrecision and Planet.jpl_low_precision() evaluate the eight heliocentric approximate planetary models carried by pykep/kep3 3.0.1: Mercury, Venus, Earth, Mars, Jupiter, Saturn, Uranus, and Neptune. Body lookup is ASCII case-insensitive. The returned name is the lowercase body followed by (jpl_lp).

Inputs are MJD2000 day counts and must satisfy the open interval -73048 < epoch < 18263, approximately 1800–2050. Cartesian output is [x, y, z, vx, vy, vz] in metres and metres per second, heliocentric and referred to the mean ecliptic and equinox of J2000. The underlying JPL table uses Julian ephemeris date/JDTDB; this library does not perform time-scale conversion. The earth coefficients describe the Earth–Moon barycentre, as in the source table.

The coefficients and rates come from the JPL Solar System Dynamics approximate-position table. JPL describes these as lower-accuracy fitted formulae and warns against using them outside their fitted interval. Its nominal 1800–2050 errors are:

These models are suitable for approximate mission design, not precision navigation. Use a high-precision integrated ephemeris when those error bounds are too large.

The Rust provider exposes true-anomaly, mean-anomaly, and both modified equinoctial element forms. Python scalar and NumPy state, element, period, and optional-acceleration batches call the same provider. Batch order is preserved; workers=0 uses the shared pool, one is serial, and larger values select an exact cached worker count.

VSOP2013

Vsop2013 and Planet.vsop2013() implement the analytical theory for Mercury, Venus, the Earth–Moon barycentre (earth_moon), Mars, Jupiter, Saturn, Uranus, Neptune, and Pluto. The coefficient source is the IMCCE VSOP2013 solution; the exact adaptation and license chain are recorded beside the embedded data and in ADR 0003.

Input is an MJD2000 day count interpreted as TDB. The evaluator applies T = (mjd2000 - 0.5) / 365250, because the theory is measured in thousands of Julian years from J2000 at JD 2451545.0, twelve hours after the MJD2000 origin. Output is heliocentric ICRF [x, y, z, vx, vy, vz] in metres and metres per second. The provider does not convert UTC, TAI, TT, or TDB.

The theory was fitted to INPOP10a over 1890–2000. IMCCE also publishes comparison errors over −4000 to +8000, but accuracy degrades with distance from the fit interval and especially for Pluto. There is no artificial hard date cutoff; callers must choose a time span appropriate to their accuracy requirements. Over the fit interval, the published largest heliocentric longitude/latitude/distance differences are:

The vsop2013 Cargo feature is enabled by default. It embeds 4.3 MiB of coefficients and supports thresholds greater than or equal to 1e-9; the upstream default is 1e-5. Smaller thresholds are rejected because the remaining 2.5 million terms are intentionally not embedded. Disable default features to omit the data and retain the Keplerian and JPL low-precision providers. Vsop2013::available() and the corresponding Python static method make the build configuration queryable.

Evaluated dynamics

Phase 11 provides three stateless six-state models:

• KeplerDynamics for inertial two-body Cartesian motion;
• Cr3bpDynamics for the circular restricted three-body problem;
• BcpDynamics for the time-dependent bicircular problem.

All states are ordered [x, y, z, vx, vy, vz]. Kepler inputs use one consistent dimensional unit system: if position and time are SI, mu is in m³/s². CR3BP and BCP use the upstream nondimensional rotating-frame convention. Their primaries are fixed at (-mu, 0, 0) and (1 - mu, 0, 0). The BCP Sun is at rho_sun [cos(omega_sun t), sin(omega_sun t), 0].

mu is the secondary-to-total mass ratio and accepts the mathematical range [0, 1]; the customary ordering has mu <= 0.5. BCP parameters are [mu, mu_sun, rho_sun, omega_sun]. mu_sun must be non-negative and rho_sun positive. The constants module supplies the upstream Earth–Moon–Sun defaults.

Evaluation and propagation

Each model exposes evaluate, propagate, and propagate_with_stm. The propagators use the Phase 10 DOP853 facade. State-transition matrices are row-major, with output state components in rows and initial state components in columns. Analytic state and parameter Jacobians implement the generic sensitivity contract.

For comparisons with the pinned C++ Taylor-adaptive implementation, the test profile uses relative and absolute tolerances of 2e-13 and a maximum step of 0.01 nondimensional time units. Across the committed sample grid this supports state tolerances of 3e-11 for Kepler/BCP and 2e-9 for the longer, close-approach CR3BP trajectory. These are validation tolerances, not a guarantee for arbitrary trajectories; callers must choose tolerances based on their scale, duration, and invariant drift requirements.

CR3BP also exposes the positive effective potential

U = (x² + y²)/2 + (1-mu)/r1 + mu/r2

and the Jacobi constant C = 2 U - |v|².

Python boundary

Python exposes kepler_rhs, cr3bp_rhs, bcp_rhs, the CR3BP potential and Jacobi functions, and model-specific propagation functions. The propagation calls release the GIL and accept initial_time, scalar tolerances, and an optional maximum step. Variational variants return (state, stm) without exposing heyoka or any internal expression graph.

Exact collisions with the active model bodies are SingularGeometryError. Invalid mass or distance parameters are ValueError. A failure to finish an otherwise valid propagation is IntegrationError; tests keep those model and integrator failure categories separate.

Unlike upstream, no compiled-expression cache is needed: the evaluated model types are zero-sized and integration configuration is local to each call.

Zero-order-hold dynamics

Phase 12 provides a common ControlSchedule<C> and four evaluated models:

• ZohKeplerDynamics: [x,y,z,vx,vy,vz,mass], normalized mu = 1;
• ZohCr3bpDynamics: the same seven-state layout in the CR3BP synodic frame;
• ZohEquinoctialDynamics: [p,f,g,h,k,L,mass], normalized mu = 1;
• ZohSolarSailDynamics: six Cartesian states with cone/clock attitude.

The first three control rows are [thrust, i1, i2, i3]. The direction is Cartesian for Kepler/CR3BP and radial-transverse-normal for equinoctial dynamics. Kepler/equinoctial constants contain mass-flow coefficient c; CR3BP constants are [c, mu]. The preserved upstream mass equation is dm/dt = -c thrust exp(-1 / (mass 1e16)).

Solar-sail rows are [alpha, beta] in radians. Its constant c scales ideal sail acceleration as c cos(alpha)^2 / r².

Schedule and switch contract

A schedule has S + 1 strictly increasing finite boundaries and exactly S finite control rows. Construction performs all dimension, finiteness, and monotonicity checks.

For lookup, segment i owns [t_i, t_(i+1)); the final boundary belongs to the final segment. Forward propagation ends each integration exactly at a switch and starts the next segment with its new control. Backward propagation uses the interval to the left of each encountered switch. This directional rule avoids evaluating a segment across a discontinuity.

Each segment is integrated exactly once. The active control is copied into a fixed-size parameter array before the solve, so no allocation or schedule scan occurs in an RHS call. control_at uses binary search for occasional external lookup.

Sensitivities

ZohSensitivitySeeds carries an arbitrary compile-time seed width through the schedule in one augmented solve per segment. It contains:

• initial state seeds;
• one control seed matrix per segment;
• constant-parameter seeds.

State sensitivities are continuous at switches. A control column for a future segment remains exactly zero until that segment becomes active. Runtime is linear in the segment count for a fixed seed width; the implementation does not repropagate every prior segment for every control.

Model Jacobians use fixed-size, allocation-free central differentiation of the evaluated source equations. They are checked against C++ Taylor variations and end-to-end schedule finite differences. The longest CR3BP reference uses a 2e-5 relative/absolute variation comparison; nominal state parity remains 3e-10. This looser derivative tolerance is explicit and will be revisited if later leg gradients require analytic Jacobians.

Python API

Python exposes all four RHS functions and propagate_zoh_* functions. Boundaries and controls are ordinary nested sequences. Propagation releases the GIL and accepts backward, scalar tolerances, and maximum_step. Malformed row widths and grids are rejected before native integration.

Pontryagin low-thrust dynamics

Phase 13 provides evaluated Cartesian and modified-equinoctial canonical dynamics for indirect low-thrust optimization. The implementation is native Rust and uses the common DOP853 integration and sensitivity contracts; it does not expose or depend on a symbolic-expression runtime.

Orders and units

The Cartesian augmented state is

[x, y, z, vx, vy, vz, m, lx, ly, lz, lvx, lvy, lvz, lm]

The modified-equinoctial augmented state is

[p, f, g, h, k, L, m, lp, lf, lg, lh, lk, lL, lm]

The equinoctial convention is prograde and thrust direction is expressed in the radial, transverse, normal (RTN) frame. Cartesian thrust direction is inertial. p, position, velocity, mass, time, mu, thrust, and exhaust velocity must form one consistent unit system. Angles are radians. Costate units follow from the caller’s chosen normalization.

Mass-optimal models use the parameter order

[mu, maximum_thrust, exhaust_velocity, barrier, lambda0]

barrier is the positive logarithmic regularization coefficient. Time-optimal models use

[mu, maximum_thrust, exhaust_velocity]

and have full throttle with the upstream implicit lambda0 = 1. Maximum thrust may be zero, which is useful for canonical coordinate checks; all other physical scale parameters and lambda0 must be positive.

Control and Hamiltonian

The public control evaluators return OptimalControl, containing throttle, direction, and the switching function rho. Mass optimality uses the regularized upstream control

u = 2 eps / (rho + 2 eps + sqrt(rho² + 4 eps²))

evaluated with an algebraically equivalent branch that avoids cancellation. Time optimality returns u = 1. The Hamiltonian functions evaluate the minimized autonomous Hamiltonian using that control.

The minimizing direction is undefined when the relevant primer norm is exactly zero. Rust reports PykepError::SingularGeometry, and Python reports SingularGeometryError; no arbitrary direction or NaN vector is returned. Cartesian position at the origin, non-positive mass or equinoctial semilatus rectum, and a zero equinoctial radial denominator are also explicit errors.

Models and sensitivities

The four zero-sized model types are:

• CartesianMassOptimal
• CartesianTimeOptimal
• EquinoctialMassOptimal
• EquinoctialTimeOptimal

Each implements DynamicsModel<14, P> and DifferentiableDynamicsModel<14, P>. Canonical costate rates use forward-mode derivatives of the minimized Hamiltonian. Full state and parameter Jacobians use centered fixed-size differences and are stored as jacobian[output][input]. This path performs no heap allocation in the model right-hand side.

Python requires the native Optimality.Mass or Optimality.Time enum. Unvalidated strings are not accepted. The pontryagin_*_rhs, control, Hamiltonian, and propagation functions use the same pykep-core implementation.

Validation

The committed C++ oracle covers both coordinates and both objectives, their eight upstream variational arguments (seven initial costates and lambda0 in mass mode), and the upstream dimensional 100-day trajectory. Independent checks cover Hamiltonian conservation, propagated central differences, Jacobian orientation, zero-primer errors, normalized controls, and canonical Hamiltonian agreement after transforming costates between coordinates with the analytic element Jacobian.

At 2e-12 DOP853 scalar tolerances and a 0.01 maximum step for normalized cases, nominal trajectories agree with the pinned C++/heyoka oracle within a scaled 3e-10. The dimensional case uses a 12-hour maximum step and a scaled 3e-9. First-order variations use a scaled 2e-4 tolerance because the generic Jacobian contract is numerical.

Sims–Flanagan low-thrust legs

pykep-core::leg provides the fixed-duration SimsFlanaganLeg and the variable-duration SimsFlanaganAlphaLeg. Both are immutable after validated construction, use native Rust two-body propagation, and have no C or C++ runtime dependency.

Units and ordering

All inputs must use one self-consistent unit system:

• endpoint state is [x, y, z, vx, vy, vz];
• endpoint mass and maximum thrust use compatible mass and force units;
• time of flight, segment durations, and exhaust velocity use compatible time units;
• mu has length cubed per time squared;
• throttle vectors are dimensionless and appear in chronological order.

An endpoint is represented by SpacecraftEndpoint, which rejects non-finite state components, zero radius, and non-positive mass. SimsFlanaganSettings rejects negative time of flight or maximum thrust, non-positive exhaust velocity or mu, non-finite values, and cuts outside [0, 1]. At least one three-component throttle is required.

Transcription and cut

For a leg with N segments, the number propagated forward from departure is

floor(N * cut)

and the remaining segments are propagated backward from arrival. A cut of zero therefore starts at the departure endpoint and propagates every segment backward; a cut of one propagates every segment forward.

Each segment applies its finite impulse at the segment midpoint. The fixed-duration leg uses duration time_of_flight / N for every segment. The velocity increment has magnitude maximum_thrust * duration * |throttle| / mass; the mass update follows the Tsiolkovsky exponential using exhaust_velocity. Backward propagation reverses both the impulse and mass update.

The alpha leg accepts direct non-negative segment durations. As in the upstream class, their sum is not required to equal time_of_flight. SimsFlanaganAlphaLeg::from_time_weights is the explicit alternative that requires positive weights and normalizes them to the configured time of flight.

Constraints and derivatives

mismatch_constraints() returns seven values in this order:

[forward_rx - backward_rx,
 forward_ry - backward_ry,
 forward_rz - backward_rz,
 forward_vx - backward_vx,
 forward_vy - backward_vy,
 forward_vz - backward_vz,
 forward_mass - backward_mass]

throttle_constraints() returns N values, dot(throttle_i, throttle_i) - 1. A zero value is the unit-throttle limit; negative values are inside it.

The fixed leg exposes an analytic mismatch_jacobian() with output-by-input rows:

throttle_jacobian() has shape N × 3N in the same flattened control order. At zero throttle the mass-direction derivative uses the defined zero subgradient, so the returned matrix remains finite.

The upstream alpha class does not expose an analytic mismatch gradient, and neither does this port. Its exact throttle Jacobian remains available.

Python

pykep_rust.SimsFlanaganLeg and pykep_rust.SimsFlanaganAlphaLeg expose the same immutable evaluations. Python states and controls are converted once at construction; all astrodynamics calculations call pykep-core. The fixed-leg mismatch_jacobian() returns a tuple of the three row matrices described above.

Construction and evaluation raise ValueError for invalid values or shapes and RuntimeError for propagation failures.

Generic zero-order-hold leg

pykep_core::leg::ZohLeg transcribes a transfer between two fixed endpoint states as non-uniform segments with continuous, piecewise-constant controls. It is generic over any Rust model implementing ZeroOrderHoldModel; aliases cover every built-in ZOH dynamics family:

The units and frame conventions of each row are those of its model in zero-order-hold.md. c is the normalized mass-flow coefficient for the three low-thrust models and the normalized lightness coefficient for the solar sail.

Grid, controls, and cut

For S segments, construction requires exactly S + 1 finite, strictly-increasing time-grid nodes and S finite control vectors. Controls are stored in chronological order and own the half-open interval [time_grid[i], time_grid[i+1]); the final endpoint belongs to the final segment.

The cut uses the same convention as the upstream leg:

forward_segments = floor(S * cut)
backward_segments = S - forward_segments

The forward state starts at the initial endpoint and the backward state starts at the final endpoint. The mismatch is their component-wise difference at the cut. Cuts zero and one are supported without special placeholder states.

The leg is immutable after construction. It rejects invalid endpoint states, model constants, dimensions, grids, cuts, tolerances, maximum steps, and maximum-step magnitudes before evaluation.

Sensitivity layout

mismatch_jacobian() returns four output-by-input matrices:

Each segment propagates a local state-transition and active-control matrix. The leg composes those matrices from each side of the cut and includes the dynamics jump at every switching time. Constant model parameters are fixed leg configuration and therefore are not a returned derivative group.

The four built-in models compute their local RHS Jacobians with fixed-size centered differences using a relative step of 3e-6. Consequently, integrator tolerances such as 1e-12 do not imply 1e-12 derivative accuracy: the pinned end-to-end validation uses scaled tolerances up to 3e-5. See zero-order-hold.md for the model-level contract and validation evidence.

Integration and failures

IntegratorOptions applies independently to every segment, including relative and absolute tolerances, optional initial and maximum step sizes, maximum steps, and maximum rejections. A propagation failure is reported as an IntegrationFailure containing the direction, chronological segment index, and attempted time interval.

state_history(samples_per_segment) uses one DOP853 dense solve per segment for uniformly spaced states including both endpoints. The selected backend’s decreasing-time dense interpolation can panic on a rounded step boundary, so backward segments are evaluated through the equivalent increasing coordinate tau = segment_start - time rather than handed to that backend path. At least two samples are required. Backward histories list the final segment first, matching propagation order. evaluate_zoh_mismatch_batch() evaluates validated Rust legs in input order.

Python

pykep_rust.ZohLeg accepts a ZohModel value: Kepler, Cr3bp, Equinoctial, or SolarSail. The constructor validates dynamic state, control, and constant dimensions for that selection. Mismatch, Jacobian, and history evaluation release the Python GIL.

The ZOH-leg batch methods clone the small immutable leg descriptors, release the GIL once, and return results in input order. workers=0 uses the shared pool, workers=1 is serial, and a larger value selects a cached pool of that exact size.

Python API contract

Install the native wheel and import the collision-safe package:

import numpy as np
import pykep_rust as pk

epoch = pk.Epoch.from_iso("2030-01")
state = pk.classical_to_cartesian(
    [7.0e6, 0.01, 0.4, 1.0, 0.5, 0.2], pk.MU_EARTH
)
future = pk.propagate_lagrangian(state, 60.0, pk.MU_EARTH)

batch = np.asarray([state, future], dtype=np.float64)
times = np.asarray([60.0, -60.0], dtype=np.float64)
propagated = pk.propagate_lagrangian_batch(batch, times, pk.MU_EARTH)
assert propagated.shape == (2, 6)

parallel = pk.propagate_lagrangian_batch(
    batch, times, pk.MU_EARTH, workers=4
)
assert np.array_equal(parallel, propagated)

Units, shapes, and defaults

• Angles are radians. Epoch and ephemeris arguments are MJD2000 days. Propagation durations are seconds when the supplied gravitational parameter uses SI units; normalized inputs produce normalized output.
• Cartesian states are [x, y, z, vx, vy, vz]. Classical elements are [a, e, i, Ω, ω, ν]. Modified equinoctial elements are [p, f, g, h, k, L].
• Scalar vector inputs accept finite Python sequences. Element and propagation batches require two-dimensional float64 arrays with shape N × 6; one-dimensional epoch/time batches require float64 arrays. Strided and read-only arrays are accepted and outputs are newly owned.
• Every parallel batch accepts workers=0 for the shared pool, workers=1 for a serial native loop, or an exact positive worker count. Results and the first reported error retain input order.
• Jacobians are row-major, output-by-input. State-transition matrices are 6 × 6. Sims–Flanagan and ZOH leg matrix shapes are documented with their transcription in the low-thrust guides.
• Adaptive propagators default to relative and absolute tolerance 1e-12, no maximum step, and an implementation limit of 100,000 accepted/rejected steps where exposed. ZOH legs default to cut=0.5.

Errors and ownership

NaN, infinity, invalid dimensions, non-positive physical parameters, and invalid time grids are rejected before numerical work. Invalid input uses ValueError (or TypeError when a NumPy dtype/rank cannot satisfy the typed buffer). Numerical failures derive from PykepError: ConvergenceError, SingularGeometryError, IntegrationError, and UnsupportedCapabilityError.

Objects copy constructor inputs and can outlive those Python sequences. Planet and immutable leg objects can be reused from multiple Python threads. Batch conversion, anomaly, Stumpff, ephemeris, mission, propagation, dynamics, and leg workloads release the GIL while executing the native loop. Python wrappers only validate and convert data; all astrodynamics formulas live in pykep-core.

See batch-processing.md for the complete batch matrix, array shapes, worker-pool contract, SpOC 4 motivation, and guidance for avoiding nested parallelism.

The complete callable surface and defaults are statically described by the shipped py.typed marker and _pykep_rust.pyi. See python-migration.md for the upstream name and behavior mapping, conventions.md for numerical conventions, and the model-specific dynamics and leg guides for parameter layouts.

Ordered parallel batch processing

pykep-core and pykep_rust deliberately extend the scalar-oriented pykep surface with ordered parallel batches. This is not an upstream compatibility claim. Experience optimizing large Lambert-transfer populations during the SpOC 4 competition, captured in the dietmarwo/pykep-lambert comparison project, motivated moving the reusable pattern into this library.

The extension addresses two independent costs:

1. one Python-to-native transition replaces thousands of scalar calls; and
2. independent native evaluations can use several CPU cores.

Every batch returns results in input order. A successful row is numerically the same scalar computation; batching does not change the physical model, units, frame, tolerances, Lambert branch order, or error criteria.

Worker contract

APIs with a workers argument use one shared contract:

Explicit pools are cached by worker count, so repeated optimizer generations do not pay thread-startup cost. Explicit counts above 1024 are rejected. Python releases the GIL around native work.

Parallel evaluation collects row results before returning an error. This keeps both successful output and the reported first error deterministic in input order, even if a later failing row finishes first on another worker. An empty batch returns an empty result with the documented output shape.

Do not add another thread pool around a parallel batch. In an already parallel optimizer or task scheduler, pass workers=1 at the inner level to avoid oversubscription. Cheap arithmetic often benefits from one native batch because it removes Python-call overhead but may not benefit from multiple threads. Measure representative batch sizes.

Rust interfaces

The following astrodynamics operations have named Rust batch APIs:

• propagate_lagrangian_batch, propagate_universal_batch, and propagate_keplerian_batch;
• propagate_lagrangian_with_stm_batch and propagate_lagrangian_grid_parallel;
• solve_lambert_batch, using one LambertRequest per problem;
• Ephemeris::states_parallel, accelerations_parallel, periods_parallel, and elements_parallel;
• all anomaly-conversion functions with a _batch suffix;
• dot_batch, norm_batch, normalize_batch, cross_batch, and skew_batch;
• evaluate_zoh_mismatch_batch_parallel.

Other immutable scalar Rust operations can use the same public pykep_core::batch::try_map executor. This avoids duplicating dozens of near-identical request structures while preserving typed scalar functions as the only numerical implementation:

#![allow(unused)]
fn main() {
use pykep_core::astro::transfers::hohmann;
use pykep_core::batch::try_map;

let radii = [(1.0, 2.0), (1.2, 2.4)];
let transfers = try_map(&radii, 2, |&(r1, r2)| hohmann(r1, r2, 1.0))?;
assert_eq!(transfers.len(), radii.len());
Ok::<(), pykep_core::PykepError>(())
}

try_map is appropriate only for independent rows. Algorithms that carry state from one time to the next remain scalar/sequential unless their public contract explicitly defines independent initial-value problems.

Python interfaces

Fixed-size numeric batches use float64 NumPy arrays. Ragged schedules, encoding vectors, and batches of immutable leg objects use Python sequences. The shipped _pykep_rust.pyi is the authoritative signature reference.

Foundations and representations

Vector input has shape N × 3; element/state input has shape N × 6. Vector outputs have shape N or N × 3, skew matrices N × 3 × 3, element outputs N × 6, and element Jacobians N × 6 × 6.

Propagation, Lambert, and ephemerides

Propagation batches pair row states[i] with times[i] or final_times[i]. STM batches return (states[N,6], stms[N,6,6]). lambert_problem_batch pairs initial_positions[N,3], final_positions[N,3], and times[N], while sharing mu, direction, and maximum-revolution configuration. It returns N normal LambertProblem objects; each object retains the scalar zero/left/right branch ordering.

states = np.tile(
    np.asarray([1.0, 0.0, 0.0, 0.0, 1.0, 0.0]), (32_768, 1)
)
times = np.linspace(0.0, 1.0, len(states))
propagated = pk.propagate_lagrangian_batch(
    states, times, 1.0, workers=8
)
assert propagated.shape == states.shape

Mission utilities

Hohmann and bi-elliptic batches return (delta_v[N], time[N], impulses[N,K]). Flyby constraint output is N × 2, its Jacobian is N × 2 × 6, and outgoing velocity is N × 3. MIMA variants return (mass[N], acceleration[N]).

Controlled dynamics and legs

ZOH schedule batches accept one state row, boundary vector, and control matrix per independent schedule. Shared physical constants and integrator settings remain scalar arguments. Pontryagin batches share one optimality mode and parameter vector.

Where the primitives are useful

These batches support existing computations without adding new physical models:

• Lambert departure/arrival/time-of-flight sweeps;
• pork-chop or transfer-table calculations assembled from ephemeris and Lambert batches;
• optimizer populations for CMA-ES, differential evolution, PGPE, MODE, or similar methods;
• Monte Carlo evaluation of independent uncertain initial conditions;
• state/STM ensembles for covariance or sensitivity workflows;
• independent flyby candidates, MIMA estimates, and low-thrust legs.

Closest-approach, eclipse, access, conjunction, and other event physics are not implemented merely because the underlying state batches make such external calculations easier. Use an appropriate validated model for those quantities.

Validation

Rust tests compare named batches with scalar functions for serial, global, and explicit worker pools and verify deterministic error order. Python tests exercise every batch family, compare values with scalar calls, validate shapes, and verify the runtime exports against the typed stub. The companion pykep-lambert project remains the end-to-end performance comparison for the SpOC 4–motivated propagation/Lambert workload.

Python migration from kep3

The native package is named pykep_rust, so it can be installed beside pykep/kep3 without import collisions. Phase 16 deliberately does not add a compatibility facade: a partial set of legacy spellings would hide important differences in ownership, dynamics construction, and error behavior. A separate facade can be considered after the deferred ecosystem modules have clear support decisions.

The comparison below is against the public Python surface of the pinned upstream version recorded in UPSTREAM_NOTICE.md.

Equivalent capabilities with descriptive names

Intentionally different

Deferred or unsupported

The following upstream areas are not part of the completed native numerical core and are not emulated:

• SPICE kernels, TLE parsing, and Python-defined UDPLA providers (PY-EXT-001);
• plotting helpers, trajectory-optimization UDPs, gym utilities, and optional ecosystem integrations (PY-ECOSYSTEM-001);
• symbolic heyoka dynamics/integrator construction and arbitrary user-supplied expression graphs (TA-SYMBOLIC-001);
• mima_from_hop and mima2_from_hop (MIMA-HOP-001);
• deprecated compatibility aliases and nested namespace layouts.

Use has_acceleration() and the VSOP2013 availability/threshold queries when code depends on an optional provider capability. Unsupported provider operations raise UnsupportedCapabilityError; unsupported ecosystem modules are absent rather than silently approximated.

Implementation status

Evidence-backed status as of 2026-07-26:

“Implemented” means the public contract is documented, validation is explicit, the committed C++ golden data passes except for documented numerical improvements, independent properties pass, and Rust/Python tests call the same core implementation. It does not imply that later modules exist.

The Python wheel uses the collision-safe pykep_rust import, ships a complete stub and py.typed, and has no C++ runtime dependency. Clean-wheel CI covers CPython 3.11–3.13 on Linux, macOS, and Windows. The upstream migration matrix records renames, deliberate contract changes, deferrals, and unsupported ecosystem modules.

The runnable example matrix covers every major public module in Rust and through the installed Python extension. Each example states units, expected behavior, runtime orientation, and required features; CI compiles all Rust examples and executes all Python scripts.

Phase 18 is a reproducible release candidate. The core is independently packageable, the binding implementation remains unpublished, performance regression/Miri/fuzz/Valgrind gates are maintained, and local crate/wheel/sdist smokes are required. Publication, a permanent name freeze, registry-download testing, and tagging remain explicitly blocked on release-owner metadata, external API review, trusted publishing, and authority.

Source map

The port baseline is pykep/kep3 3.0.1 at commit 53b1ca3ce5f8c223f96819b2ea9ba16c3719e63e. A checked box means the Rust module, C++ golden parity, independent validation, Python binding, and documentation required by the definition of done are complete.

Header-only numerical sources

• include/kep3/core_astro/constants.hpp → constants (Phase 2)
• include/kep3/core_astro/convert_julian_dates.hpp → time::julian (Phase 2)
• include/kep3/core_astro/kepler_equations.hpp → math::kepler_equations (Phase 2)
• include/kep3/core_astro/special_functions.hpp → math::stumpff (Phase 2)
• include/kep3/core_astro/convert_anomalies.hpp → astro::anomalies (Phase 3)

Translation units

• src/linalg.cpp → math::linalg (Phase 2)
• src/epoch.cpp → time::epoch (Phase 3)
• src/core_astro/ic2par2ic.cpp → astro::elements::classical (Phase 4)
• src/core_astro/mee2par2mee.cpp → astro::elements::equinoctial (Phase 4)
• src/core_astro/ic2mee2ic.cpp → astro::elements::equinoctial (Phase 4)
• src/core_astro/propagate_lagrangian.cpp → astro::propagation::lagrangian (Phase 5)
• src/core_astro/stm.cpp → astro::propagation::stm (Phase 5)
• src/core_astro/basic_transfers.cpp → astro::transfers::basic (Phase 6)
• src/core_astro/encodings.cpp → astro::encodings (Phase 6)
• src/core_astro/flyby.cpp → astro::flyby (Phase 6)
• src/lambert_problem.cpp → astro::lambert (Phase 6)
• src/core_astro/mima.cpp → astro::mima (Phase 6)
• src/planet.cpp → ephemeris (Phase 7)
• src/udpla/keplerian.cpp → ephemeris::keplerian (Phase 7)
• src/udpla/jpl_lp.cpp → ephemeris::jpl_lp (Phase 8)
• src/udpla/vsop2013.cpp → ephemeris::vsop2013 (Phase 9)
• heyoka integration requirements → integration facade and ADR 0004 (Phase 10; infrastructure rather than a source translation)
• src/ta/kep.cpp → dynamics::KeplerDynamics (Phase 11)
• src/ta/cr3bp.cpp → dynamics::Cr3bpDynamics (Phase 11)
• src/ta/bcp.cpp → dynamics::BcpDynamics (Phase 11)
• src/ta/zoh_kep.cpp → dynamics::zoh::ZohKeplerDynamics (Phase 12)
• src/ta/zoh_cr3bp.cpp → dynamics::zoh::ZohCr3bpDynamics (Phase 12)
• src/ta/zoh_eq.cpp → dynamics::zoh::ZohEquinoctialDynamics (Phase 12)
• src/ta/zoh_ss.cpp → dynamics::zoh::ZohSolarSailDynamics (Phase 12)
• src/ta/pontryagin_cartesian.cpp → dynamics::pontryagin::{CartesianMassOptimal, CartesianTimeOptimal} (Phase 13)
• src/ta/pontryagin_equinoctial.cpp → dynamics::pontryagin::{EquinoctialMassOptimal, EquinoctialTimeOptimal} (Phase 13)
• src/leg/sf_checks.cpp → leg::sims_flanagan validation (Phase 14)
• src/leg/sims_flanagan.cpp → leg::sims_flanagan (Phase 14)
• src/leg/sims_flanagan_alpha.cpp → leg::sims_flanagan::SimsFlanaganAlphaLeg (Phase 14)
• src/leg/zoh.cpp → leg::zoh (Phase 15)

The C++-specific visibility, serialization, and type-erasure support headers are reviewed for semantics but are not port targets.

Explicitly unavailable upstream ecosystem areas

These are not unchecked source-map rows. They are technically outside the native numerical-core product and have stable internal tracking identifiers:

• PY-EXT-001: SPICE kernels, TLE parsing, and Python-defined UDPLA providers require external data/runtime and dynamic Python callback contracts that are absent from the C/C++-free core.
• TA-SYMBOLIC-001: arbitrary heyoka expression graphs and Taylor-integrator objects cannot be reproduced without exposing a symbolic/JIT runtime; native evaluated dynamics and DOP853 propagation are the supported contract.
• PY-ECOSYSTEM-001: plotting, trajectory-optimization UDPs, gym utilities, and third-party integrations belong to their Python ecosystems rather than the numerical core.
• MIMA-HOP-001: mima_from_hop and mima2_from_hop depend on an upstream higher-order-propagation object that the native API intentionally does not expose. The Python migration matrix gives user-visible alternatives. A future change must resolve the corresponding tracking item and add a source-map row before claiming support.

Numerical validation

The evidence hierarchy for each numerical API is:

1. direct unit cases and invalid-input tests;
2. golden values generated by the pinned C++ implementation;
3. independent mathematical properties or external reference data;
4. Rust/Python cross-interface parity;
5. equivalent release-mode benchmarks.

Golden files use schema-versioned JSON. Floating-point values are encoded as C99 hexadecimal strings, with NaN, +Infinity, and -Infinity used only when recording upstream non-finite behavior. This preserves every finite binary64 value without relying on a JSON parser’s decimal conversion.

The initial C++ baseline was built in release mode with GCC 14.3.0. All 29 upstream C++ test executables passed on 2026-07-25. Benchmark executables were built but their results are not presented as Rust comparisons until equivalent Rust workloads exist.

Randomized oracle cases use a named PCG32 seed recorded in each data file. Discovered failures are promoted to fixed regression cases.

Phase 3 adds 9 epoch cases and 134 anomaly cases from the pinned C++ implementation. The anomaly set includes every conversion direction, elliptic and hyperbolic boundary regimes, angles outside one revolution, and 64 deterministic PCG32 solver samples. The Rust layer additionally rejects invalid calendars and true anomalies outside the physical hyperbolic asymptote instead of propagating non-finite values.

Phase 4 adds 206 direct element/Cartesian cases and 16 analytic Jacobians from the pinned implementation. Jacobian files explicitly record row-major output-by-input order. Independent tests cover 2,000 deterministic elliptic and hyperbolic state round trips, both equinoctial pole conventions, finite-difference derivatives, and inverse-Jacobian identities. NumPy batches are compared row-for-row with the scalar Python API.

Phase 5 adds 28 propagation cases spanning zero and negative duration, circular, elliptic, hyperbolic, near-parabolic, and many-period trajectories. Each case records both Lagrange-coefficient and universal-variable states plus the Lagrangian and Reynolds 6 by 6 STMs. The source file SHA-256 is 63e2340aa8e4f65a4ae831b7e1c8eb3a28930e2e1ed336d9324e0dbf44469e57. Independent tests check energy and angular-momentum conservation, time reversal, central finite differences, and STM composition. Scalar propagation uses only fixed-size stack values; the core path performs zero heap allocations. Batch APIs allocate exactly their returned output storage after copying Python-owned input before releasing the GIL.

Phase 6 records transfer, encoding, flyby, MIMA, and 13 Lambert solutions across three geometries in phase6-v1.json (SHA-256 d2a81311264ecc94a6caee15854273abb4df5f32f682fe5340ffeb247e3dd3e3). Lambert ordering is zero revolution followed by left/right pairs. Every returned branch is independently propagated to its requested endpoint; encoding pairs round-trip, and the flyby Jacobian is checked by central differences. MIMA2 is checked against the published upstream reference case.

Phase 7 adds six Keplerian states over negative, reference, near-reference, and long-span MJD2000 epochs. phase7-v1.json has SHA-256 7eda8fb03796ab7d70cdb0e94c422773422ed682aff9310cd06512c0f6396a54. Independent tests cover period recurrence, all supported element representations, ordered batches, explicit unsupported capabilities, and concurrent read-only evaluation through shared ownership.

Phase 8 adds true-anomaly elements, Cartesian states, and physical metadata for all eight JPL low-precision bodies at five epochs spanning the open 1800–2050 validity interval. phase8-v1.json contains 40 cases and has SHA-256 0a3408893b5c04fdfddb452408057faa92d38a46775290c32eb6d2393683e3da. Independent tests cover case-insensitive lookup, the exact supported-name set, ordered batches, safe-radius validation, and both excluded interval boundaries. The provider retains the source table’s slightly negative fitted Earth inclination near the ends of the interval while keeping the general public classical-element converter’s canonical [0, π] inclination contract.

Phase 9 adds 54 VSOP2013 states covering all nine bodies at six epochs around the 1890–2000 fit interval and the J2000/MJD2000 half-day offset. Two additional cases verify coefficient selection at the default 1e-5 and coarse 0.5 thresholds. phase9-v1.json has SHA-256 9d5af01df18acb17fcd1b2356af6f2cc08f1cf75e10864270d0c30732f8b00d8. The oracle uses pykep 3.0.1 linked to heyoka 7.10.0. At the embedded high-precision floor of 1e-9, the maximum observed Rust/C++ difference is 0.185 m in position and 6.6e-8 m/s in velocity. Independent checks cover case-insensitive names, feature reporting, threshold errors, clone determinism, ordered Python batches, and a build/test run without the optional coefficient feature.

Phase 10 selects a pure-Rust DOP853 backend through a pykep-owned facade. Decision tests compare Kepler states and the 6 by 6 STM to the independent analytic propagator, verify a parameter-sensitivity column by central differences, bound 100-orbit energy drift, and exercise a CR3BP close approach with Jacobi drift and backward reversal checks. Separate tests cover rejected steps, step exhaustion, bit-for-bit repeated solves, dense interpolation, terminal root location, malformed grids, non-finite values, and physical singularities. The selected final-state path retains no internal trajectory and performs no per-step heap allocation for fixed-size model states.

The candidate harness and matching C++ benchmark use the same six-state Kepler initial condition, final time, and 1e-12 scalar tolerances. Nominal and variational timings, allocation limitations, and the dense-output maximum-step caveat are recorded in ADR 0004. Those tests validate integration machinery independently of the production models added in Phase 11.

Phase 11 adds five sampled states and the final 6 by 6 STM for each of Kepler, CR3BP, and BCP in phase11-v1.json. The BCP case uses a nonzero Sun mass and a nonzero initial epoch, while the CR3BP case reproduces the representative upstream trajectory. The file has SHA-256 1c2b67eb203da62baf921db9f811b0ad0f763cbaa03f0ca80d6319870b0da807. The oracle uses pykep 3.0.1 and heyoka 7.10.0 at a requested tolerance of 1e-16; Rust comparison settings and achieved tolerances are documented in dynamics.md. Separate tests evaluate the source equations directly, verify the triangular equilibrium, preserve the Jacobi constant, make zero-Sun BCP reduce to CR3BP, check every state/parameter Jacobian column by central differences, and distinguish body singularities from solver failures.

Phase 12 records the final state and all first-order state/control variations for one upstream regression case from each of the four ZOH systems. phase12-v1.json has SHA-256 294bab93355628cd22991b83563c6f93deefc19dd93911dfc4b654a101764f22. Single-segment nominal tolerances range from 2e-12 to 3e-10; variation tolerances range from 2e-7 to 2e-5, reflecting the fixed-size numerical Jacobians documented in zero-order-hold.md. Independent tests cover exact switch ownership, malformed grids, zero-control reductions, manual segment-by-segment equivalence, forward/backward reversal, and activation of per-segment sensitivity columns.

Phase 13 records mass- and time-optimal Cartesian and modified-equinoctial states and all first-order variations with respect to the seven initial costates and upstream lambda0 variational argument. It also records the upstream dimensional 100-day equinoctial case. phase13-v1.json has SHA-256 a1276c4c35c7ad60b481d81d69f8df64b542bad1d5cc6571738db820cb0e2c3d. Normalized nominal trajectories use a scaled 3e-10 bound, the dimensional case uses 3e-9, and variations use 2e-4 to account for the generic centered numerical Jacobians. Independent tests cover minimized-Hamiltonian conservation, propagated central differences, output-by-input Jacobian orientation, explicit zero-primer errors, control normalization, and canonical Hamiltonian agreement across the analytic coordinate/costate transform.

Phase 14 records fixed-duration Sims–Flanagan mismatch and throttle constraints, all three analytic mismatch-Jacobian groups, and throttle Jacobians for one physical five-segment case and normalized four-segment cases at cuts zero, one half, and one. It also records equal and irregular direct durations for the alpha variant. phase14-v1.json has SHA-256 0bd6adddc72d850f1de5e4ab95a436128d9b8181f8de2b3bbbc67300e004d542. The Rust analytic gradients agree with the pinned C++ values to scaled bounds of 2e-10 to 3e-10 and with independent scale-adjusted central differences. Separate tests cover one segment, odd/even splits, zero and unit-limit throttle, all cut boundaries, normalized weights, invalid propulsion/gravity/mass values, dimension mismatches, and non-finite inputs.

Phase 15 records generic ZOH-leg mismatches and all four Jacobian groups for Kepler, CR3BP, modified-equinoctial, and ideal solar-sail dynamics. phase15-v1.json has SHA-256 9d8d424c2af10ceee497f74f62acb30c53924f78f7be8a528dbe44bf14767935. Nominal mismatches agree with the pinned C++/heyoka oracle within a scaled 3e-9 bound; endpoint, chronological-control, and time-grid derivatives agree within 3e-5, reflecting the fixed-size numerical model Jacobians. Every Kepler derivative column is also checked against independent scale-adjusted central differences. Separate tests cover cut zero and one, strict time grids, dimensions, finite values, solver options, state histories, ordered batches, and maximum-step exhaustion with direction, segment index, and time interval in the reported failure.

Phase 16 audits the complete Python surface against the pinned upstream exports and records every equivalent, renamed, intentionally different, deferred, and unsupported area in python-migration.md. Pytest checks every runtime export against the shipped stub, including parameter names/order, defaults, return-annotation presence, and class-member documentation, then exercises strided/read-only arrays, wrong dtype and rank, wrong shape, NaN/infinity, typed exceptions, repeatability, copied constructor inputs, and concurrent reuse. Mypy checks representative code against the packaged stub. A clean wheel matrix installs and imports CPython 3.11–3.13 wheels on Linux, macOS, and Windows; local shared-library inspection verifies that the extension has no C++ runtime dependency.

Phase 17 adds deterministic Rust/Python pairs for epochs and anomalies, elements and propagation, Lambert branches, ephemerides, gravity assists, Sims–Flanagan gradients, and CR3BP/ZOH dynamics, plus a NumPy batch example. Every Rust binary is compiled under workspace test and clippy gates and run locally in release mode. Pytest discovers and runs every Python script against the installed release extension. Ten Rustdoc examples cover every major module landing page and representative Lambert/ZOH types and compile with warnings denied.

Phase 18 adds a protocol-matched 100-sample Rust/C++ performance distribution, coarse CI regression thresholds, batch scaling, and allocation/cache/ vectorization profiling before any optimization. Miri runs suitable structural and parser tests, Valgrind checks the release harness, and bounded libFuzzer/AddressSanitizer campaigns cover epoch parsing, element conversion, Lambert inputs, and Reynolds-STM overflow boundaries. stabilization.md records exact results, limits, tool constraints, and the absence of algorithm changes.

Performance methodology

Benchmarks use release mode and Criterion. They contain no file or console I/O inside the measured loop. Raw Criterion state is build output under target/ and is not committed.

An orientation run on 2026-07-25 used Rust 1.97.1 on an AMD Ryzen 9 9950X under Linux 6.8.0-136. One Criterion process reported:

The Phase 5 orientation run used identical scalar inputs in Rust and the pinned C++ oracle, with both built in release mode on the same machine:

Phase 6 Rust medians from the same orientation run were 6.258 ns for Hohmann, 10.923 ns for flyby constraints, 34.567 ns for flyby delta-v, 235.76 ns for a zero-revolution Lambert problem, and 1.263 µs for the seven-solution multi-revolution case. Keplerian ephemeris evaluation measured 80.299 ns for one scalar epoch and 22.575 µs for an ordered 256-epoch batch in the Phase 7 orientation run. JPL low-precision Earth evaluation measured 95.198 ns for one scalar epoch and 31.611 µs for an ordered 256-epoch batch in the Phase 8 orientation run. The Phase 9 VSOP2013 spike measured 11.615 µs default-threshold initialization, 339.98 ns per default scalar state, 89.554 µs per ordered 256-state batch, and 37.157 µs per 1e-9 scalar state. The matching C++/heyoka orientation harness measured 63.96 ms and 552.53 ms first-time JIT initialization at 1e-5 and 1e-9, then 166 ns and 5.908 µs per scalar call. The feature increased the release benchmark executable from 2.6 MiB to 7.0 MiB. See ADR 0003 for the data/cache decision and caveats. The Phase 10 integration decision harness measured the selected DOP853 facade at 11.865 µs for the representative nominal Kepler solve and 85.663 µs for the state plus 6 by 6 STM. The equivalent warmed C++/heyoka solves measured 3.722 µs and 84.440 µs; cloning the cached C++ nominal integrator once cost 0.467 ms. A same-profile candidate tool measured 15.431 µs for the selected facade and 7.369 µs for ode_solvers, whose missing root/sensitivity facilities and per-step allocations outweighed the nominal timing advantage. See ADR 0004 for configuration, ranges, and risks. The Phase 11 evaluated-model orientation measured Kepler, CR3BP, and BCP right-hand sides at 10.828 ns, 32.833 ns, and 63.386 ns. A representative CR3BP propagation measured 7.808 µs and its state-plus-STM propagation measured 45.368 µs. The same initial state, final time, parameter, and 1e-12 tolerance in warmed C++/heyoka measured 2.885 µs and 95.017 µs median, respectively. The Phase 12 ZOH Kepler RHS measured 8.928 ns. A 32-segment alternating control schedule measured 36.089 µs in Rust and 10.268 µs in the warmed C++/heyoka harness with the same state, boundaries, controls, and 1e-12 tolerance. The Rust path deliberately starts an independent DOP853 solve at each switch; the timing includes those 32 restarts and confirms there is no segment-count-dependent control search inside integration. The Phase 13 Cartesian mass-optimal RHS measured 52.384 ns and its normalized 1.2345-time-unit propagation measured 142.52 µs in Rust. The warmed C++/heyoka integrator measured a 11.716 µs median for the same state, parameters, final time, and 1e-12 tolerance. Rust uses a 0.01 maximum step to enforce the recorded oracle tolerance; the C++ Taylor solve has no equivalent maximum-step restriction, so this orientation identifies an integration-performance target rather than like-for-like algorithm speed. The Phase 14 physical five-segment Sims–Flanagan case measured 1.002 µs for mismatch evaluation and 4.241 µs for its complete analytic mismatch Jacobian in Rust. The same warmed C++ configuration measured 1.037 µs and 12.748 µs, respectively. Both sides used the same endpoint states, masses, controls, duration, propulsion parameters, gravity parameter, and cut = 0.6; neither timing includes construction or validation. The Phase 15 20-segment normalized Kepler ZOH leg measured a Criterion 23.67 µs median for mismatch evaluation and 493.55 µs for the complete endpoint/control/time-grid Jacobian. The matching warmed C++/heyoka harness averaged 6.865 µs and 182.17 µs. Both use the same states, chronological controls, time grid, constants, cut, and 1e-12 tolerance, with no maximum step. Rust integrates each segment independently and computes fixed-size numerical dynamics Jacobians, making both paths visible optimization targets. The Phase 16 release-wheel wrapper harness used 20,000 items and the median of nine samples on the same machine. A Python scalar loop around stumpff_c measured 38.9 ns/item versus 14.6 ns/item through stumpff_c_batch, a 2.67× throughput improvement. Scalar Lagrange calls measured 0.72 µs/item versus 0.09 µs/item for the N × 6 NumPy batch, a 7.82× improvement. These include Python input/output conversion and demonstrate that throughput-sensitive code should use explicit batches; they are not Rust-core timings.

These are not cross-language speed claims. CPU frequency was not fixed and the run is not a substitute for distributions collected under controlled affinity and power settings. The Julian arithmetic result is small enough that compiler optimization and timer resolution dominate its interpretation. The C++ column is an elapsed average from five million calls rather than a Criterion distribution; it is included as the required same-input orientation baseline, not as a statistically controlled language comparison.

Run the maintained harness with:

cargo bench -p pykep-core --bench foundation
cargo bench -p pykep-core --bench elements
cargo bench -p pykep-core --bench propagation
cargo bench -p pykep-core --bench mission
cargo bench -p pykep-core --bench integration
cargo bench -p pykep-core --bench dynamics
cargo bench -p pykep-core --bench legs
python python/benchmarks/wrapper_overhead.py
cargo run --release -p pykep-lambert-optimization-benchmark

The same harness includes scalar elliptic/hyperbolic anomaly solvers and a 64-value batch-equivalent loop. This keeps branch-heavy iterative work separate from the foundation arithmetic measurements. Element benchmarks separate scalar classical/equinoctial conversion, analytic Jacobian evaluation, and a 64-state batch-equivalent loop. Propagation benchmarks separate elliptic and hyperbolic Lagrange coefficients, universal variables, analytic STM evaluation, and a 1,024-state throughput loop. Scalar propagation and STM APIs operate entirely on fixed-size arrays and perform no heap allocation. Mission benchmarks cover basic transfers, flyby constraints/delta-v, and both zero- and multi-revolution Lambert solution construction. The same harness measures scalar and 256-epoch Keplerian and JPL low-precision ephemeris evaluation separately, plus VSOP2013 initialization, scalar, high-precision, and batch paths. Integration benchmarks separate nominal six-state DOP853 propagation from an augmented state plus STM solve. The final-state output callback does not retain internal steps. Dynamics benchmarks separate raw evaluated right-hand sides from CR3BP nominal and variational propagation, ZOH schedules, and Pontryagin state/costate propagation. Leg benchmarks separate Sims–Flanagan mismatch and analytic-gradient evaluation and generic ZOH mismatch and sensitivity evaluation for the same configurations used by their C++ harnesses.

The standalone Lambert optimization benchmark ports the fixed easy.kttsp leg from pykep-lambert. It measures deterministic objective throughput and then applies the native fcmaes-core CMA-ES and BiteOpt implementations to wait time and time of flight. Its source revision, physical constants, decision bounds, penalty, optimizer budget, and seed are printed with every run; see tools/lambert-optimization-benchmark/README.md for the complete protocol.

C++ comparisons are added only when both sides execute identical input data, validation policy, branch families, tolerances, and output work. Initialization and batch throughput are reported separately from warm scalar latency.

Phase 18 adds a protocol-matched 100-sample Rust/C++ distribution, bootstrap median confidence intervals, CI regression limits, allocation/cache/ vectorization profiles, and five-point Python batch scaling. Full results and environment limitations are in stabilization.md.

Stabilization and release-candidate evidence

This document records the Phase 18 checks that complement numerical parity and invariant testing. It does not claim that registry publication, cross-platform hosted CI, external API review, or a release tag has occurred.

Matched Rust/C++ distribution

The representative five-segment Sims–Flanagan case uses byte-for-byte equivalent inputs and output work in the Rust release tool and an external C++ oracle harness. Both ran sequentially on CPU 0 with one thread, a three-second warmup per workload, 100 samples, 10,000 mismatch calls/sample, and 1,000 gradient calls/sample. Rust used 1.97.1 with the workspace release profile (opt-level=3, thin LTO, one codegen unit); C++ used GCC 14.3.0 and CMake Release. The host is an AMD Ryzen 9 9950X under Linux 6.8.0-136.

The host governor reported powersave; CPU frequency was not fixed. Hardware performance counters were unavailable because perf_event_paranoid=4. Accordingly, these are same-session orientation distributions, not universal language-speed claims. The deterministic analysis script reports a percentile-bootstrap 95% interval for the median:

Reproduce the Rust samples and table with:

cargo run --release -p pykep-release-benchmark > rust.csv
python tools/analyze_benchmark.py Rust=rust.csv C++=cpp.csv

The C++ oracle remains outside the public repository by policy.

Regression policy

CI runs an 11-sample/100-ms-warmup smoke protocol. Median limits are 10 µs for mismatch and 45 µs for the gradient, about 9.9× and 10.6× the local baselines. They intentionally catch order-of-magnitude regressions while leaving room for shared-runner scheduling and CPU differences. They are not used for cross-language assertions or small optimization claims.

Allocation, cache, vectorization, and scaling

Valgrind 3.22.0 found zero memory errors and zero definitely/indirectly/ possibly lost bytes; one 544-byte runtime block remained reachable at exit. DHAT quick-protocol aggregates were:

The gradient returns dynamic output matrices and is the clear allocation target. Cachegrind over both quick workloads recorded 335,788,370 instructions, 133,876,521 data references, 3,775 L1 data misses, 5,590 last-level misses, and a 2.2% branch-mispredict rate. LLVM loop-vectorization remarks and emitted IR showed no packed-double vector loop in this workload. No algorithm was changed: the active parity and invariant suite is green, and the profile points to API output storage rather than numerical formulas.

Release-wheel Lagrange batch scaling (median of 11) was:

Dynamic analysis and fuzzing

• Miri nightly passed selected fixed-layout, derivative-consistency, error-formatting, and epoch-parser tests. Native tests remain authoritative for exact floating-point references because Miri software floating point can differ by a few ULPs.
• Valgrind Memcheck passed the release benchmark with zero errors.
• Ten-second libFuzzer/AddressSanitizer campaigns completed without a crash: 13,016,052 epoch-parser inputs, 10,115,059 element-conversion inputs, and 4,588,894 Lambert inputs. The 2026-07-26 review follow-up added a Reynolds-STM overflow target, which completed 7,197,911 inputs without a crash.
• RustSec and cargo-deny pass. paste 1.0.15 remains an allowed unmaintained transitive dependency under RUSTSEC-2024-0436, documented in ADR 0004.

The fuzz harness is excluded from released artifacts. Its native compiler and sanitizer are development tools, not core build/runtime dependencies.

Public API documentation

The Rust core denies missing documentation for public items, and rustdoc builds with warnings denied. The Python audit inventories all 145 exported names; each is represented in the checked type stub and directly exercised by the integration suite. Public module, object, method, and descriptor docstrings are checked alongside that inventory. The stub audit compares runtime parameter names, order, and defaults and requires return annotations.

The 2026-07-26 review follow-up measured 376/376 documented Rust items (100.0%) and ten compiled examples, up from one. Rustdoc’s per-item example metric is 5.5% because examples are intentionally placed on major module and type landing pages rather than repeated on every field and accessor. Core coverage is 91.01% regions, 92.26% functions, and 90.26% lines. The formerly weak Lambert module rose from 74.86% to 88.07% line coverage; ephemeris/mod.rs is 97.48% and error.rs is 100%.

The migration matrix records every unavailable upstream name with a technical reason and tracking identifier. External review is still required before declaring these public names frozen.

Local artifact consumption

cargo package produced and verified an 80-file, 2.3 MiB compressed pykep-core-0.1.0.crate. Its normal/build graph contains only Rust crates and no cc/CMake/native-library dependency. The extracted archive—not the workspace source—compiled and ran from a fresh Cargo binary project.

Maturin produced a 2.9 MiB CPython 3.12 manylinux wheel and a 2.4 MiB sdist. The wheel and the wheel built from the sdist both imported from separate fresh virtual environments. The packaged stub and py.typed marker were present. ldd reported only libgcc_s, libm, libc, and the ELF loader; there was no C++/kep3/Boost/heyoka runtime. This is the local Linux result. The hosted matrix is configured to build and consume wheels for CPython 3.11–3.13 on Linux, macOS, and Windows, but is not represented here as an executed remote run.

External blockers

Public names are documented but not declared permanently frozen: external API review is still required. The independent repository URL and owner are now recorded, while private security contact, registry-side trusted publishers, registry uploads, downloads of published artifacts, and the release tag all require release-owner authority. See RELEASE.md for the exact boundary.

Development

The workspace MSRV is Rust 1.88.0. The normal local quality gate is:

cargo fmt --all -- --check
cargo clippy --workspace --all-targets --all-features -- -D warnings
cargo test --workspace --all-features --locked
RUSTDOCFLAGS="-D warnings" cargo doc --workspace --all-features --no-deps
cargo +1.88.0 check --workspace --locked

Build the GitHub Pages documentation with the same mdBook 0.5.4 release used by .github/workflows/docs.yml:

python tools/check_markdown_links.py
mdbook build

Run benchmarks separately so timing work is never hidden in the test suite:

cargo bench -p pykep-core
cargo run --release -p pykep-release-benchmark -- --quick --check
python python/benchmarks/wrapper_overhead.py --scaling

The standalone Phase 10 candidate comparison is reproducible with:

cargo run --release --manifest-path tools/phase10-candidates/Cargo.toml

For Python integration:

python -m venv .venv
.venv/bin/python -m pip install --upgrade pip
.venv/bin/python -m pip install "maturin[patchelf]>=1.7,<2" pytest
env -u CONDA_PREFIX VIRTUAL_ENV="$PWD/.venv" \
  PATH="$PWD/.venv/bin:$PATH" \
  .venv/bin/maturin develop --release \
  --manifest-path crates/pykep-py/Cargo.toml
.venv/bin/python -m pytest

Coverage:

cargo llvm-cov -p pykep-core --all-features --all-targets --summary-only
cargo llvm-cov --workspace --all-targets --summary-only

Release-candidate dynamic checks:

cargo +nightly miri test -p pykep-core --lib --no-default-features fixed_matrix_operations
cargo +nightly fuzz run epoch_parser -- -max_total_time=30
cargo +nightly fuzz run element_conversions -- -max_total_time=30
cargo +nightly fuzz run lambert_inputs -- -max_total_time=30
cargo +nightly fuzz run reynolds_stm -- -max_total_time=30
valgrind --error-exitcode=1 --leak-check=full \
  target/release/pykep-release-benchmark --quick

The full protocol, profiler commands/results, and Miri floating-point scope are documented in stabilization.md. Packaging and clean-artifact consumption are documented in RELEASE.md.

Build state belongs in target/ or .venv/ and is ignored. Development-only C++ oracle tools and internal planning notes are not part of this standalone repository.