empyrean

Uncertainty-first orbit propagation, ephemeris, orbit determination, and event detection for asteroids and comets, powered by automatic differentiation

pip install empyrean

Current release: 0.11.0rc1 (release candidate); --pre is required until 0.11.0 is final.

A plain install pulls empyrean together with the B612 Foundation's pre-packaged SPICE kernels (~740 MB — see the table below). After installation, the first call to empyrean.initialize() downloads a small remainder (the moon_pa Moon-orientation kernel and the bias.dat star-catalog debiasing table — about 50 MB) that isn't available on PyPI.

Wheels are published for CPython >= 3.10 as a single abi3 stable-ABI wheel per architecture — one wheel covers CPython 3.10 and every newer version — across four platforms: macOS arm64, macOS x86_64, manylinux_2_28 x86_64, and manylinux_2_28 aarch64. There is no source distribution, so the install will not resolve on other platforms — use the other distribution channels in the meantime.

Each wheel bundles the libempyrean engine built for its own release, and the binding checks that pairing when the library opens: it reads the engine's ABI version and compares it against the one the wheel was built against, failing immediately — naming both numbers — if they differ. The version is per release, so an engine from a different release is rejected rather than used through a layout that may have moved. Nothing to configure; the bundled engine always matches. It matters only if you point the loader at an engine of your own.

Documentation lives at https://docs.empyrean-dynamics.com (a plain-text index for coding agents at /llms.txt, everything in one file at /llms-full.txt). Every release is compared across its Python, Rust, C and CLI channels and against JPL Horizons and ASSIST in the validation report, and nine runnable scenarios carry reference values to check a result against.

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, and the aberrated (light-time corrected) barycentric state per row — with sky-plane and aberrated-state covariances when the input orbit carries one.
• Orbit determination — batch-first: one determine() call fits every object in an ADES set and returns per-object results. Gauss, Herget, and systematic-ranging (admissible region + Manifold of Variations) IOD, plus a co-orbital lane that recovers Earth co-orbitals of the 2010 TK7 / 2020 XL5 class → N-body differential correction over optical and radar (delay / Doppler) observations solved jointly rather than in sequence, with span-grouped Jacobian reuse and outlier rejection. Long comet arcs deliver as full-arc fits. Solves the state — escalating to the Marsden A1/A2/A3 non-gravitational coefficients on a poor fit — plus, on the refine path, the cometary outgassing time delay DT, SRP area-to-mass, and continuous-thrust Δv corrections, all differentiated analytically, returned in a tagged solved covariance. Every deselected observation carries a typed reason. Optional post-OD H–G photometry fit recovers absolute magnitude H with an honest σ. 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

import empyrean
from empyrean import Epochs

empyrean.download_data()   # SPICE kernels, first run only
empyrean.initialize()

# Query SBDB for Apophis and propagate through its 2029 Earth flyby
orbits = empyrean.query_sbdb(["Apophis"])
epochs = Epochs.from_mjd([65000.0], scale="tdb")
result = empyrean.propagate(orbits, epochs)

# Event timeline
for i in range(len(result.events.summary)):
    ev = result.events.summary
    print(f"{ev.event_type.to_pylist()[i]:25s} "
          f"{ev.body.to_pylist()[i]:8s} "
          f"MJD {ev.epochs.mjd_tdb()[i]:.2f}")

Every time you hand to empyrean is an Epochs, and every Epochs states its scale. A bare list or array is refused: 61000.5 read as UTC and 61000.5 read as TDB are about 69 seconds apart, easily enough to move an encounter geometry, so which one you mean is stated rather than defaulted. Build them with Epochs.from_mjd(values, scale="tdb") / scale="utc", Epochs.from_jd, or Epochs.from_iso([...]) for ISO-8601 UTC timestamps.

Every time empyrean returns is an Epochs as well, in the scale you supplied. A table's primary time column is epochs and any other time column is <role>_epochs (start_epochs, ca_epochs, perihelion_epochs); one table carries one scale. Instants you asked for come back as your own values, and instants the engine found (close approaches, impacts, events, rise and set) are converted into your scale on the way out. Orbit-determination results follow the observations and come back in UTC. Two things keep a scale fixed by name instead: files written by the io writers and the CLI carry epoch_mjd_tdb columns in TDB, and the error-location diagnostic reports epoch_mjd_tdb.

Orbit determination

determine is batch-first. It groups the observations by ADES object identifier (permID, else provID, else trkSub) and fits every object, so a multi-object file is one call. Radar is optional; when supplied it is grouped by the same identifier and fitted jointly with the optical astrometry.

obs, radar = empyrean.read_ades("observations.psv")   # (optical, radar)
fits = empyrean.determine(obs, radar=radar)            # every object in the file

print(len(fits), "object(s),", len(fits.delivered), "delivered")
print(fits.orbits)      # one row per DELIVERED object, carrying object_id
print(fits.summary)     # one row per INPUT object, delivered or failed
print(fits.residuals)   # every delivered fit's rows, tagged with object_id

for object_id, failure in fits.failures.items():
    print(f"{object_id}: {failure.kind}: {failure.message}")

result = fits["2024 YR4"]                              # DetermineResult
print(
    f"converged={result.converged}, "
    f"RMS={result.summary.rms_ra_arcsec:.2f}\" RA / "
    f"{result.summary.rms_dec_arcsec:.2f}\" Dec"
)

A failed object never aborts the batch and never disappears: it gets a .summary row with status="failed", NaN measurements (never 0.0, which would read as a value at the floor), the reason in error, and a typed DetermineFailure in .failures — branch on its kind ("iod" / "od" / "radar_only" / "observer_construction" / "earth_orientation_coverage" / …), not on the message text. Seed orbits passed as initial_orbits that matched no observation group come back in .unmatched_orbit_ids rather than being dropped.

Fitting one object is the one-entry case of the same call, not a different call: fits.single() unwraps it and refuses — naming the objects — rather than choosing among several.

result = empyrean.determine(obs).single()              # one object in, one fit out

fits.orbits is a CartesianOrbits table carrying object_id, covariance, and the fitted non-grav / SRP slots — feed it straight back into propagate, generate_ephemeris, or compute_impact_probabilities.

