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 records current parity and independent checks; the stabilization evidence preserves the historical pre-0.1.0 coverage, performance, and release-candidate snapshot.

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, taylor-integration.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: historical pre-0.1.0 matched distributions, profiling, regression limits, Miri/fuzz/Valgrind evidence, and resolved release blockers;
• status.md, source-map.md, and validation.md: limitations, provenance, and validation evidence;
• development.md: contributor commands and repository policy;
• add-ode-system.md: complete implementation, Taylor, test, benchmark, Python, and documentation checklist for a new dynamics model.

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.

Complete trajectory-optimization tutorial

For an end-to-end application, see the GTOC1 “Save the Earth” Rust tutorial. It combines pykep-core ephemerides, propagation, multi-revolution Lambert solutions, gravity assists, and Sims–Flanagan legs with fcmaes-core optimization. The tutorial covers a reproducible fixed planet sequence, multi-fidelity sequence discovery, split-brain planet-order search, strict validation, and the distinction between a VSOP2013 model score and an official GTOC1 result.

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. Nominal propagate uses the fixed-system Taylor backend, matching upstream pykep’s ta families; propagate_with_method can select DOP853 explicitly. State-transition matrices are row-major, with output state components in rows and initial state components in columns. propagate_with_stm retains the generic DOP853 direct-variational sensitivity contract; Taylor sensitivities currently use centered complete propagations and remain available through propagate_with_stm_method. See High-accuracy Taylor integration.

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.

High-accuracy Taylor integration

pykep-core has two pure-Rust adaptive integration backends:

• Taylor is the default for the eleven built-in dynamics types, matching upstream pykep’s ta model families. It is most useful near 1e-14 and below, or when long-term invariant drift matters more than minimum low-accuracy wall time.
• Dop853 accepts arbitrary DynamicsModel implementations, supports event location, and remains the general-purpose choice.

This is a focused implementation of Taylor-series integration for pykep’s fixed systems. It is not a port of heyoka’s symbolic engine, LLVM compiler, event machinery, batch mode, or arbitrary-precision support.

What the backend computes

For a state expanded around the current time,

\[ x(t+h)=\sum_{k=0}^{p}x_kh^k, \]

the model recurrence constructs coefficients through order p. Products use Cauchy convolution; real powers, square roots, exponentials, sine, and cosine use coefficient recurrences. The controller:

1. selects order 8–24 from the requested tolerance;
2. scales the final two coefficients by atol + rtol * max(abs(state), 1);
3. chooses a conservative step from both coefficient estimates;
4. applies initial and maximum step caps;
5. clips the final step exactly to the requested time.

Backward propagation uses the same coefficients with a negative step. Every accepted state is checked for finite values and revalidated against the model domain. A zero or non-advancing floating-point step is reported as underflow.

The implementation follows the recurrence and step-selection ideas described by Jorba and Zou, while being written independently for fixed Rust arrays: A software package for the numerical integration of ODEs by means of high-order Taylor methods.

Selecting the backend

KeplerDynamics::propagate(), the corresponding CR3BP and BCP convenience methods, and AdaptiveIntegrator::default() select Taylor. The method can also be stated explicitly:

#![allow(unused)]
fn main() {
use pykep_core::dynamics::KeplerDynamics;
use pykep_core::integration::{IntegrationMethod, IntegratorOptions};

let result = KeplerDynamics.propagate_with_method(
    0.0,
    [0.8, -0.2, 0.1, 0.03, 1.0, 0.02],
    0.5,
    1.0,
    IntegratorOptions {
        relative_tolerance: 1e-14,
        absolute_tolerance: 1e-14,
        ..IntegratorOptions::default()
    },
    IntegrationMethod::Taylor,
)?;
assert_eq!(result.time, 0.5);
Ok::<(), pykep_core::PykepError>(())
}

The lower-level Taylor and AdaptiveIntegrator facades also provide final, dense-grid, and seeded-sensitivity propagation. ZOH schedules have parallel *_with_method functions.

Taylor event location is not implemented. Use Dop853::propagate_until_event() when a terminal crossing is required. Direct DOP853 propagation is also the supported path for arbitrary third-party DynamicsModel implementations.

Supported models

The nine upstream equation families produce eleven concrete Rust types:

Singular points stay errors. Examples include collisions, zero Pontryagin primer norm, non-positive mass, and a zero equinoctial radial denominator. The series backend does not smooth or silently cross these non-analytic model boundaries.

Measured crossover

The fixed benchmark uses the eccentric nondimensional state [0.5, 0, 0, 0, sqrt(3), 0], one revolution, release mode, and analytical Lagrange propagation as the error reference.

These are measurements on one development host, not portable promises. The 100- and 1,000-revolution rows and all raw counters are in data/taylor-kepler-benchmark.csv. Taylor coefficient sweeps and DOP853 RHS calls are different work units, so wall time at matched accuracy—not the work-counter ratio—is the meaningful comparison.

Reproduce the data with:

cargo run --release -p pykep-taylor-benchmark \
  > docs/data/taylor-kepler-benchmark.csv

Incremental coefficient evaluators

All eleven built-in models advance only the next required Taylor coefficient. Kepler and ZOH Kepler use compact hand-written recurrences. The other nine models lazily build immutable scalar-expression tapes, eliminate structural common subexpressions, and then reuse those tapes for every propagation. Trigonometric, exponential, explicit-time, and real-power nodes all have incremental coefficient recurrences.

The table records the eight final migrations against the previous repeated full-series evaluator at commit c9ecd17. Each row is an end-to-end propagation at 1e-12, pinned to one core, with 1,000 iterations per sample and the median of seven warmed samples:

