empyrean

Safe Rust wrapper over libempyrean — uncertainty-first orbit propagation, ephemeris, orbit determination, and event detection for asteroids and comets, powered by automatic differentiation

The idiomatic Rust API over the libempyrean C ABI. Every C function exposed in the cdylib has a typed, Result<_, Error>-returning wrapper here. RAII handles the underlying allocations so callers never juggle raw FFI pointers.

[dependencies]
empyrean = "0.9.0"

What it does

Propagation — N-body (Sun, planets, Moon, Pluto) with EIH general relativity, Sun J2 and Earth J2–J4 zonal harmonics, 16 asteroid perturbers, and the Marsden non-gravitational model — selectable across Approximate / Basic / Standard force-model tiers (Standard is the default). GR15 and DOP853 integrators. Optional finite-burn thrust arcs — constant-RTN, velocity-tangent, or inertial-fixed steering, with per-arc Δv targeting corrections — layer on as a continuous-thrust force input.
Uncertainty — First-order (Jet1) state transition matrices; second-order (Jet2) state transition tensors; unscented sigma-point and Monte Carlo sampling; an adaptive Auto mode that escalates the method automatically through close approaches and relaxes it elsewhere. Optional per-epoch tagged-covariance readback.
Ephemeris — RA/Dec, rates, photometry (H–G, H–G₁G₂, H–G₁₂), light time, phase angle, solar elongation, local horizon. Each row carries the 6×6 sky-plane covariance over (ρ, RA, Dec) and their rates, and the aberrated barycentric ICRF state at the photon-emission epoch with its own 6×6 covariance — both present when the input orbit carries a state covariance.
Orbit determination — Gauss, Herget, and systematic-ranging (admissible region + Manifold of Variations) IOD → N-body differential correction over optical and radar (delay / Doppler) observations, with STM caching and outlier rejection. Solves beyond the six-element state for the Marsden A1/A2/A3 non-gravitational block, the cometary outgassing time delay DT, the SRP area-to-mass ratio AMRAT, and thrust Δv-correction segments — each partial supplied analytically by the hyperdual integrator — and returns a tagged solved covariance that names every fitted parameter, plus an event-aware trust verdict on the delivered covariance. Optional post-fit photometry recovers H and the phase-function slope. Validated against find_orb and JPL SBDB.
Events — Close approach (start/end), periapsis, gravitational capture (start/end), shadow entry/exit, atmospheric entry/exit, impact, and possible impact.

Quick start

use empyrean::{Context, Epoch, PropagationConfig};

let ctx = Context::from_data_dir(None)?;

// Query SBDB for Apophis and propagate through its 2029 Earth flyby.
let orbits = empyrean::query_sbdb(&["Apophis"], None)?.orbits;
let epochs = vec![Epoch::from_mjd_tdb(65000.0)];
let result = ctx.propagate(&orbits, &epochs, &PropagationConfig::default())?;

println!("{} states, {} events", result.states.len(), result.events.len());
# Ok::<(), empyrean::Error>(())

Orbit determination

determine runs a full IOD (Gauss / Herget / systematic ranging) → N-body differential correction; refine is a Bayesian update against a prior orbit; evaluate returns residuals without fitting. The fitted result.orbit is a re-feedable [Orbit] carrying state, covariance, and any fitted non-gravitational parameters — pass it straight back into propagate, generate_ephemeris, or compute_impact_probabilities.

# use empyrean::{Context, ODConfig};
# let ctx = Context::from_data_dir(None)?;
let obs = ctx.read_ades("observations.psv")?;   // optical + radar
let result = ctx.determine(&obs, None, &ODConfig::default())?;

println!(
    "converged={}, RMS = {:.2}\" RA / {:.2}\" Dec",
    result.converged,
    result.summary.rms_ra_arcsec,
    result.summary.rms_dec_arcsec,
);
# Ok::<(), empyrean::Error>(())

Every residual row carries per-observation diagnostics: χ² with its survival probability, along/cross-track residuals with the full symmetric 2×2 covariance, and influence measures including the D-optimality information loss on removal (+∞ marks an observation whose removal makes the normal matrix singular). Radar rows carry a typed delay / Doppler block — observed − predicted in seconds / hertz, with χ², survival probability, and the combined observed+predicted variance. result.covariance_trust is an event-aware verdict on the delivered covariance: Trusted, EncounterIntervenes (naming the intervening close approach or high-nonlinearity crossing, and whether a second-order state-only correction can recover it), or WeaklyDeterminedHighN for wider-than-state fits. None means no trust gate ran — absence of a verdict is not trust.