Underneath, the fit is the engine's: optical and radar are solved jointly rather than in sequence, so a hard object's delay and Doppler tighten the same covariance the astrometry does; the co-orbital IOD lane recovers Earth co-orbitals of the 2010 TK7 / 2020 XL5 class; and long comet arcs deliver as full-arc fits. Two ODConfig switches turn those lanes off when you want the historical behaviour: coorbital_enabled=False (leaving it on does not route ordinary objects through the lane), and allow_arc_truncation=False, which makes an arc that genuinely cannot be fitted as one piece fail rather than deliver its reconcilable part with the remainder tagged outside_arc. Per-observation rejection is orthogonal and still runs.

Fit summary

fits.summary is a FitSummary quivr table with one row per input object, so a partially successful batch is readable rather than silently shorter than its input. The column names match the fit_summary.parquet / fit_summary.csv files the CLI writes, so a table read back off disk and this one describe a fit identically.

Both acceptability verdicts are columns on it, and every extrapolation gate sits beside the value it measured and the threshold it was compared against — so a False says which axis failed and by how much:

s = fits.summary

print(s.object_id.to_pylist(), s.status.to_pylist())   # "delivered" / "failed"
print(s.fit_acceptable.to_pylist())
print(s.extrapolation_acceptable.to_pylist())

def col(name):     # value columns are nullable — NaN where not computable
    return s.column(name).to_numpy(zero_copy_only=False)

print(col("selection_fraction"),    col("selection_fraction_threshold"))
print(col("selected_arc_fraction"), col("selected_arc_fraction_threshold"))
print(col("trailing_gap_days"),     col("trailing_gap_threshold_days"))
print(col("fractional_sigma_a"),    col("fractional_sigma_a_threshold"))

The same verdict in full structured form is result.acceptability, an AcceptabilityReport whose every *_ok flag is paired with its *_value and *_threshold — plus radar_fit_ok, a tri-state that is None when no radar contributed, which is never the same as False. Tune the bounds with AcceptabilityThresholds on ODConfig.

Per-observation residuals

result.observations (and, batch-wide, fits.residuals) is an ObservationResults table carrying the whole 35-column residual surface, not a projection of it: the obs_id / object_id join keys, observatory code, catalog and epoch; the RA/Dec residuals with their effective combined covariance (residual_cov_ra / residual_cov_dec / residual_cov_corr); the complete rejection block; the influence diagnostics; the along/cross-track decomposition with its full 2×2 covariance; and the radar block. Radar rows carry the delay / Doppler residual (observed − predicted in seconds / hertz) with its χ², dof, survival probability, and combined variance instead of RA/Dec.

Nothing is deselected silently — rejection_reason names the layer that dropped each row (chi_squared, sigma_clip, cooks_distance, adaptive, cmc2003, unsupported_observatory, outside_arc, non_finite_chi2, missing_jacobian, …) alongside the criterion value, the static threshold, and the effective threshold it was tested against. influence_information_loss is the D-optimality information loss on removal (+inf marks an indispensable observation).

print(result.observations.rejected_only().rejection_reason.to_pylist())
print(result.observations.worst_chi2(5).obs_id.to_pylist())
print(fits.residuals.select_station("F51").rms_combined_arcsec)

result.covariance_trust is an event-aware verdict on the delivered covariance: trusted, encounter_intervenes (naming the intervening close-approach or high-nonlinearity event and whether a second-order state-only correction can recover it), or weakly_determined_high_n. It is None when no trust gate ran — absence of a verdict is not trust.

How the solver stopped

Every fit says how its solve ended, so a converged orbit and one the solver merely stopped producing steps for are no longer the same success:

print(result.termination)            # "gradient_tolerance", "stalled_delivered", ...
print(result.gn_step_qnorm)          # undamped Gauss-Newton step qnorm, or None
print(result.accepted_steps)         # 0 is a real reading, not "unknown"
print(result.final_solve_iterations, result.mu_final)
print(result.stall_delivery)         # StallDelivery, or None

termination is a stable lowercase tag — "gradient_tolerance", "step_tolerance", "cost_tolerance", "max_iterations", "damping_exhausted", "inner_trials_exhausted", "stalled_delivered", "schur_step_tolerance" — or "unrecognized" for a stop this engine build cannot name. None means nothing was reported, which is never the same claim as "unrecognized": a stop that fired but cannot be named is not no stop at all.

gn_step_qnorm is the undamped Gauss-Newton step's quadratic form at the delivered iterate — the quantity convergence_tol bounds, and the only step norm comparable to it. update_norm is the μ-damped last accepted step and is not comparable to convergence_tol; read termination and gn_step_qnorm for convergence. update_norm of exactly 0.0 beside accepted_steps == 0 means the solver latched a criterion at its starting point and never moved — the ordinary outcome of refitting an already-stationary orbit, and previously indistinguishable from a step of size zero having been taken. stall_delivery is populated exactly when termination is "stalled_delivered", carrying the solver's own verdict underneath and the numbers that earned the delivery.

Iterating on one arc: Session

Session is the stateful counterpart to the one-shot call — it owns one observation set, its mask state, and the fit history, so "mask a night, refit, compare" is three calls rather than three rebuilds. Construct it from a path, PSV text, or an ADESObservations table (no PSV round trip):

import empyrean

obs, _radar = empyrean.read_ades("apophis.psv")

sess = empyrean.Session.from_observations(obs)     # or Session("apophis.psv")
sess.refine()                                      # initial fit -> history[0]

for i, code in enumerate(obs.stn.to_pylist()):
    if code == "T05":
        sess.mask(i)
refit = sess.refine()                              # refit without T05

print(sess.n_observations, sess.n_masked, sess.n_active)

diff = sess.diff(0)                                # current vs history[0]
print(f"Δχ²/dof = {diff.reduced_chi2_delta:+.3f}")
print(f"Δn_obs  = {diff.n_observations_delta:+d}")
print(f"stopped on {refit.termination}")

mask / unmask / unmask_all / is_masked move the mask; n_observations / n_masked / n_active and history_len are properties; history(i) reads any earlier DetermineResult back whole. diff(prior_idx) returns a SessionDiff carrying reduced_chi2_delta, iterations_delta, n_observations_delta, and the two update_norm_* values — read those last two as a damping-trajectory diagnostic only, never as convergence, for the reason above.