These development-host measurements establish that every migration cleared the significance gate; they are not portable latency guarantees. Matching final-state checksums and the validation layers below guard the optimization. Criterion entries under taylor/* keep all eight paths visible in routine performance runs.

Sensitivities and validation boundaries

Taylor::propagate_with_sensitivities() currently uses centered differences of complete Taylor propagations, one plus/minus pair per seed direction. This is deterministic and matches the direct DOP853 variational result and heyoka Kepler STM within the declared 3e-7 test tolerance, but its cost grows as 2W + 1 propagations. For a wide STM, DOP853’s directly integrated variational equations are normally cheaper. The STM and ZOH sensitivity convenience methods retain DOP853 as their default for this reason; request Taylor explicitly when algorithm-family parity is more important than cost.

Validation has three independent layers:

• closed-form tests for series products, real powers, trigonometric functions, exponentials, and square roots;
• fixed DOP853 comparisons for all eleven concrete dynamics types, including all four Pontryagin variants;
• a development-only official-heyoka harness for Kepler, CR3BP, BCP, and a Kepler STM.

The heyoka fixture records version 7.10.1 and is regenerated with tools/heyoka-cross-validation/generate.py. It is intentionally not a normal CI dependency. The larger ZOH and Pontryagin families use the existing pinned upstream pykep fixtures plus DOP853 comparisons; they are not presented as direct heyoka cross-checks.

Why no LLVM or new crate?

The equations are already fixed and hand-written in pykep-core. A symbolic JIT would add compilation latency and a large dependency surface without changing the set of systems this crate needs. Keeping the series engine private also avoids publishing an extension trait before it has a real external consumer.

The four Pontryagin kernels use a smaller form of generation without LLVM. During lazy tape construction, the scalar Hamiltonian is differentiated in reverse graph order to generate the seven analytic costate expressions. Direction and throttle nodes are explicit stop-gradients, preserving the Pontryagin envelope convention used by the evaluated equations. For example, structural common-subexpression elimination reduces each Cartesian mass-optimal graph from 425–426 forward-dual operations to 153 operations. Runtime integration only evaluates immutable coefficient graphs; it does not repeat symbolic differentiation.

If another project later needs the arithmetic and the internal model contract survives unchanged, the module can be extracted into a workspace crate. That is an API and release decision, not a prerequisite for using Taylor today.

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.

The established propagate_schedule* functions continue to use DOP853. The corresponding *_with_method functions accept IntegrationMethod::Taylor for all four built-in models. Each segment still ends exactly at its switching boundary; backend selection does not change switch ownership.

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 DOP853 implementation does not repropagate every prior segment for every control.

Taylor’s current seeded-sensitivity implementation uses centered complete propagations and therefore costs 2W + 1 nominal solves per segment. Prefer DOP853 for wide ZOH Jacobians until direct Taylor variational recurrences are implemented.

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 integration and sensitivity contracts; it does not expose or depend on a symbolic-expression runtime. The four built-in types can be passed to the default Taylor AdaptiveIntegrator, to Taylor directly, or to Dop853 when events or direct variational propagation are needed.

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.

In the Taylor recurrence the minimizing control evolves with the complete state/costate series. When forming canonical costate rates, the control is held fixed with respect to the physical-state derivative, as required by the envelope theorem. Zero primer norm remains an explicit non-analytic error.

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.

mismatch_constraints() preserves the established DOP853 propagation. Built-in models additionally provide mismatch_constraints_with_method(IntegrationMethod), so repeated derivative-free objective evaluations can select the accelerated Taylor backend explicitly.

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. state_history_with_method(samples_per_segment, IntegrationMethod) selects DOP853 or Taylor for built-in models. DOP853 backward segments are evaluated through the equivalent increasing coordinate tau = segment_start - time to avoid a decreasing-time dense-output boundary issue in the numerical backend; Taylor supports the decreasing grid directly. 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.

The evaluated model propagators currently keep their established DOP853 backend in Python. Runtime selection of the new fixed-system Taylor backend is Rust-only in this release; Python does not expose an integrator object or a backend-name string whose compatibility would need to be frozen.

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-08-03:

“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.

The synchronized 0.1.4 release is published as pykep-core on crates.io and pykep-rust on PyPI. The internal pykep-py implementation crate remains publish = false by design. Tag-gated trusted-publishing workflows, clean crate/wheel/source-distribution consumption, docs.rs, and GitHub Pages are part of the release process. Performance regression, Miri, fuzz, Valgrind, dependency, MSRV, cross-platform wheel, and API-documentation checks remain maintained quality gates rather than claims of formal verification.

Version 0.1.4 includes incremental Taylor coefficient evaluation for all eleven built-in models and Rust-only method selection for nominal ZOH-leg mismatch/history calculations. Existing no-suffix ZOH-leg methods, Jacobians, and the Python ZOH-leg surface retain their DOP853 behavior. See the performance decision guide for the measured boundary between the ahead-of-time Rust backend and warmed upstream heyoka.

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)
• high-order Taylor recurrence/controller → integration::Taylor (independent implementation for the eleven fixed model types; not a symbolic-runtime 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 user-defined Taylor-integrator objects remain out of scope. Native evaluated dynamics support DOP853 and the eleven built-in types additionally support the fixed-system Taylor backend.
• 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-by-phase evidence

Phase 3: epochs and anomalies

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: elements and Jacobians

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: two-body propagation

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: mission-design utilities

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: Keplerian ephemerides

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: JPL low-precision ephemerides

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: VSOP2013 ephemerides

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: adaptive integration

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: evaluated dynamics

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: zero-order-hold dynamics

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: Pontryagin dynamics

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: Sims–Flanagan legs

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: generic ZOH legs

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: Python API

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: runnable examples

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: release stabilization

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.

Taylor extension

The Taylor extension adds closed-form recurrence tests, coefficient-level comparisons between the optimized incremental evaluators and an independent full-series reference through order 24, matched DOP853 trajectories for all eleven built-in model types, and a separate official heyoka 7.10.1 fixture for Kepler, CR3BP, BCP, and the Kepler STM.

taylor-heyoka-v1.json is regenerated by the optional tools/heyoka-cross-validation/generate.py harness. The committed one/100/1,000-revolution accuracy and timing protocol is data/taylor-kepler-benchmark.csv; analytical Lagrange propagation is its state-error reference. The ZOH and Pontryagin Taylor checks use their existing pinned upstream fixtures plus same-problem DOP853 comparisons and are not mislabelled as direct heyoka runs.

Downstream release rehearsals

Before changing the built-in nominal default to Taylor, two downstream applications were rehearsed against the local crate through a temporary Cargo registry patch:

• fcmaes-rust/tutorials/gtoc1 passed all five release tests and reproduced its stored score, mismatch, flyby margin, and sampled solar distance exactly. A disposable physical-units check propagated all seven stored Lambert arcs with analytical Lagrange, Taylor, and DOP853 propagation at rtol = 1e-12, atol = 1e-6. Taylor’s maximum normalized state error over the arcs was 1.1e-12, versus 3.8e-10 for DOP853.
• pykep-lambert passed all eleven release tests, including scalar/batch identity and every optimizer mode. Its standard 32,768-transfer workload retained identical checksums and the same feasible optimizer solution. Five warmed runs remained consistent with its published medians; this application uses analytical Keplerian ephemerides and Lambert solves, so the Taylor default does not alter its normal computation path.

Neither downstream manifest was changed during that pre-release rehearsal; both pinned pykep-core = "=0.1.2" at the time. They were subsequently moved to published versions and retested: the GTOC1 tutorial now pins 0.1.4, while pykep-lambert pins 0.1.3 because its analytical ephemeris/Lambert workload does not exercise the 0.1.4 Taylor and ZOH-leg 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:

Mission-design kernels

Phase 6 measured each operation separately in the same Rust orientation run:

Ephemerides

Scalar and ordered-batch measurements were kept separate:

The matching C++/heyoka VSOP2013 harness measured first-time JIT costs of 63.96 ms at 1e-5 and 552.53 ms at 1e-9. Warm scalar calls took 166 ns and 5.908 µs, respectively. Enabling VSOP2013 increased the Rust release benchmark executable from 2.6 MiB to 7.0 MiB. ADR 0003 explains the data and cache decision.

Integration and dynamics

DOP853 selection

The Phase 10 decision harness produced these orientation measurements:

Cloning the cached C++ nominal integrator once cost 0.467 ms. In a same-profile candidate test, the selected facade took 15.431 µs and ode_solvers took 7.369 µs. The latter lacked the required root and sensitivity facilities and allocated per step, so nominal timing alone did not decide the dependency. ADR 0004 records the configuration, ranges, and risks.

Taylor integration

The fixed-system Taylor backend uses a matched-accuracy protocol because coefficient sweeps are not comparable to DOP853 right-hand-side calls. For one eccentric nondimensional revolution, Taylor’s speed relative to DOP853 was:

Taylor produced smaller final-state errors at every point. This is deliberately a narrow conclusion: Taylor is a high-accuracy option, not a low-accuracy replacement. The committed CSV also contains the 100- and 1,000-revolution rows; see High-accuracy Taylor integration.

All eleven built-in Taylor models now use incremental coefficient evaluators. Across the representative eight-model migration benchmark, replacing repeated full-series evaluation reduced warmed end-to-end propagation time by 6.45× to 28.97×. The linked Taylor guide gives per-model results and validation limits.

Evaluated models and controlled dynamics

The Phase 11 comparisons use the same initial state, final time, parameter, and 1e-12 tolerance.

The Phase 12 Rust schedule deliberately starts an independent DOP853 solve at each switch. Its timing includes all 32 restarts and confirms that integration does not perform a segment-count-dependent control search.

The Phase 13 propagation covers 1.2345 normalized time units. Rust uses a 0.01 maximum step to meet the recorded oracle tolerance, whereas the C++ Taylor solve has no equivalent limit. Treat this as an identified performance target, not as a like-for-like algorithm comparison.

Low-thrust legs

The Phase 14 comparison uses the same endpoints, masses, controls, duration, propulsion parameters, gravity parameter, and cut = 0.6. Construction and validation are outside both timings.

The Phase 15 comparison uses the same states, chronological controls, time grid, constants, cut, and 1e-12 tolerance, with no maximum step. Rust integrates every segment independently and uses fixed-size numerical dynamics Jacobians. Both operations therefore remain visible optimization targets.

Python batch throughput

The Phase 16 release-wheel harness used 20,000 items and the median of nine samples:

These measurements include Python input/output conversion; they are not Rust-core timings. They show why throughput-sensitive Python code should use the explicit batch APIs.

Interpreting the results

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.

Reproducing the measurements

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
cargo run --release -p pykep-taylor-benchmark

Each benchmark group has a distinct scope:

The foundation anomaly loop keeps branch-heavy iterative work separate from arithmetic kernels. Scalar propagation and STM APIs use fixed-size arrays and perform no heap allocation. Ephemeris measurements distinguish initialization, scalar, high-precision, and batch paths. The integration final-state callback does not retain internal steps.

Standalone Lambert optimization

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.

Cross-language comparison policy

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.

Release stabilization evidence

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 is the historical Phase 18, pre-0.1.0 evidence snapshot that complements numerical parity and invariant testing. Its artifact versions and then-open blockers are intentionally preserved as the record of that release candidate; they do not describe the current release state. The synchronized 0.1.4 artifacts are now published. See Implementation status and the release procedure for the current state and procedure.

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.

Historical external blockers

At the time of this snapshot, public names were not declared permanently frozen and external API review was still required. The repository URL and owner had been recorded, while the private security contact, registry-side trusted publishers, registry uploads, download checks, and release tag still required release-owner authority. Those publication blockers were later resolved; the paragraph remains to explain what the Phase 18 evidence did and did not establish.

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

Build and verify the complete PDF edition after the HTML book is current:

python scripts/build_book_pdf.py

The PDF builder needs Chrome or Chromium plus mutool. It creates a navigable bookmark tree, keeps chapter links inside the book, converts links to repository artifacts into portable GitHub URLs, and writes docs/pykep-rust-book.pdf.

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.

For the model contract and the complete definition of done for a new dynamics family, see Adding an ODE system.

Adding an ODE system

This guide describes the complete path for adding a first-order ordinary differential equation to pykep-core, including the evaluated Rust model, sensitivities, the optional fixed-system Taylor backend, batches, Python bindings, tests, benchmarks, and documentation.

The central contract is

dy/dt = f(t, y, p)

where the state has a compile-time dimension N and the constant parameter vector has a compile-time dimension P.

Decide the scope first

There are two materially different extension levels.

External or DOP853-only model

Any downstream crate can implement DynamicsModel<N, P> and pass the model to Dop853. This is the stable public extension point. It supports nominal propagation, dense output, events, and, after implementing DifferentiableDynamicsModel, direct variational sensitivities.

An external model cannot implement Taylor support. TaylorCoefficientModel is private and TaylorDynamicsModel is deliberately sealed. This prevents the crate from promising a symbolic/coefficient API before it has a stable third-party contract.

Built-in pykep model

A model shipped by this repository normally needs the complete product surface:

• evaluated Rust right-hand side;
• domain validation and typed errors;
• state and parameter Jacobians when sensitivities are meaningful;
• Taylor coefficients and the TaylorDynamicsModel marker;
• convenience propagation methods;
• ordered scalar and parallel-batch APIs where users will evaluate many independent cases;
• Python scalar and batch bindings if the model is part of the Python product;
• independent validation, performance evidence, Rustdoc, mdBook documentation, source-map/status updates, and a changelog entry.

Do not add a model to the default AdaptiveIntegrator surface until its Taylor implementation exists and has been validated. A deliberately DOP853-only built-in model must use Dop853 directly and document that limitation.

1. Write down the numerical contract

Before coding, record:

1. The model name and source equation or upstream reference.
2. N, the state dimension, and the exact order of every state component.
3. P, the parameter dimension, and the exact order of every parameter.
4. Units, normalization, coordinate frame, epoch/time convention, and sign conventions.
5. The valid domain of every state and parameter.
6. Singular geometries and non-analytic switching surfaces.
7. Whether the equations depend explicitly on time.
8. Conserved quantities, symmetries, closed-form cases, equilibria, or reversible cases that can serve as independent checks.
9. Whether analytic state and parameter Jacobians are available.
10. Whether users need dense output, events, sensitivities, ZOH control, batches, or Python access.

Treat state and parameter ordering as API. Changing it later is a breaking change even if the Rust type remains [f64; N].

2. Add the evaluated Rust model

Place a substantial new family in crates/pykep-core/src/dynamics/<system>.rs and export the module from crates/pykep-core/src/dynamics.rs. A small closely related model can live beside its family.

Prefer a zero-sized, copyable model:

#![allow(unused)]
fn main() {
use pykep_core::integration::{DynamicsModel, Dop853, InitialValueProblem};
use pykep_core::integration::IntegratorOptions;
use pykep_core::{PykepError, Result};

#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)]
pub struct OscillatorDynamics;

impl DynamicsModel<2, 1> for OscillatorDynamics {
    const NAME: &'static str = "harmonic oscillator";

    fn validate(
        &self,
        time: f64,
        state: &[f64; 2],
        parameters: &[f64; 1],
    ) -> Result<()> {
        if !time.is_finite() {
            return Err(PykepError::NonFiniteInput { parameter: "time" });
        }
        if state.iter().any(|value| !value.is_finite()) {
            return Err(PykepError::NonFiniteInput { parameter: "state" });
        }
        if !parameters[0].is_finite() {
            return Err(PykepError::NonFiniteInput {
                parameter: "angular_frequency",
            });
        }
        if parameters[0] <= 0.0 {
            return Err(PykepError::InvalidInput {
                parameter: "angular_frequency",
                reason: "must be greater than zero".into(),
            });
        }
        Ok(())
    }

    fn rhs(
        &self,
        time: f64,
        state: &[f64; 2],
        parameters: &[f64; 1],
        derivative: &mut [f64; 2],
    ) -> Result<()> {
        self.validate(time, state, parameters)?;
        let omega = parameters[0];
        *derivative = [state[1], -omega * omega * state[0]];
        if derivative.iter().all(|value| value.is_finite()) {
            Ok(())
        } else {
            Err(PykepError::NumericalOverflow {
                operation: Self::NAME,
            })
        }
    }
}
}

Inside pykep-core, reuse the crate-private finite-input and finite-output helpers instead of duplicating them. External crates should construct the public PykepError variants as above.

The implementation rules are:

• rhs writes every derivative component into caller-owned storage.
• rhs must not allocate.
• Validate before divisions, roots, logarithms, normalizations, or other domain-sensitive operations.
• Do not silently clamp, normalize, or cross a physical singularity.
• Use InvalidInput for finite values outside the declared domain, NonFiniteInput for NaN/infinity, SingularGeometry for mathematically undefined geometry, and NumericalOverflow for a non-finite result from otherwise finite inputs.
• Keep model-domain errors distinct from IntegrationFailure.
• Return a stable, descriptive NAME; integration errors include it.
• Do not keep mutable propagation state in the model object. Constant model configuration belongs in the parameter array or in an immutable model value.

An external model is immediately usable with DOP853:

#![allow(unused)]
fn main() {
let result = Dop853.propagate(
    &OscillatorDynamics,
    InitialValueProblem::new(0.0, [1.0, 0.0], 10.0, [2.0]),
    IntegratorOptions::default(),
)?;
}

Taylor event location is not implemented. Event-driven propagation must use Dop853, even for a built-in Taylor-capable model.

3. Add Jacobians and sensitivity support

Implement DifferentiableDynamicsModel<N, P> when state-transition matrices, parameter sensitivities, or gradients are part of the model’s use case:

#![allow(unused)]
fn main() {
use pykep_core::integration::DifferentiableDynamicsModel;

impl DifferentiableDynamicsModel<2, 1> for OscillatorDynamics {
    fn jacobians(
        &self,
        time: f64,
        state: &[f64; 2],
        parameters: &[f64; 1],
        state_jacobian: &mut [[f64; 2]; 2],
        parameter_jacobian: &mut [[f64; 1]; 2],
    ) -> Result<()> {
        self.validate(time, state, parameters)?;
        let omega = parameters[0];
        *state_jacobian = [[0.0, 1.0], [-omega * omega, 0.0]];
        *parameter_jacobian = [[0.0], [-2.0 * omega * state[0]]];
        Ok(())
    }
}
}

Matrix layout is output-by-input:

state_jacobian[i][j]     = d f_i / d state_j
parameter_jacobian[i][k] = d f_i / d parameter_k

Always overwrite both complete output matrices. Never assume the caller provided zeros.

Prefer analytic Jacobians. If a family uses numerical Jacobians, reuse the existing scale-aware helper and list every strictly positive state or parameter index. Numerical differentiation must not step across a mass, radius, barrier, or other one-sided domain.

DOP853 integrates these Jacobians directly. The current Taylor sensitivity path uses centered differences of complete Taylor propagations and scales as 2W + 1 propagations for W seed directions. Keep DOP853 as the default for wide sensitivity matrices unless measurement supports a different choice.

4. Add convenience APIs

For a public built-in model, mirror the shape of the existing dynamics APIs:

• evaluate;
• propagate;
• propagate_with_method;
• propagate_with_stm or a general sensitivity method;
• propagate_with_stm_method;
• model-specific invariants or controls, if they are already part of the underlying physics.

The no-suffix propagate method may use AdaptiveIntegrator::default() only after Taylor support is complete. The method-selecting variant should accept IntegrationMethod, and documentation must say which method is the default.

Avoid inventing derived physics merely to make an API appear symmetrical. Expose only quantities actually defined by the implemented system.

5. Add the fixed-system Taylor evaluator

This section applies only to models shipped inside pykep-core.

Add the model to crates/pykep-core/src/integration/taylor/systems.rs by implementing the private coefficient contract and the public sealed marker:

#![allow(unused)]
fn main() {
impl TaylorCoefficientModel<2, 1> for OscillatorDynamics {
    fn coefficients(
        &self,
        time: f64,
        state: &[f64; 2],
        parameters: &[f64; 1],
        order: usize,
        jet: &mut [[f64; MAX_ORDER + 1]; 2],
    ) -> Result<()> {
        self.validate(time, state, parameters)?;
        oscillator_tape().coefficients(time, state, parameters, order, jet);
        Ok(())
    }
}

impl TaylorDynamicsModel<2, 1> for OscillatorDynamics {}
}

There are two implementation strategies.

Specialized recurrence

Use a hand-written recurrence for a small equation where doing so clearly reduces work. The Kepler and ZOH Kepler implementations are the templates.

The recurrence must:

1. Clear jet.
2. Copy the initial state into coefficient zero.
3. Advance only coefficient n + 1 from already available coefficients 0..=n.
4. Divide the right-hand-side coefficient by n + 1.
5. Stop at the requested order, never beyond MAX_ORDER.

Provide an independent coefficient reference in tests. Hand-written recurrences without coefficient-level tests are not acceptable.

Cached expression tape

Use TapeBuilder for a larger fixed expression:

#![allow(unused)]
fn main() {
use std::sync::OnceLock;

fn oscillator_tape() -> &'static IncrementalTape<2> {
    static TAPE: OnceLock<IncrementalTape<2>> = OnceLock::new();
    TAPE.get_or_init(build_oscillator_tape)
}

fn build_oscillator_tape() -> IncrementalTape<2> {
    let builder = TapeBuilder::new();
    let outputs = {
        let position = builder.state(0);
        let velocity = builder.state(1);
        let omega = builder.parameter(0);
        [velocity, -omega * omega * position].map(Expression::index)
    };
    builder.finish(outputs)
}
}

The lexical block is important: finish consumes the builder after all borrowed expressions have been reduced to output indices.

The tape currently supports:

• constants, time, state components, and parameters;
• addition, subtraction, multiplication, division, and negation;
• real powers and square roots;
• exponentials;
• sine, cosine, and paired sin_cos;
• stop_gradient;
• reverse symbolic gradients.

Use builder.time() for explicit time dependence. Treating the initial time as a constant produces incorrect higher-order coefficients.

TapeBuilder performs structural common-subexpression elimination. Reuse the same algebraic form where practical, but do not obscure the equations solely to reduce the operation count.

If the system needs a new elementary operation, adding only an evaluated operation is insufficient. Extend all of these together:

1. the Operation enum;
2. the builder and Expression API;
3. constant folding;
4. the incremental Taylor recurrence;
5. reverse differentiation, when gradients can traverse the node;
6. companion/workspace metadata, if required;
7. direct elementary-function and derivative tests in tape.rs.

Piecewise expressions and optimal control

A Taylor series is valid only on an analytic branch. For algebraically equivalent stable formulas, keep one OnceLock tape per branch and select the branch from coefficient-zero state/parameter values at the start of each step. Test every branch and its boundary policy.

For minimized Pontryagin Hamiltonians, controls may evolve as Taylor series while remaining excluded from the Hamiltonian partial derivative. Wrap those direction or throttle expressions with stop_gradient, matching the envelope convention. Omitting this changes the costate equations.

Do not smooth a real discontinuity merely to make Taylor integration possible. Document the boundary and use segmented propagation or DOP853 when the physical model requires it.

6. Add ZOH support when applicable

If controls are constant over user-supplied segments, implement the internal ZeroOrderHoldModel<N, C, K, P> mapping:

parameters = parameters(control[C], constants[K])

Also map control and constant sensitivity seeds into parameter seeds. Define which side owns an exact switch time, validate the complete schedule once, and test forward and backward segment traversal.

Add the model to every relevant ZOH dispatch point:

• scalar segment propagation;
• dense/history output;
• sensitivities;
• leg support, if the state/control semantics match the leg;
• Python ZohModel selection, if exposed there.

Do not hide a control lookup inside rhs; a validated schedule should split the integration into constant-parameter segments.

7. Add ordered batch APIs

Provide a batch when callers will naturally evaluate many independent states, controls, epochs, or final times. A batch is an explicit API extension, not an implicit change to scalar behavior.

Use pykep_core::batch::try_map so that:

• workers = 0 uses Rayon’s global pool;
• workers = 1 is deterministic serial execution;
• larger values use a cached fixed-size pool;
• outputs preserve input order;
• if multiple items fail, the earliest failing input is reported.

Test empty input, one item, multiple worker counts, mismatched input lengths, deterministic output order, and deterministic error order. Every batch result must match the scalar API for the same item.

Avoid nested parallelism. If the objective is parallelized at a higher level, call the scalar or workers = 1 model API inside each objective evaluation.

8. Add Python bindings

Python exposure is optional for a low-level internal model and expected for a user-facing pykep model.

Implement bindings in crates/pykep-py/src/dynamics.rs:

1. Parse dynamically sized Python inputs into fixed Rust arrays before releasing the GIL.
2. Call the same pykep-core scalar implementation; do not duplicate the equations in the binding.
3. Convert PykepError with the shared to_python mapping.
4. Release the GIL with Python::detach for propagation and batch work.
5. Use NumPy N × state_dimension arrays for batch states.
6. Accept workers on parallel batches and preserve input order.
7. Register every function in dynamics::register.

Then update:

• python/pykep_rust/__init__.py imports and __all__;
• python/pykep_rust/_pykep_rust.pyi signatures and docstrings;
• python/tests/test_smoke.py for scalar values and error mapping;
• python/tests/test_parallel_batch.py for scalar/batch parity, shapes, worker counts, and deterministic failures;
• a Python example when the usage is not obvious.

Use consistent names:

<system>_rhs
<system>_rhs_batch
propagate_<system>
propagate_<system>_batch
propagate_<system>_with_stm
propagate_<system>_with_stm_batch

If the Rust model supports both DOP853 and Taylor but Python intentionally exposes only the default, state that explicitly. Do not expose internal tape objects.

9. Build the test pyramid

No single reference is sufficient. Add the following layers in order.

Evaluated right-hand side

• One or more hand-computed or authoritative reference states.
• Exact zero/control-free/equilibrium reductions where available.
• Explicit-time cases at more than one time.
• Forward-frame and sign-convention checks.
• Every invalid parameter range.
• Every singular geometry.
• NaN and infinity for time, state, and parameters.
• Finite inputs that would overflow the result.

Put small local tests beside the implementation. Put source-parity and cross-family tests in a dedicated crates/pykep-core/tests/phase*_*.rs integration test.

Jacobians

Compare every state and parameter column with independent, scale-adjusted central differences over several non-singular states. Use perturbations based on each variable’s scale, not one absolute epsilon for all columns.

Also test:

• the documented output-by-input layout;
• analytically zero columns;
• one-sided positive domains;
• consistency between the Jacobian and propagated sensitivities.

Propagation

Use at least two independent checks:

• a closed-form solution, invariant, upstream trajectory, or independently generated fixture;
• a same-problem DOP853 comparison at a tighter reference tolerance.

Cover forward and backward time, zero duration, dense sampling, step limits, rejected steps when relevant, and a duration long enough to expose drift or branch errors. Validation tolerances must be justified by the state scale and reference accuracy; they are not general user guarantees.

Taylor coefficients

Keep a test-only full-series form of the right-hand side and compare its jet with the optimized recurrence/tape at representative orders such as [8, 15, MAX_ORDER].

The test should verify:

• every state component and every coefficient;
• explicit-time coefficients;
• trigonometric/exponential/power paths used by the model;
• every stable algebraic branch;
• an operation-count ceiling for expression tapes.

Choose a scaled relative comparison:

abs(incremental - reference) <= tolerance * max(abs(reference), 1)

High-order coefficients of poorly scaled systems can be ill-conditioned. Use a realistic state, record why a looser coefficient tolerance is needed, and retain the independent propagated-state comparison.

Cross-backend and regression tests

Add the model to the central Taylor-versus-DOP853 tests. If an upstream heyoka/C++ implementation exists, generate committed fixtures with a pinned upstream version and keep the generator outside the normal runtime dependency graph.

Compare physical outputs, not only internal work counters. Taylor coefficient sweeps and DOP853 RHS evaluations are not equivalent units.

Batch and Python tests

For every scalar Python function, test the corresponding batch when one exists. Verify dtype/shape, empty batches, mismatched lengths, exception classes, worker behavior, and equality to scalar calls.

10. Add benchmarks and apply a performance gate

Add Criterion entries in crates/pykep-core/benches/dynamics.rs for:

• one RHS evaluation;
• one representative DOP853 propagation;
• one Taylor propagation for a Taylor-capable built-in;
• one state/parameter sensitivity propagation when it is important;
• representative scalar and batch throughput.

Benchmark release builds. Warm lazy OnceLock tapes before recording steady state, but report first-call setup separately if it is material.

For an optimization, declare the benchmark, correctness tolerance, and minimum meaningful gain before changing the implementation. Pin:

• initial time, state, parameters, and final time;
• integration tolerances and maximum step;
• release profile and dependency versions;
• iteration/sample counts;
• CPU affinity when comparing small kernels;
• final-state checksum or independent accuracy measure.

Keep an optimization only when the gain is larger than run-to-run noise and all accuracy checks remain satisfied. Record development-host timings as evidence, not portable latency promises.

11. Document the model

Every public type and method needs Rustdoc describing:

• state and parameter order;
• units and reference frame;
• default integration backend;
• singularities and domain restrictions;
• output layout;
• # Errors;
• a short runnable example for a non-obvious API.

Update the relevant mdBook pages:

• docs/dynamics.md, docs/zero-order-hold.md, or docs/pontryagin.md;
• docs/taylor-integration.md when Taylor support changes;
• docs/python-api.md and docs/python-migration.md for Python;
• docs/batch-processing.md for a new batch family;
• docs/validation.md with the independent evidence and tolerances;
• docs/performance.md with the benchmark protocol and interpretation;
• docs/status.md and docs/source-map.md;
• docs/SUMMARY.md for any new page;
• CHANGELOG.md under Unreleased.

Update model counts such as “eleven built-in models” wherever the new model changes them. Search for the old count rather than fixing only one page:

rg -n "eleven|11 built-in|supported models|TaylorDynamicsModel" \
  README.md CHANGELOG.md docs crates

If the model ports an upstream equation, add the exact source file, upstream version/commit, and Rust destination to docs/source-map.md.

12. Run the complete quality gate

From the public repository:

cargo fmt --all -- --check
cargo clippy --workspace --all-targets --all-features -- -D warnings
cargo test --workspace --all-features --locked
cargo test --workspace --all-features --locked --release
cargo test -p pykep-core --no-default-features --locked
RUSTDOCFLAGS="-D warnings" cargo doc --workspace --all-features --no-deps
cargo +1.88.0 check --workspace --locked
python tools/check_markdown_links.py
mdbook build

For Python changes, build the extension in an isolated environment and run the Python suite:

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

Also run the relevant Criterion benchmark and coverage report:

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

Use Miri for new unsafe-sensitive structure or recurrence code if applicable. The crate forbids unsafe code, but Miri can still catch invalid assumptions in dependency-free state manipulation. Fuzz new parsers or dynamically shaped input boundaries; do not fuzz a fixed smooth RHS merely to increase a metric.

Definition of done

A built-in ODE system is complete only when:

• State, parameter, unit, frame, and domain contracts are written down.
• DynamicsModel is allocation-free, validated, and fully documented.
• Errors use the stable public taxonomy.
• Jacobians are implemented and independently checked, or their absence is explicitly justified.
• Taylor support is implemented and coefficient-tested, or the model is explicitly documented as DOP853-only.
• Propagation matches at least two independent references or invariants.
• Forward, backward, zero-duration, singular, and non-finite cases pass.
• Every analytic branch and control switch policy is tested.
• Scalar and batch APIs agree for all worker modes.
• Python functions, stubs, exports, error mapping, and tests are complete when Python exposure is in scope.
• RHS, propagation, sensitivity, and batch benchmarks are recorded where relevant.
• Rustdoc, mdBook, status, source map, validation, performance, and changelog are updated.
• Formatting, Clippy, debug/release/no-default tests, Rustdoc, MSRV, Markdown links, mdBook, Python tests, and relevant coverage checks pass.
• The public worktree contains only intended changes and the work is committed.

ADR 0001: Fixed numerical types

• Status: accepted
• Date: 2026-07-25

Context

The common pykep values are three-vectors, six-element Cartesian states or orbital elements, and small fixed matrices. Dynamic allocation and a general-purpose public matrix type would add cost and couple the API to a dependency without improving these fixed-shape contracts.

Decision

Parity work starts with documented aliases over stack-allocated Rust arrays:

• Vector3 = [f64; 3];
• CartesianState = [f64; 6];
• Elements6 = [f64; 6];
• Matrix3 = [[f64; 3]; 3];
• Matrix6 = [[f64; 6]; 6].

Small matrix helpers use const generics internally and in deliberately generic public operations. Variable-sized solution families, controls, grids, and batches use owned vectors. No external matrix type appears in the public API.

Consequences

The common hot path is allocation-free and Python conversion remains straightforward. Aliases do not prevent semantic mix-ups, so public functions must use meaningful parameter names and named result structs where multiple arrays would be ambiguous. Newtypes remain an option before API stabilization if later phases demonstrate enough safety benefit.

Phase 4 exercised that option for semantically distinct element sets: ClassicalElements and ModifiedEquinoctialElements provide named fields and lossless [f64; 6] conversions. Cartesian states retain the fixed alias because their ordering is unambiguous in every consuming API.

ADR 0002: Numerical error taxonomy

• Status: accepted
• Date: 2026-07-25

Context

The C++ implementation mixes exceptions, NaN sentinels, and unchecked floating-point behavior. A mission-analysis library must not turn failed validation or convergence into a plausible output.

Decision

All fallible public numerical operations return pykep_core::Result<T> with a non-exhaustive PykepError. Stable categories distinguish:

• invalid values and non-finite values;
• singular geometry;
• convergence failure;
• dimension mismatch;
• unsupported capabilities;
• floating-point overflow;
• integration failure.

Input validation happens at public boundaries. Private, already-validated helpers may avoid repeated checks in iterative hot paths. Panics are reserved for violated internal invariants, not caller input.

The Python layer maps invalid values and dimensions to ValueError, numerical overflow to OverflowError, and exposes package exceptions for convergence, singular geometry, unsupported capabilities, and integration failures.

Consequences

This intentionally differs from upstream functions that return NaN for an invalid anomaly domain or accidentally return finite limit values for NaN inputs. Error strings retain useful context, but callers should match the category rather than punctuation.

ADR 0003: embedded thresholded VSOP2013 evaluator

• Status: accepted
• Date: 2026-07-25

Context

The upstream pykep provider asks heyoka to construct VSOP2013 expressions and JIT-compile a function for each (body, threshold) pair. That is not suitable for a C/C++-free Rust crate, reproducible offline builds, or bounded global state.

The authoritative IMCCE solution contains 2,607,947 terms. In heyoka’s packed layout, retaining every term would require about 100 MiB before compression. The pykep tests evaluate at thresholds down to 1e-9; 112,270 terms survive that floor and occupy 4.3 MiB in the committed encoding. The ordinary pykep default is 1e-5, for which only 1,858 terms survive across all bodies.

The coefficient source in heyoka 7.10.0 is MPL-2.0-covered code generated from the IMCCE VSOP2013 release. The public data notice pins the exact heyoka commit, generator, retained threshold, and binary hash.

Decision

The default vsop2013 Cargo feature embeds the original integer multipliers and binary64 coefficients down to 1e-9. A provider validates and decodes only its selected planet, applies the requested threshold once, and owns the result through Arc. Evaluation is a direct Rust series evaluator followed by the published equinoctial-to-Cartesian transformation and ICRF rotation.

Thresholds below 1e-9 return an explicit validation error. Disabling default features removes the coefficient asset; construction then returns an explicit unsupported-capability error. Availability and the threshold floor are queryable in Rust and Python.

There is no process-global provider cache. Default-threshold initialization measured about 11.6 µs, so a cache would add synchronization and lifetime complexity without justifying itself. Cloning an existing provider shares its immutable decoded data.

No source checkout, generator, network access, C++, LLVM, or heyoka is needed to build or run the public crate.

Measurements

An orientation run on an AMD Ryzen 9 9950X measured:

The Rust and C++ timing harnesses use the same Earth–Moon provider and epoch. These are orientation results without fixed CPU frequency, not release performance guarantees.

Across 54 golden states at 1e-9, the largest absolute difference from the JIT-compiled upstream path was 0.185 m in position and 6.6e-8 m/s in velocity. The difference is consistent with floating-point summation and trigonometric evaluation order.

Consequences

• Startup, builds, and runtime dependencies are bounded and reproducible.
• Keplerian and JPL low-precision providers remain available without the feature.
• The full sub-1e-9 theory is intentionally not embedded because its data and evaluation cost are disproportionate to the upstream-tested public use.
• The direct evaluator is slower than optimized JIT code at 1e-9; Phase 18 may profile argument precomputation, SIMD-friendly term grouping, and batch evaluation while golden tests remain active.

ADR 0004: adaptive integration backend

• Status: accepted
• Date: 2026-07-25

Context

The upstream ta modules return heyoka expression graphs and cached Taylor-adaptive integrators. Rust cannot preserve that implementation API without a C++/LLVM runtime. The remaining physical models and legs instead need an evaluated-model interface with:

• six-, seven-, fourteen-, and augmented state dimensions up to 126;
• constant parameters that can change exactly between ZOH segments;
• scalar relative and absolute tolerances down to the upstream test regimes;
• initial/maximum step and accepted/rejected-step limits;
• forward and backward propagation;
• deterministic exact-final-time stops;
• dense evaluation and terminal zero-crossing events;
• state and parameter Jacobians plus arbitrary first-order sensitivity seeds;
• no heap allocation per accepted step on nominal fixed-size hot paths.

Kepler and CR3BP were used as representative smooth and close-approach systems. ZOH switching is treated as a known discontinuity: a leg ends one solve exactly at a switch and starts the next with new parameters. It does not ask a continuous interpolant to cross a discontinuity.

Candidates

The research snapshot was taken on 2026-07-25.

differential-equations 0.6.1 provides pure-Rust DOP853, seventh-order dense interpolation, backward solves, terminal root finding, scalar/component tolerances, and explicit step and rejection limits. It accepts stack-allocated arrays and user-defined state types. Its Apache-2.0 license is compatible, and a no-default-feature build passes Rust 1.88 even though the crate does not declare an MSRV.

Its simba dependency currently brings paste 1.0.15. RustSec advisory RUSTSEC-2024-0436 marks that macro crate unmaintained but reports no vulnerability and no safe upgrade. The exact advisory is narrowly allowlisted in deny.toml and remains visible in cargo audit output.

ode_solvers 0.6.2 provides pure-Rust DOP853 and dense output under Apache-2.0. The same-input spike was faster, but the implementation allocates vectors and intermediate states inside each solve/step. It has a post-step stop callback rather than interpolated root location and no first-order sensitivity interface.

diffsol 0.16.1 has the richest event, reset, dense-output, forward, and adjoint sensitivity system. It was screened out before the timing final because its implicit-solver and matrix abstractions are disproportionate for these small non-stiff explicit systems. Its optional DSL/JIT path also solves a problem this port explicitly avoids.

A local DOP853 implementation would give maximum workspace control but would make pykep responsible for a second adaptive-solver implementation, dense coefficients, root location, step controller, and long-term maintenance. C++ FFI was excluded by the release requirements.

Decision

Use differential-equations 0.6.1 DOP853 behind the pykep-owned integration facade. No backend type appears in a public signature.

DynamicsModel<N, P> evaluates into caller-owned arrays and validates the physical domain. DifferentiableDynamicsModel<N, P> adds row-major state and parameter Jacobians. InitialValueProblem owns time, state, and constant parameters. SensitivityProblem supplies arbitrary dstate/dseed and dparameters/dseed matrices, integrating

dS/dt = (df/dstate) S + (df/dparameters) (dparameters/dseed).

This directly represents STMs, parameter variations, and the state/control columns needed by ZOH legs. The facade provides final-state, dense-sampling, terminal-event, and sensitivity propagations with stable pykep errors and work counters.

The ordinary and variational final-state paths use fixed-size state storage and an output callback that does not retain internal steps. Dense and event APIs allocate returned samples intentionally; the event adapter currently retains accepted steps because the selected crate’s event wrapper owns its output strategy. Events are not on the remaining leg hot path.

Measurements

All values are orientation measurements on the same AMD Ryzen 9 9950X host. CPU frequency was not pinned.

The Git-tracked candidate tool used 30 samples of 200 complete Kepler propagations with rtol = atol = 1e-12:

The maintained workspace Criterion profile measured the selected nominal path at 11.865 µs (95% estimate 11.849–11.884 µs) and the state plus 6 × 6 STM path at 85.663 µs (85.260–86.102 µs).

The matching C++/heyoka harness measured 3.722 µs mean nominal steady-state propagation and 84.440 µs mean for the STM. Cloning the cached C++ nominal integrator once cost 0.467 ms; the steady-state result excludes that clone.

Validation and risks

• A 100-orbit eccentric Kepler solve kept absolute specific-energy drift below 2e-10 at 1e-13 tolerances.
• Kepler final states and the STM match the independent analytic propagator; parameter columns match central finite differences.
• A CR3BP close approach preserves its Jacobi constant and reverses within scale-aware tolerances.
• Rejected steps, exhausted step limits, backward integration, dense samples, terminal roots, non-finite values, and bit-for-bit repeatability are tested.
• DOP853 is not a Taylor method. Equal numeric tolerance values do not imply equal internal steps or bitwise parity with heyoka.
• Nominal integration is about three times slower than warmed C++/heyoka in this case. The variational path is approximately equal. Phase 18 will profile before changing the solver.
• Dense interpolation error depends on internal step length as well as local solve tolerance. High-accuracy dense/event callers should set maximum_step; ZOH switches always use exact final-time solves.
• Close singularities may exhaust rejections or step size. Errors retain the model and integration context instead of returning a plausible state.

Consequences

The backend supports every remaining ta and generic ZOH-leg requirement without C++, LLVM, a JIT, global mutable caches, or a hidden thread pool. Physical Kepler, CR3BP, BCP, ZOH, and Pontryagin models remain Phase 11–13 work; completing this gate does not claim those APIs are implemented.

The dependency is deliberately hidden so it can be upgraded or replaced without changing model, problem, result, or Python contracts. The transitive paste maintenance advisory must be reconsidered on every backend upgrade; the allowlist must not be broadened to a vulnerability.

ADR 0005: Add Taylor as the built-in nominal backend

• Status: accepted
• Date: 2026-07-28
• Supersedes: ADR 0004’s default for built-in nominal dynamics; DOP853 remains the general backend

Context

ADR 0004 selected DOP853 because it supplied final state, dense output, events, and direct variational equations behind a small pure-Rust dependency. The same evidence showed a nominal high-accuracy performance gap relative to the upstream warmed heyoka implementation.

The eleven pykep dynamics types come from only nine fixed equation families. They do not require arbitrary symbolic expressions, LLVM, events in Taylor mode, heyoka batch mode, or arbitrary precision. A disposable eccentric Kepler prototype passed a predeclared matched-error gate at tight tolerances.

Decision

Add a private fixed-capacity Taylor-series engine inside pykep-core. Implement coefficient kernels only for the built-in models. Publish Taylor, IntegrationMethod, and AdaptiveIntegrator, but keep the coefficient trait private until a real third-party system requires it.

Make Taylor the default of AdaptiveIntegrator and of the nominal Kepler/CR3BP/BCP convenience methods. This follows upstream pykep, where these TA-derived equation families are exposed as heyoka Taylor-adaptive integrators. Keep direct Dop853 calls, generic user dynamics, event location, and the existing direct-variational sensitivity convenience methods on DOP853. Add explicit *_with_method variants for built-in dynamics and ZOH schedules.

Use tolerance-selected order 8–24 and component-scaled step selection from the final two coefficients. Preserve an optimized Kepler recurrence where full-series generic evaluation would erase the measured crossover.

For the first release, dense grids clip Taylor steps at sample times and seeded sensitivities use centered complete propagations. Document their cost instead of claiming the not-yet-implemented polynomial interpolation and direct Taylor variational optimizations.

Consequences

• No C/C++ runtime, LLVM, JIT latency, or new runtime dependency is added.
• At 1e-14 and machine epsilon the representative Kepler solve is faster and more accurate than DOP853; at 1e-9 it is slower.
• Built-in nominal dynamics default to the algorithm family used by upstream pykep; callers can still request DOP853 explicitly.
• User-defined DynamicsModel implementations continue to use DOP853.
• Wide STMs and event-driven problems should continue to use DOP853.
• Equation duplication between evaluated and coefficient kernels is accepted as localized technical debt, protected by upstream, DOP853, and official-heyoka cross-validation.

The raw measurements and validation boundaries are recorded in High-accuracy Taylor integration.