Wide-parameter fitting

Beyond the six-element state, determine and refine can solve for the Marsden A1/A2/A3 non-gravitational block, the cometary outgassing time delay DT, the SRP area-to-mass ratio AMRAT, and thrust Δv-correction segments — every partial derivative supplied analytically by the hyperdual integrator rather than finite differences. Choose the axes with SolveForParams: StateOnly, StateAndNonGrav, Auto (starts state-only and escalates the non-grav block automatically on a poor fit), or Explicit(SolveFor { .. }) for the wider axes the coarse variants can't name.

DT, AMRAT, and thrust are refine-path solves: the input orbit must carry a prior — the variance that opens the parameter. Request an axis without its prior and the fit errors loudly; it never hands back a zeroed or defaulted column.

Every wide fit reports a SolvedCovariance whose fitted-parameter identities travel with the matrix. Read a parameter's variance by its slot (marsden_slot, dt_slot, amrat_slot, thrust_slots) rather than by guessing column order — width alone is ambiguous (a 9×9 is Marsden-only or one thrust segment).

use empyrean::{Context, ODConfig, SolveFor, SolveForParams};

let ctx = Context::from_data_dir(None)?;
let obs = ctx.read_ades("comet_67p.psv")?;

// First solve state + Marsden A1/A2/A3.
let fit = ctx.determine(&obs, None, &ODConfig {
    solve_for: SolveForParams::StateAndNonGrav,
    ..Default::default()
})?;

// Refine, additionally solving the outgassing time delay DT. Opening DT
// requires a prior on it — its variance (days²) — carried on the orbit.
// Ask for DT without the prior and refine errors, never a zeroed column.
let prior = fit.orbit
    .with_non_grav_dt(Some(30.0))
    .with_non_grav_dt_variance(Some(100.0));
let refined = ctx.refine(&prior, &obs, &ODConfig {
    solve_for: SolveForParams::Explicit(SolveFor {
        marsden: true,
        dt: true,
        ..Default::default()
    }),
    ..Default::default()
})?;

// The solved covariance names its columns — read σ(DT) by slot.
if let Some(cov) = &refined.solved_covariance {
    if let Some(k) = cov.dt_slot {
        println!(
            "ΔDT = {:?} d,  σ(DT) = {:.3} d",
            refined.dt_delta,
            cov.matrix[k][k].sqrt(),
        );
    }
}
# Ok::<(), empyrean::Error>(())

Post-fit photometry

Attach a PhotometryConfig to ODConfig::photometry and the pipeline recovers absolute magnitude H and the phase-function slope from the observation magnitudes after the orbit is solved. The photometric fit has no astrometric partials, so it never touches the state. In Auto it climbs a model ladder — H-only → HG12 → HG1G2 (Muinonen et al. 2010) — admitting the richest model the arc's phase-angle coverage supports and reporting the one it actually fit on model_used (never Auto). H carries an honest 1σ through the fitted covariance; the per-model gate decisions come back in gates. Magnitudes whose band has no adopted V-band conversion are excluded and counted — n_mags_dropped_unconvertible, with the distinct offending band codes in dropped_bands — and the observations' astrometry is unaffected.

use empyrean::{Context, ODConfig, PhotometryConfig};

let ctx = Context::from_data_dir(None)?;
let obs = ctx.read_ades("observations.psv")?;

// Fit the orbit, then fit H/G from the magnitudes (Auto ladder:
// H-only -> HG12 -> HG1G2).
let fit = ctx.determine(&obs, None, &ODConfig {
    photometry: Some(PhotometryConfig::default()),
    ..Default::default()
})?;

if let Some(phot) = &fit.photometry {
    let sigma_h = phot.covariance.map(|c| c[0][0].sqrt());
    println!(
        "H = {:.2} ± {:.2}  ({:?}, {} mags, α span {:.1}°)",
        phot.h,
        sigma_h.unwrap_or(f64::NAN),
        phot.model_used,
        phot.n_mags_used,
        phot.alpha_span_deg,
    );
}
# Ok::<(), empyrean::Error>(())

Ephemeris