Weighting

ODConfig ships production defaults: the VFCC2017 weighting preset (Vereš, Farnocchia, Chesley & Chamberlin 2017 station floors) plus per-night 1/√N de-weighting, and EFCC2020 catalog debiasing (Eggl, Farnocchia, Chamberlin & Chesley 2020). WeightingConfig replaces or extends the layer chain — first-match-wins, so an entry in additional_layers overrides the preset for its stations and the preset is the fallback.

from empyrean import (
    ODConfig, WeightingConfig, WeightingLayer, WeightingLayerKind, WeightingPreset,
)

config = ODConfig(
    weighting=WeightingConfig(
        preset=WeightingPreset.VFCC2017,
        additional_layers=[
            WeightingLayer(
                kind=WeightingLayerKind.OBSERVATORY_RULE,
                obs_code="F51",
                sigma=(0.1, 0.1),
            ),
            # additional_layers REPLACES the default list — re-include the
            # nightly layer to keep production behavior.
            WeightingLayer(kind=WeightingLayerKind.NIGHTLY_DEWEIGHTING),
        ],
    )
)

Wide-parameter fitting

A fit solves the 6-element state by default, escalating to the Marsden A1/A2/A3 non-gravitational coefficients on a poor fit. SolveFor on ODConfig.solve_for_flags requests an explicit wider solve: beyond state + Marsden, determine and refine can also solve for the cometary outgassing time delay dt, the solar-radiation-pressure area-to-mass ratio amrat, and per-segment thrust Δv corrections (thrust) — each differentiated analytically by the same hyperdual integrator that drives the dynamics.

Each axis takes a disposition, not a flag, because "not solved" is two different answers:

A considered axis is not a safety margin. Under an uncorrelated prior the consider correction strictly widens the posterior, but the fits that need it are the ones with cross terms between the considered axis and the solved ones — and there the correction is sign-indefinite, so the posterior can come back tighter. Report it as an unestimated error source folded through its measurement partials, never as conservatism (Schmidt–Kalman consider analysis; Tapley, Byron D., Schutz, Bob E., and Born, George H., Statistical Orbit Determination, Elsevier Academic Press, 2004, ch. 6).

False cannot say which of the last two was meant, so a bool is refused by name rather than coerced. result.dispositions reports the partition the fit actually ran — the partition resolved against the orbit, not the one requested, so an AUTO escalation is readable after the fact. It is also the only place a considered axis appears, since a solved covariance's slot tags record what occupied a column and a considered axis occupies none. That is what tells you whether re-attaching a prior to an axis would double-count it: a considered axis already has its uncertainty inside the delivered 6×6, a fixed one does not. Same covariance, opposite conclusions.

result.warnings is the other half of that honesty: covariance the fit was handed and deliberately did not use, delivered as payload rather than written to a log, because a dropped prior cross term changes how the σ for that slot should be read. It is empty on a fit that used everything it was given, which is the common case.

dt, amrat, and thrust are refine-path solves: the seed orbit must carry the prior that opens each axis, so run them through refine. Each prior rides the orbit's packed covariance joint (orbit.covariance) on the axis it opens — the DT axis for DT (days²), the Marsden block for the non-gravitational coefficients, the AMRAT axis for AMRAT — so a fitted orbit already carries what reopens the same axis on a refit. Requesting an axis whose prior is absent is rejected loudly — the fit never returns a zeroed or defaulted column.

from empyrean import ODConfig, ParamDisposition, SolveFor

# Solve state + Marsden A1/A2/A3 + the outgassing time delay DT. The
# seed orbit's packed covariance joint carries the priors that open each
# axis — a Marsden block (opens Marsden) and a DT-axis variance (opens DT).
config = ODConfig(solve_for_flags=SolveFor(marsden="solved", dt="solved"))
result = empyrean.refine(orbit, obs, config=config)

print(result.dt_delta)      # fitted ΔDT (days); None if DT was not solved
print(result.amrat_delta)   # fitted ΔAMRAT (m²/kg); None if not solved
print(result.dispositions)  # what the fit did with each axis

# The enum is equivalent to the string form, and is what `dispositions`
# reports back:
ParamDisposition.parse("considered") is ParamDisposition.CONSIDERED

print(result.warnings)      # covariance the fit declined to use; often []

# Per-segment thrust dispositions are positional with the orbit's
# declared Δv-correction segments — a considered or fixed burn sits
# between solved ones as readily as after them, so a count could not
# say which burn is which. At most three entries (the engine's
# MAX_THRUST_SEGMENTS, on empyrean.od.result); a longer list raises
# rather than being truncated to a shorter fit that drops a burn
# silently. A bool raises too, by name — it cannot say which of
# "considered" and "fixed" was meant.
SolveFor(thrust=["solved", "fixed", "solved"])

Tagged solved covariance

A wide fit returns a SolvedCovariance on result.solved_covariance whose fitted-parameter identities travel with the matrix. Read a parameter's variance by its slot — never by guessing column order:

sc = result.solved_covariance          # None for a state-only fit
if sc is not None and sc.dt_slot is not None:
    dt_var = sc.matrix[sc.dt_slot, sc.dt_slot]   # DT variance (days²)
    print(f"σ(DT) = {dt_var ** 0.5:.4f} days")
# sc.marsden_slot / sc.amrat_slot / sc.thrust_slots locate the rest;
# canonical layout is [state 6 | Marsden 3 | DT 1 | AMRAT 1 | thrust 3×k].

Carrying the joint onward