# use empyrean::{Context, EphemerisConfig, Epoch};
# let ctx = Context::from_data_dir(None)?;
# let orbits = empyrean::query_sbdb(&["Apophis"], None)?.orbits;
let epochs = vec![Epoch::from_mjd_tdb(65000.0)];
let observers = ctx.get_observers(&["W84", "F51"], &epochs)?;
let eph = ctx.generate_ephemeris(&orbits, &observers, &EphemerisConfig::default())?;

for entry in &eph.entries {
    println!("RA {:.4}°  Dec {:.4}°  V {:.2}", entry.ra_deg, entry.dec_deg, entry.mag);
}
# Ok::<(), empyrean::Error>(())

Beyond the printed astrometry, each EphemerisEntry carries the 6×6 sky-plane covariance over (ρ, RA, Dec) and their rates (AU / degree units), and the aberrated — light-time corrected — barycentric ICRF Cartesian state at the photon-emission epoch with its own 6×6 covariance; both covariances are None when the input orbit carried no state covariance. Non-fatal generation warnings (an Earth-orientation kernel coverage gap handled by the analytic IAU 2006 fallback, a row whose observation-sensitivity chain was skipped) come back on EphemerisResult::warnings — empty when the run had nothing to report.

Uncertainty

First-order (the default) propagates the covariance with the state-transition matrix — accurate when the orbit is approximately linear over the uncertainty region. Second-order adds the state-transition tensor for the curvature that linear covariance misses near a close approach.

# use empyrean::{Context, Epoch, PropagationConfig, UncertaintyMethod};
# let ctx = Context::from_data_dir(None)?;
# let orbits = empyrean::query_sbdb(&["Apophis"], None)?.orbits;
# let epochs = vec![Epoch::from_mjd_tdb(65000.0)];
let config = PropagationConfig {
    uncertainty_method: UncertaintyMethod::SecondOrder,
    ..Default::default()
};
let result = ctx.propagate(&orbits, &epochs, &config)?;
# Ok::<(), empyrean::Error>(())

Continuous thrust

Model finite burns / low-thrust arcs by attaching a ThrustParams to an orbit before propagation. Each ThrustArc carries its own thrust, mass, specific impulse, steering law (constant-RTN, velocity-tangent, or inertial-fixed), and central body; the burn perturbs the trajectory through the same differentiated dynamics as gravity and the non-gravitational forces.

use empyrean::{Context, Epoch, Origin, PropagationConfig, SteeringLaw, ThrustArc, ThrustParams};

let ctx = Context::from_data_dir(None)?;
let orbit = empyrean::query_sbdb(&["Apophis"], None)?.orbits.remove(0);

// One finite burn: 1 N over MJD 65000–65010 on a 500 kg spacecraft,
// mass depleting at Isp = 3000 s, steered at constant RTN angles
// relative to the Sun. `sharpness` sets the tanh on/off transition.
let arc = ThrustArc::new(
    65000.0,                                                   // start_mjd_tdb
    65010.0,                                                   // end_mjd_tdb
    1.0,                                                       // thrust_n (N)
    500.0,                                                     // mass_kg
    100.0,                                                     // sharpness (1/day)
    SteeringLaw::ConstantRTN { alpha_rad: 0.0, beta_rad: 0.0 },
    Origin::SUN,                                               // RTN frame reference
)
.with_isp(Some(3000.0));

// Attach to the orbit and propagate. Add per-arc Δv targeting
// corrections with `ThrustParams::new(arcs).with_dv_corrections(..)`.
let orbit = orbit.with_thrust(Some(ThrustParams::new(vec![arc])));
let epochs = vec![Epoch::from_mjd_tdb(65020.0)];
let result = ctx.propagate(&[orbit], &epochs, &PropagationConfig::default())?;
println!("{} states", result.states.len());
# Ok::<(), empyrean::Error>(())

System handles

Assembling the force model (planets, Moon, asteroid perturbers, harmonics, relativistic corrections) has a fixed per-call cost. A [BuiltSystem] assembles it once for a frozen {force model, frame, encounter-timescale divisor} key and reuses it across many propagations — the build-once, propagate-many pattern for short-arc campaigns. It is Send + Sync, so &handle can be shared across threads. A call whose config disagrees with the frozen key, or that pairs the handle with a different data instance, is rejected loudly by axis — never silently rebuilt against the wrong dynamics.