The fitted orbit carries the whole matrix, not just its diagonal blocks, so it feeds straight back into propagation and every downstream entry point. It rides on orbit.covariance — a PackedJoints column whose per-row joint(i) is a PackedJoint holding the layout, the packed lower triangle, and the tags (kind, quality, functional, provenance, and an asserted joint's asserted_source). The joint is the one authority; the flat 6×6 state covariance is joint.state_block(), a derived read-only view of its leading block, never a second, independently settable home.

joint = result.orbit.covariance.joint(0)   # the fit's own PackedJoint, or None
if joint is not None:
    block = joint.state_block()             # the derived 6×6 state view
    sigma_x = block[0, 0] ** 0.5
    width = joint.width                     # 6 + solved parameters

orbit.covariance is a sub-table column, never None on a parent table — quivr returns a table of parent length regardless. A row with no covariance is a per-row absence: test orbit.covariance.row_is_empty(i) or let joint(i) return None, never orbit.covariance is None.

Carrying it onward is passing the orbit through whole — there are no halves to reassemble. propagate, compute_impact_probabilities, compute_b_planes, generate_ephemeris, and the determine / evaluate / refine seed path all condition on the joint when the orbit carries it. An impact probability computed against a block-diagonal covariance materially understates the tails: it asserts an independence the fit never found, and the propagated joint has non-zero state↔parameter columns even from a block-diagonal input, because propagation itself generates the correlation — so a second leg handed only the 6×6 reports a tighter uncertainty than the first leg supports.

Moving a covariance between representations, frames or origins goes through transform_orbits, whose Jacobian transports the whole joint; the bare-state transform_coordinates refuses a covariance-bearing input rather than drop it silently, raising a PushError (PushErrorVariant.TRANSFORM_COVARIANCE_UNSUPPORTED) that points at transform_orbits.

Propagated states carry the same joint per output epoch. With tagged_covariance=True, result.tagged_covariance_series(i) and the flat result.tagged_covariance table hold one PackedJoint per epoch on their own covariance column; tc.joint is the authority and tc.matrix() its derived state 6×6. It is populated on every uncertainty method that produces a joint, including the sampled ones, which recover the state↔parameter columns from the propagated cloud.

result = empyrean.propagate(orbits, epochs, tagged_covariance=True)

tc = result.tagged_covariance_series(0)[-1]   # last epoch of orbit 0
joint = tc.joint                              # the PackedJoint; tc.matrix() is its 6×6

The joint crosses the compiled-extension boundary as a flat wire dict: joint.to_wire() serializes it and PackedJoint.from_wire(d) reads it back (PackedJoints.to_wire_list() / from_wire_list(rows) for a whole column). A malformed joint — a bad layout, a triangle of the wrong length, a retired discriminant — is refused there by a PushError carrying the PushErrorVariant that names the cause, never accepted and quietly turned into a wrong answer.

Post-OD photometry

Attach a PhotometryConfig to recover the absolute magnitude H and a phase-function slope from the observation magnitudes after the orbit is solved — the 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 reports the model it actually fitted on model_used. H comes back with an honest 1σ from the fit covariance. Each reported magnitude is converted to Johnson V before the fit with the MPC Magnitude Bands table; PhotometryConfig.band_conversions supplies entries the table does not carry, as {code: (offset_to_v, sigma)}. The LSST codes Lu Lg Lr Li Lz Ly carry no MPC conversion and are refused by name until supplied — their magnitudes are counted (n_mags_dropped_unconvertible) and their codes listed (dropped_bands), and the observations' astrometry is unaffected. Each converted band's offset is reported on the result's per-band offset_applied.

from empyrean import ODConfig, PhotometryConfig

config = ODConfig(photometry=PhotometryConfig())   # AUTO ladder, MPC table

# Supply a V conversion for an LSST band so its magnitudes enter the fit:
config = ODConfig(
    photometry=PhotometryConfig(band_conversions={"Lu": (-1.41, 0.30)})
)
result = empyrean.determine(obs, config=config).single()

phot = result.photometry               # None if photometry was not requested
if phot is not None and phot.covariance is not None:
    sigma_h = phot.covariance[0, 0] ** 0.5
    print(f"H = {phot.h:.2f} ± {sigma_h:.2f} mag  (model {phot.model_used.value})")

determine folds the fit onto the result orbit's photometric column — including the 3×3 covariance over (H, slope1, slope2) — so a fitted orbit handed straight to generate_ephemeris reports the H uncertainty in mag_sigma:

[ \sigma_V = \sqrt{\sigma^2_\text{photo} + \sigma^2_\text{state}} ]

with (\sigma^2_\text{photo} = J \Sigma_p J^\top) contracting the full 3×3 against (J = [\partial V/\partial H, \partial V/\partial \text{slope}_1, \partial V/\partial \text{slope}_2]).

Because (V = H + 5\log_{10}(r\Delta) + \phi(\alpha)) gives (\partial V/\partial H \equiv 1) exactly, an orbit carrying no state covariance and a covariance of the H-only shape (\mathrm{diag}(\sigma_H^2, 0, 0)) reports (\sigma_V = \sigma_H):

import numpy as np
from empyrean import PhotometricParams

orbits = orbits.set_column(
    "photometric",
    PhotometricParams.from_kwargs(
        model=["hg"],
        h=[19.7],
        g=[0.15],
        covariance=[np.diag([0.3**2, 0.0, 0.0]).reshape(9).tolist()],
    ),
)

Slope variances and (H)–slope covariances do not drop out: they contract against (\partial V/\partial\text{slope}), which vanishes only at zero phase angle, so a covariance carrying them reports (\sigma_V > \sigma_H). SBDB is the common case — a queried orbit carries (\mathrm{diag}(\sigma_H^2, \sigma_G^2, 0)), so for Apophis ((\sigma_H = 0.19), (\sigma_G = 0.11)) the reported (\sigma_V) grows with phase angle and is strictly above 0.19 at every nonzero phase angle. The column is a row-major flattened 3×3 over ((H, \text{slope}_1, \text{slope}_2)); it must be symmetric, and it is null when absent, never a block of zeros:

orbits = orbits.set_column(
    "photometric",
    PhotometricParams.from_kwargs(
        model=["hg"],
        h=[19.7],
        g=[0.15],
        covariance=[np.diag([0.19**2, 0.11**2, 0.0]).reshape(9).tolist()],
    ),
)

The two terms are combined as independent. They are not strictly independent — a fitted (\sigma_H) is conditional on the fitted state, because the photometric fit holds the geometry exact — and no joint state↔photometry covariance is computed anywhere in the stack, so there is no cross term to add. The resulting (\sigma_V) is mildly conservative, which is the safe direction.

Ephemeris

from empyrean import Frame, Origin

observers = empyrean.get_observer_states(["W84", "F51"], epochs)
eph = empyrean.generate_ephemeris(orbits, observers)

# Observer states default to the ICRF / SSB construction basis, which is
# what ephemeris generation and orbit determination require and which is
# returned untransformed. Pass frame= / origin= when you want the sites
# somewhere else — e.g. heliocentric ecliptic, for geometry plots:
sites = empyrean.get_observer_states(
    ["W84", "F51"], epochs, frame=Frame.ECLIPTICJ2000, origin=Origin.SUN
)

print(eph.ephemeris.coordinates.lon.to_numpy())   # RA (degrees)
print(eph.ephemeris.coordinates.lat.to_numpy())   # Dec (degrees)
print(eph.ephemeris.mag.to_numpy(zero_copy_only=False))  # apparent V (nullable)

# Orbits carrying a covariance also get, per row, the 6×6 sky-plane
# covariance over (rho, RA, Dec + rates) in AU / degree units, and the
# aberrated (light-time corrected) barycentric ICRF state at the
# photon-emission epoch with its own 6×6 covariance:
print(eph.ephemeris.coordinates.covariance.to_matrix().shape)      # (N, 6, 6)
print(eph.ephemeris.aberrated_state.covariance.to_matrix().shape)  # (N, 6, 6)

That sky covariance is the marginal over whatever force-model parameter uncertainty the orbit declares, not the conditional: an orbit declaring none is unaffected, and any Marsden 3×3, DT / AMRAT variance or Δv block widens every sky σ.

eph.warnings lists non-fatal generation warnings — e.g. an Earth-orientation kernel coverage gap handled by the analytic IAU 2006 fallback, or rows whose sensitivity chain was skipped — naming the affected orbit / observatory / epoch. Empty when the run had nothing to report.

Self-perturbing targets

All sixteen SB441-N16 bodies (1 Ceres, 2 Pallas, 4 Vesta, 7 Iris, …) are simultaneously members of the Standard force model and legitimate objects to propagate. ephemeris_overlap_policy says what the engine does when a target coincides with its own perturber: SUBSTITUTE_SPK (the default) returns the body's authoritative SPK states and integrates nothing — so there is no dense trajectory, no STM, and no sensitivity chain, and ephemeris generation for such a body fails outright. EXCLUDE_AND_INTEGRATE drops the overlapped perturber, integrates the caller's own initial conditions, and reports the overlap. It is what generating an ephemeris for an SB441-N16 body at Standard tier requires.

from empyrean import EphemerisConfig, EphemerisOverlapPolicy, PropagationConfig

config = EphemerisConfig(
    propagation=PropagationConfig(
        ephemeris_overlap_policy=EphemerisOverlapPolicy.EXCLUDE_AND_INTEGRATE,
    )
)
ceres = empyrean.query_sbdb(["Ceres"])
eph = empyrean.generate_ephemeris(ceres, observers, config=config)

Naming the body in excluded_perturbers is the other escape. Use one or the other: overlap detection is what decides whether the policy is consulted at all, and it is suppressed whenever an explicit exclusion list is present.

Uncertainty

from empyrean import UncertaintyMethod

# Second-order: populates STM (6x6) and STT (6x6x6)
result = empyrean.propagate(
    orbits, epochs,
    uncertainty_method=UncertaintyMethod.SECOND_ORDER,
)
print(result.sensitivity.stms_array().shape)   # (N, 6, 6)
print(result.sensitivity.stts_array().shape)   # (N, 6, 6, 6)

Reading back the mixture

Under GAUSSIAN_MIXTURE — and inside AUTO's close-approach windows — the engine splits the input Gaussian and retains the resulting components at every close approach where the splitter actually fired. Those components come back on PropagationResult.mixtures, a flat MixtureChains table with one row per (orbit, CA epoch, component). They are the mixture itself, not its moment collapse: you can evaluate (\sum_k w_k,\mathcal{N}(x \mid \mu_k, \Sigma_k)) directly at the CA epoch.

from empyrean import GaussianMixture

result = empyrean.propagate(orbits, epochs, uncertainty_method=GaussianMixture())

# One list per retained CA epoch of orbit 0.
for group in result.mixture_chains(0):
    print(len(group), "components,", sum(c.weight for c in group), "retained weight")
    print(group[0].mean, group[0].covariance.shape, group[0].frame, group[0].origin)

Four limits on what is retained, each a real property of the engine's retention rather than a marshaling shortfall:

• Depth-0 only. Only the initial split is retained; recursive AGM calls (depth > 0) are not captured.
• Only CA epochs where AGM fired. An orbit that never triggered a split contributes no rows — a FIRST_ORDER run yields an empty table, never zero-filled placeholder rows.
• Component covariance is the linear map. Each component's covariance is (\Phi \Sigma_k \Phi^\top); the second-order mean correction is intentionally omitted.
• Retained weights may sum to less than 1. A sub-Gaussian whose own sub-propagation missed the close approach contributes no component, and the deficit is not recorded anywhere. Do not assume the weights normalize — sum them and check.

Note the name: empyrean.propagation.mixtures.MixtureComponent is the basis-tagged read-back component, a different type from the top-level empyrean.MixtureComponent, which is the split_gaussian primitive at (t_0).

Continuous thrust

Model finite burns / low-thrust arcs by passing one ThrustParams per orbit through propagate's thrust_arcs keyword (None for the ballistic orbits). 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.

import empyrean
from empyrean import Origin
from empyrean.orbits.thrust import ConstantRTN, ThrustArc, ThrustParams

# 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.
arc = ThrustArc(
    start_mjd=65000.0,
    end_mjd=65010.0,
    scale="tdb",
    thrust_n=1.0,
    mass_kg=500.0,
    steering=ConstantRTN(alpha_rad=0.0, beta_rad=0.0),
    sharpness=100.0,
    central_body=Origin.SUN,
    isp_s=3000.0,
)

# One entry per orbit, positionally aligned with `orbits`. Add per-arc Δv
# targeting corrections with ThrustParams(arcs=[arc], dv_corrections=[...]).
result = empyrean.propagate(orbits, epochs, thrust_arcs=[ThrustParams(arcs=[arc])])

System handles

Assembling the force model has a fixed per-call cost. build_system assembles it once for a frozen {force model, frame, encounter-timescale divisor} key and returns a BuiltSystem you reuse across many propagations — the build-once, propagate-many pattern for short-arc campaigns. Its propagate / generate_ephemeris release the GIL, so the handle can be shared across threads. A call that disagrees with the frozen key is rejected loudly, never silently rebuilt; rebuild the handle after any initialize() / data reload.

The same handle serves orbit determination: determine / evaluate / refine mirror the module-level empyrean.determine / empyrean.evaluate / empyrean.refine argument for argument, which is the measure-and-extend loop — assemble once, refine many. Build that handle with od_system(force_model) rather than build_system: a fit picks its own integration frame and encounter-timescale divisor, and the frame is not the one the propagation examples freeze, so a frame="icrf" handle is refused on the OD path with a frame-mismatch error. ODConfig.auto_force_model is refused too — it lets the fit re-pick its own tier part-way through, which no frozen handle can follow.

import empyrean
from empyrean import ForceModelTier, Frame

# Build once for the Standard model in the ecliptic frame. force_model and
# frame accept the enums or their string / int forms.
system = empyrean.build_system(ForceModelTier.STANDARD, Frame.ECLIPTICJ2000)

result = system.propagate(orbits, epochs)

# describe() is the reproducibility record: the force-model menu plus the
# identity of every loaded kernel (SHA-256 for file-backed kernels; the
# model name for built-in fields).
desc = system.describe()
print(len(desc.perturber_origins), "perturbers,", len(desc.kernels), "kernels")

# For fitting, freeze the recipe a fit actually runs under and reuse one
# handle across the whole measure-and-extend loop.
od = empyrean.od_system(ForceModelTier.STANDARD)
for arc in arcs:
    fit = od.refine(orbit, arc)  # one force-model assembly total

Impact probability and B-plane geometry

For each detected close approach, you can ask the propagator for an impact-probability assessment or a full B-plane breakdown — and run several uncertainty methods side-by-side on the same encounter:

import pyarrow.compute as pc

from empyrean import UncertaintyMethod

end_epoch = empyrean.Epochs.from_mjd([63000.0], scale="tdb")

ips = empyrean.compute_impact_probabilities(
    orbits,
    end_epoch=end_epoch,
    methods=[UncertaintyMethod.FIRST_ORDER, UncertaintyMethod.SECOND_ORDER],
)
ips.epochs.scale                    # "tdb"
second = ips.where(pc.field("method") == "second_order")
second.ip_second_order.to_numpy(zero_copy_only=False)   # nullable
ips.ip_linear.to_numpy()            # always populated

bps = empyrean.compute_b_planes(orbits, end_epoch, [UncertaintyMethod.SECOND_ORDER])
print(bps.b_dot_t_km.to_numpy())    # B·T (km)
print(bps.b_dot_r_km.to_numpy())    # B·R (km)
print(bps.semi_major_3sig_km.to_numpy(zero_copy_only=False))  # 3σ semi-major

Returns typed ImpactProbabilities and BPlanes quivr tables — one row per (method × orbit × body) encounter, with the closest-approach time as an embedded Epochs sub-table so .to_utc() / .to_tdb() just works.

Each entry carries its own parameters — methods in one call never share them. Auto is the adaptive method: it escalates the covariance treatment over each close-approach window and relaxes it elsewhere, choosing among first-order, second-order, and the adaptive Gaussian mixture. Pass the Auto dataclass instead of the enum to tune its κ band edges and mixture knobs, and those are the values the engine runs with:

from empyrean import Auto, MonteCarlo

ips = empyrean.compute_impact_probabilities(
    orbits,
    end_epoch=empyrean.Epochs.from_mjd([63000.0], scale="tdb"),
    methods=[
        Auto(threshold_first=0.05, threshold_mixture=5.0, gmm_max_depth=4),
        MonteCarlo(n_samples=100_000, seed=7),
    ],
)
# ip_agm is populated when the mixture refinement fired; a null means it
# was a no-op, so use ip_agm when finite and ip_linear otherwise. Never
# back-fill it — the null is what distinguishes "ran" from "no-op".
print(ips.ip_agm.to_numpy(zero_copy_only=False))
print(ips.agm_components.to_numpy(zero_copy_only=False))

UncertaintyMethod.AUTO (or "auto") selects the same method with the engine-default thresholds; a default-constructed Auto() is identical to it.

Each ImpactProbabilities row also carries the geodetic impact point (latitude / longitude / altitude on the body's reference ellipsoid; null when no surface projection is available), the 95% binomial confidence half-width on ip_mc, the second-order corrected mean miss distance, 1σ miss-distance 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.

Observation planning

Given an orbit that already carries a covariance, evaluate_plan ranks candidate follow-up observations by how much each one would tighten it — before you spend the telescope time:

import empyrean
from empyrean import (
    CartesianCoordinates,
    Epochs,
    PlannedObservation,
    transform_coordinates,
)

# The planner consumes a Cartesian, barycentric covariance and converts
# neither, so both are required. One call does both; the origin half is a
# pure translation, so the covariance and its metrics come across unchanged.
coords = transform_coordinates(fit.orbit.coordinates, CartesianCoordinates, origin="SSB")
orbit = fit.orbit.set_column("coordinates", coords)
t0 = float(coords.epochs.mjd_tdb()[0])

plan = empyrean.evaluate_plan(
    orbit,
    [
        PlannedObservation.optical(Epochs.from_mjd(t0 + 30.0, scale="tdb"), "F51", (0.2, 0.2)),
        PlannedObservation.optical(Epochs.from_mjd(t0 + 31.0, scale="tdb"), "568", (0.3, 0.3)),
        PlannedObservation.radar(Epochs.from_mjd(t0 + 45.0, scale="tdb"),
            radar_bandwidth_hz=1.0e5,
            radar_freq_resolution_hz=0.1,
            radar_snr=50.0,
        ),
    ],
)

plan.metrics.to_dataframe()                 # two rows: "prior", "posterior"
plan.metrics.prior().position_sigma_km[0].as_py()      # before any candidate
plan.metrics.posterior().position_sigma_km[0].as_py()  # after all of them
plan.candidates.best_by_information_gain(3).obs_code.to_pylist()

plan.metrics is a two-row PlanMetrics table with a stage column ("prior" / "posterior") plus the five covariance summary metrics — RSS position and velocity σ, the 1σ position-ellipsoid semi-axes, and log_det. Keeping it a table means plans concatenate and join like any other output; prior() / posterior() are the one-row views.

plan.candidates is a PlanCandidates quivr table — one row per candidate, carrying marginal_volume_reduction (the per-dimension generalized-variance ratio (det Σ_post / det Σ_prior)^(1/6) over the 6×6 state covariance — a D-optimality score normalized to one dimension, so it reads as a linear scale factor; the 1σ ellipsoid volume ratio is that value cubed), the fractional position-σ improvement, and the running covariance metrics after that candidate and every one folded before it. plan.ephemeris is the predicted sky position at each optical candidate's epoch, with the epoch as an embedded Epochs sub-table; an optical row's index is its row there. A radar row has no epoch of its own — its index is its rank among the radar candidates, ordered by epoch.

Rows come back in ascending epoch order, and each candidate's marginal gain is measured against the covariance that already contains every earlier one. The gains are therefore conditional: two identical observations do not score identically, and best_by_information_gain ranks contributions within one campaign rather than standalone candidate value. To compare candidates head to head, evaluate a one-candidate plan for each.

observable means different things per kind. On an optical row it is a real engine verdict — today a solar-elongation test and nothing else, since the target's absolute magnitude does not reach the planner, so the limiting-magnitude filter cannot fire. On a radar row it is always True: no radar feasibility test runs here, not even the antenna elevation limit, so True means "not assessed". Either way the filters are engine-set and not caller-configurable, and observable does not gate the fold — an unobservable candidate still contributes to posterior and to every later cumulative_* row. What does gate it is whether the engine could compute the candidate's observation partials; a candidate for which it could not reports a marginal_volume_reduction of exactly 1 and leaves the covariance untouched. observable_only() filters rows, not information; to price the observable subset, drop those candidates from the plan and evaluate again.

Not exposed in this release, recorded so the omissions are not mistaken for oversights: the non-gravitational planning variant that solves over state ⊕ (A1, A2, A3) and reports the σ(A2) tightening a radar campaign buys, batch plan evaluation across many orbits, and the encounter B-plane. An orbit carrying non-gravitational parameters is accepted and evaluated state-only — the acceleration still acts in the dynamics, but the solve-for set stays 6×6 and no σ(A2) is reported.

Observability (the night product)

compute_visibility surveys an orbit against a set of ground sites and reports, per site and per local night, the twilight epochs, the local solar midnight, the observable windows, and an anchor epoch under each policy — nearest to midnight, midnight, transit, or an explicit UTC — with its circumstances: altitude, azimuth, the Moon's altitude, azimuth and illuminated fraction, the predicted magnitude, and the sky covariance. Unlike evaluate_plan, it reads the whole ObservatoryConfig — the elevation and darkness gates and the twilight thresholds — so a VisibilityPoints row's observable is the conjunction of all four site gates (elongation, magnitude, elevation, darkness), the site-aware verdict plan evaluation cannot give.

import empyrean
from empyrean import Epochs, ObservatoryConfig, PlanningConfig, compute_visibility

config = PlanningConfig(
    observatories=[
        ObservatoryConfig(
            obs_code="F51",                 # Pan-STARRS 1, Haleakalā
            sigma_arcsec=(0.2, 0.2),
            max_apparent_mag=24.0,
            min_elongation_deg=60.0,
            min_elevation_deg=20.0,
        )
    ]
)

# `orbits` is any Orbits table carrying a covariance (e.g. a fit's
# `result.orbit`, or query_sbdb output). Survey a three-day window hourly.
t0 = float(orbits.coordinates.epochs.mjd_tdb()[0])
# The survey window is a scale-carrying input, like impact's end_epoch; the
# result epochs come back in the window's scale.
start = Epochs.from_mjd(t0, scale="tdb")
end = Epochs.from_mjd(t0 + 3.0, scale="tdb")
vis = compute_visibility(orbits, start, end, 1.0 / 24.0, config)

vis.nights.to_dataframe()      # twilight, local solar midnight, windows
vis.anchors.to_dataframe()     # one Moon-aware anchor per policy, with its joint
vis.points.to_dataframe()      # per-epoch gates + Moon geometry

compute_visibility returns a VisibilityResult with five quivr tables — points, windows, nights, anchors, and the per-site observatories — joining on orbit_id / obs_code (and night_index for a night's windows and anchors). Each anchor carries its sky covariance as a packed joint, so an anchor's on-sky uncertainty reads back the same way a fitted orbit's does.

Optical and radar candidates fill different blocks and the cross-block is null, never zero: sky-plane geometry (along_track_sigma_arcsec, position_angle_deg, …) is populated only on optical rows, and radar_mode / radar_snr / radar_range_km only on radar rows. Radar adds the line-of-sight range and range-rate that angles-only astrometry cannot supply; its measurement σ is the Cramér-Rao bound set by the waveform bandwidth and the effective SNR.

Leave radar_snr at None and the SNR is derived from a link budget over the target's physical properties instead. That path never substitutes a value it was not given — a missing property it needs is a loud refusal, and anything it derived or adjusted comes back on radar_provenance. Supplying an SNR selects the other request shape, so the link-budget inputs are refused alongside it rather than dropped:

plan = empyrean.evaluate_plan(
    orbit,
    [
        PlannedObservation.radar(Epochs.from_mjd(t0 + 45.0, scale="tdb"),
            radar_bandwidth_hz=1.0e5,
            radar_freq_resolution_hz=0.1,
            radar_target_h_mag=19.7,
            radar_target_visual_albedo=0.23,
            radar_target_radar_albedo=0.15,
            radar_integration_s=600.0,
        )
    ],
)
plan.candidates.radar_provenance.to_pylist()
# [['diameter derived from H + p_V',
#   'spin period unknown — coherent integration uncapped']]

Reading and writing files

empyrean.io writes orbits, ephemerides, events, and residuals to parquet / JSON / CSV, and reads orbits back.

empyrean.io.write_orbits_parquet("orbits.parquet", fits.orbits)
empyrean.io.write_orbits_csv("orbits.csv", fits.orbits)
empyrean.io.write_residuals_parquet("residuals.parquet", fits.residuals)

orbits = empyrean.io.read_orbits_parquet("orbits.parquet")

Neither writer projects. The residual writers marshal every column of ObservationResults across the boundary — join keys, rejection attribution, influence diagnostics, sky-motion decomposition and radar block included — for a 36-column file (the table's 35 plus a has_radar flag), identical across all three formats; rejection_reason and radar_kind are written as names rather than integer codes, and a non-computable number is a literal NaN in CSV and null in JSON.

Parquet and JSON round-trip the orbit's full packed joint; CSV carries the state 6×6 only. The packed joint rides the orbit as its covariance column — the layout, the packed lower triangle and the tags — so a fitted orbit written to parquet or JSON round-trips carrying the whole joint the fit computed rather than its diagonal blocks, and reads back ready to re-feed. CSV is a flat scalar row shape that cannot express the layout distinguishing an absent axis from a present one, so it carries the state 6×6 (the joint's leading block) and refuses a joint carrying solved-parameter axes by name rather than writing a projection that would read back as a different orbit.

Turning event detection off

Every propagation runs the built-in detector set on each accepted integrator substep. The five per-type flags on EventConfig filter what gets emitted; they do not stop the detectors running. When you want states and nothing else, detection_enabled is the switch that stops the work:

config = empyrean.PropagationConfig(
    events=empyrean.EventConfig(detection_enabled=False),
)
result = empyrean.propagate(orbits, epochs, config=config)
assert len(result.events.summary) == 0

Measured on a 64-orbit, covariance-free, two-epoch Standard-tier batch, single-threaded over a 400-day arc: 272 ms with detection on, 185 ms with it off — 1.47x.

State accuracy is unchanged — the propagated states, their covariance and the STM come back bit-for-bit identical either way, because detection is observation and not dynamics. What is gone is everything the detectors produce: no events, no close approaches, and therefore no impact probabilities.

"auto" and "gaussian_mixture" resolve themselves from that output, so pairing either with detection_enabled=False raises rather than returning a silently linear answer. Every other method is served normally.

EventConfig also carries dense_origin ("bodycentric" or "barycentric") and capture_criterion ("population", "individual" or "energy_only"); both were silently dropped at the boundary until now, and an unrecognised value for either is now refused rather than ignored.

When a batch fails, it names the orbit

A batch call takes N orbits and M epochs and fails as a whole. When the failure belongs to one orbit, the exception says which — so the offending row is read off the failure instead of found by re-running the batch one orbit at a time.

try:
    result = empyrean.propagate(orbits, epochs)
except Exception as e:
    print(e.orbit_index, e.orbit_id, e.epoch_mjd_tdb)   # 2317 '2024 YR4' 60800.5

The three attributes are set on every exception raised from an engine failure, None included, so a caller reads them directly rather than guarding each access with getattr. They are None together when the failure belongs to no single orbit — an empty epoch grid, a missing kernel, a config the whole call was refused on. orbit_id is the id from the batch you passed, so it joins directly against the rows you sent.

Data files

empyrean needs a set of SPICE kernels. Most arrive via PyPI as installation dependencies; the remainder download on first use.

From pip (installed automatically with empyrean)

empyrean bundles gm_de440.tpc (12 KB) in the wheel itself. On initialize(), empyrean stages symlinks to these files in the platform data directory (~/.local/share/empyrean/data/ on Linux, ~/Library/Application Support/empyrean/data/ on macOS; honors EMPYREAN_DATA_DIR) under the filenames the engine expects.

Downloaded by the engine when needed

Any of these can be relocated by pointing EMPYREAN_DATA_DIR at a directory holding them.

Provisioning up front

empyrean.download_data() provisions a usable data directory and does nothing else — it does not build a context or load anything. It is idempotent: files already present are kept, only what is missing is fetched, and whatever the installed B612 wheels supply is staged from disk with zero network access. It returns the directory it provisioned.

path = empyrean.download_data()            # or download_data(data_dir=...)

Running offline

empyrean.initialize(refresh=False) resolves the kernel set from the data directory alone — no downloads, no staleness checks. If anything is missing it raises FileNotFoundError naming every absent file, both in the message and as a missing_data_files list on the exception, so one run tells you exactly what to stage. There is no partial load and no quiet fall back to a smaller force model.

try:
    empyrean.initialize(refresh=False)
except FileNotFoundError as e:
    print("stage these first:", e.missing_data_files)

Setting EMPYREAN_OFFLINE=1 in the environment applies the same policy to every initialize() in the process and announces itself on stderr. It is a floor, not a switch: it can turn network access off, never back on. download_data() has no offline form — reaching the network is the whole call — so under the floor it refuses rather than provisioning anyway.

A fetch that was attempted and failed is missing data, not an I/O fault: the message leads with "Data download failed: " and names the failing kernel by URL. The two shapes are distinguishable — a populated missing_data_files is a named set the directory does not have with no fetch attempted, while the download prefix is an acquisition that was tried and did not land. The remedy for the second is connectivity, or a pin the upstream still serves.

Accuracy

Validated against JPL Horizons, find_orb, and GRSS across a curated catalog of 50 objects in 13 dynamical populations (NEOs, MBAs, Trojans, TNOs, comets, etc.), with ASSIST as an additional propagation reference on a 39-object subset. Sub-meter propagation accuracy on bounded timescales. See the validation notes.

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 shipped with this package for the full terms.

License

empyrean is dual-licensed:

• Wrapper / binding source code — the Rust API surface, C-ABI bindings, and Python wrapper sources in the main repository — is licensed under the BSD 3-Clause License.
• This Python wheel (and any other pre-compiled binary distribution of empyrean) is licensed under the proprietary Empyrean Binary License. The wheel is free to install and use (including commercial use) but may not be redistributed, modified, reverse-engineered, decompiled, or disassembled.

The BSD-3 grant covers only the binding / integration layers in the public repository. The propagation engine, orbit- determination engine, and automatic-differentiation library are proprietary closed-source components distributed only inside the compiled wheel — the wrapper sources call into them through stable internal APIs but do not contain their implementations. Cloning the repository will not let you build a working empyrean from source; install the published wheel.

Copyright © 2024–2026 Joachim Moeyens. All rights reserved.

Links

• Website: https://www.empyrean-dynamics.com
• Repository: https://github.com/Empyrean-Dynamics/empyrean
• Issues: https://github.com/Empyrean-Dynamics/empyrean/issues