# use empyrean::{Context, ForceModelTier, Frame, PropagationConfig, Epoch};
# let ctx = Context::from_data_dir(None)?;
# let orbits = empyrean::query_sbdb(&["Apophis"], None)?.orbits;
// Build once; freeze the divisor at the engine default (0.0).
let handle = ctx.built_system(ForceModelTier::Standard, Frame::EclipticJ2000, 0.0)?;

let epochs = vec![Epoch::from_mjd_tdb(65020.0)];
let result = handle.propagate(&ctx, &orbits, &epochs, &PropagationConfig::default())?;
println!("{} states", result.states.len());

// describe() reports the reproducibility record: the force-model menu
// plus the identity (SHA-256) of every loaded kernel.
let desc = handle.describe()?;
println!("{} perturbers, {} kernels", desc.perturber_origins.len(), desc.kernels.len());
# Ok::<(), empyrean::Error>(())

Impact probability and B-plane geometry

For each detected close approach you can ask for an impact-probability assessment or a full B-plane breakdown, and run several uncertainty methods side-by-side on the same encounter. Each returns one record per (method × orbit × body), tagged with its method and closest-approach epoch. Each record also carries the geodetic impact point on the body's reference ellipsoid (latitude / longitude / altitude — NaN when no surface projection is available for the encounter), the 95% binomial confidence half-width on the Monte-Carlo fraction, the second-order corrected mean miss distance with its 1σ uncertainty and skewness, the closest-approach distance gradient and 6×6 Hessian with respect to the initial state, and the adaptive Gaussian-mixture component count — fields a given method didn't compute carry NaN / 0 sentinels.

# use empyrean::{Context, Epoch, UncertaintyMethod, Origin};
# let ctx = Context::from_data_dir(None)?;
# let orbits = empyrean::query_sbdb(&["Apophis"], None)?.orbits;
let end = Epoch::from_mjd_tdb(65000.0);

let ips = ctx.compute_impact_probabilities(
    &orbits,
    end,
    &[UncertaintyMethod::FirstOrder, UncertaintyMethod::SecondOrder],
    &[Origin::EARTH, Origin::MOON],
)?;
for ip in &ips {
    println!("{:?}: miss {:.0} km", ip.body, ip.miss_distance_km);
}

let bps = ctx.compute_b_planes(&orbits, end, &[UncertaintyMethod::SecondOrder], &[Origin::EARTH])?;
for bp in &bps {
    println!("B·T {:.1} km, B·R {:.1} km", bp.b_dot_t_km, bp.b_dot_r_km);
}
# Ok::<(), empyrean::Error>(())

Runtime requirement

This crate (via empyrean-sys) loads libempyrean.{dylib,so} at run time, which is distributed separately as a binary release on GitHub and inside the published Python wheel. The path is resolved from the EMPYREAN_LIB environment variable if set, else a libempyrean.* sitting next to the loaded module, else a build-time location — an EMPYREAN_LIB_DIR override, a sibling ../target/release build, or a checksum-pinned prebuilt downloaded from the GitHub release (in that order); no system library path setup is required.

Prebuilt engine binaries are currently published for four targets — macOS arm64 (macos-aarch64), macOS x86_64 (macos-x86_64), Linux x86_64 (linux-x86_64), and Linux aarch64 (linux-aarch64); on other targets the build stops with an error unless EMPYREAN_LIB_DIR points at an engine build.

The full distribution surface (Python wheel, CLI binary, C SDK, this Rust crate) lives at the main repository — see its README for installation paths and the cross-channel quickstart.

Accuracy

Validated against JPL Horizons, ASSIST, and find_orb on 43 objects across 13 dynamical populations (NEOs, MBAs, Trojans, TNOs, comets, and more). Sub-meter propagation accuracy on bounded timescales; see the validation notes in the main repository for the comparison setup.

No guarantee of accuracy

empyrean performs numerical computations used in planetary-science and mission-planning contexts. Outputs should not be used as the sole basis for any decision — including but not limited to impact monitoring, mission planning, collision avoidance, or navigation — without independent verification. See the LICENSE file for the full terms.

License

Source code in this crate is licensed under the BSD 3-Clause License. The closed-source libempyrean runtime it loads at runtime is governed by a separate proprietary binary license; see the main repository for the dual-license breakdown.

empyrean 0.9.0