kshana 0.35.0

// SPDX-License-Identifier: AGPL-3.0-only
//! A pluggable **ephemeris provider** seam — the source of one body's position relative to another
//! in the inertial frame — abstracted so the *fidelity* of those positions can be swapped without
//! touching the deep-space light-time correction (D0.7) or the interplanetary orbit determination
//! (D2/D3) that will consume it.
//!
//! ## The split, and why it is shaped this way
//!
//! This module ships the trait ([`EphemerisProvider`]), the **kernel-free, low-precision**
//! built-in implementation ([`BuiltinEphemeris`]), which reuses Kshana's own analytic Montenbruck &
//! Gill Sun/Moon series ([`crate::ephem`]), the kernel-free whole-solar-system
//! [`AnalyticSolarSystem`], and the **kernel** implementation [`KernelEphemeris`]: JPL
//! Development Ephemeris positions (for example DE440 from `de440s.bsp`) read by the engine's own
//! NAIF (Navigation and Ancillary Information Facility) kernel reader, [`crate::naif_kernel`].
//! That reader is pure Rust written from the public NAIF Required Reading documents and is
//! validated on its own matrix row (`tests/naif_reader_spice_oracle.rs`: 600 states and 200
//! rotations on a random grid drawn after pre-registration, against the SPICE Toolkit and ANISE,
//! every one inside a bar derived from double-precision Chebyshev evaluation). That row says the
//! engine reads the kernel correctly; it says nothing about how accurate DE440 is.
//!
//! The analytic series stay the default everywhere: a kernel is used only where a caller passes a
//! [`KernelEphemeris`] (in the lunar laser-ranging datum and the lunar VLBI, very long baseline
//! interferometry, information-matrix scenarios, through their `planetary_kernel_path` keys), and
//! a kernel run records the kernel's SHA-256 in its output so the two kinds of run can never be
//! mistaken for one another.
//!
//! ANISE (the pure-Rust NAIF/SPICE reimplementation) is still kept out of this crate, and the
//! cross-validation crates that use it (`xval/anise-mars-od`, `xval/anise-lunar-od`) stay
//! workspace-excluded, for a hard CI (continuous integration) reason: ANISE and `hifitime` are
//! licensed **MPL-2.0** (Mozilla Public License 2.0) and built on **edition 2024 (Rust ≥ 1.85)**.
//! Pulling either into the main `kshana` dependency graph — even behind a cargo feature — would
//! break two otherwise-green CI gates:
//!
//! 1. `cargo deny check licenses` — MPL-2.0 sits outside the crate's strict license allow-list; and
//! 2. the **MSRV 1.75** (minimum supported Rust version) build — cargo 1.75 cannot even parse an
//!    edition-2024 manifest in the resolved dependency graph.
//!
//! [`KernelEphemeris`] needs neither: it reads the kernel with [`crate::naif_kernel`]. See
//! `xval/anise-lunar-od/Cargo.toml`'s header for the rationale applied to the lunar environment.

use crate::body::Body;
use crate::ephem::{
    ecliptic_to_icrf, moon_icrf, moon_icrf_velocity, moon_position, satellite_state,
    standish_state, sun_position, Planet, Satellite, StandishTable,
};

use crate::jd2::Jd2;

type Vec3 = [f64; 3];

/// Seconds per day, for the ephemeris-time conversion of the kernel path.
const SECONDS_PER_DAY: f64 = 86_400.0;

/// J2000.0 epoch as a Julian Date (TT), the zero of `t_tt_jc`.
const JD_J2000: f64 = 2_451_545.0;
/// Julian days per Julian century, the scale of `t_tt_jc`.
const DAYS_PER_JULIAN_CENTURY: f64 = 36_525.0;

/// A pluggable source of one body's position relative to another in the inertial frame.
///
/// Implementors return the position of `target` **relative to** `center` in the ICRF/EME2000
/// inertial frame (metres) at the requested epoch (Julian Date, **TDB**). A provider returns `None`
/// for any `(target, center)` pair it has no series or kernel for — letting a caller fall back to a
/// higher-fidelity provider ([`KernelEphemeris`], or the DE-grade `xval/anise-mars-od` path) for
/// the bodies the kernel-free [`BuiltinEphemeris`] cannot supply.
pub trait EphemerisProvider: std::fmt::Debug {
    /// Position of `target` relative to `center` in the ICRF/EME2000 inertial frame (metres) at the
    /// epoch `jd_tdb` (Julian Date, TDB). Returns `None` for a `(target, center)` pair this provider
    /// cannot supply.
    fn relative_position(&self, target: &Body, center: &Body, jd_tdb: f64) -> Option<Vec3>;
}

/// The default, kernel-free **low-precision** ephemeris provider: Kshana's own analytic Montenbruck
/// & Gill Sun/Moon series ([`crate::ephem`]).
///
/// It knows the **geocentric** Sun and Moon directions only — the closed-form series give the Sun to
/// ~0.005° and the Moon to ~0.3° about J2000 (Montenbruck & Gill, *Satellite Orbits* §3.3.2) — so it
/// can supply the Earth/Sun/Moon mutual positions and nothing else. For **Mars or any other planet**
/// it has no series and returns `None`; [`KernelEphemeris`] is the path for those, as
/// [`crate::lunar_od::AnalyticLunarEnvironment`] gives way to the out-of-crate
/// `AniseLunarEnvironment` for DE-grade lunar inputs.
///
/// The analytic series are parameterised by `t_tt_jc`, Julian centuries of **TT** since J2000. This
/// provider takes the epoch as TDB and uses TDB ≈ TT directly: the TDB−TT difference is a periodic
/// term bounded by ~1.7 ms, which moves the Sun/Moon directions by far less than the series' own
/// ~0.005°/~0.3° truncation error, so it is dropped here (the DE-grade provider carries the exact
/// TDB timescale).
///
/// Positions are returned in the geocentric mean-equator/equinox of date — a close approximation to
/// the ICRF/EME2000 frame whose precession/nutation difference is well below this model's own
/// truncation error (the same approximation [`crate::ephem`] documents).
#[derive(Debug, Clone, Copy, Default)]
pub struct BuiltinEphemeris;

impl BuiltinEphemeris {
    /// Convert a Julian Date (TDB) to `t_tt_jc`, Julian centuries of TT since J2000 — the argument
    /// of the analytic [`crate::ephem`] series. TDB ≈ TT is used (see the type-level note).
    fn t_tt_jc(jd_tdb: f64) -> f64 {
        (jd_tdb - JD_J2000) / DAYS_PER_JULIAN_CENTURY
    }
}

impl EphemerisProvider for BuiltinEphemeris {
    /// Supply the geocentric Earth/Sun/Moon mutual positions the built-in analytic series know:
    ///
    /// * any body relative to **itself** → `[0, 0, 0]` (matched on [`Body::name`]);
    /// * **Sun** relative to **Earth** → [`crate::ephem::sun_position`];
    /// * **Moon** relative to **Earth** → [`crate::ephem::moon_position`];
    /// * the reverse pairs (Earth relative to the Sun / Moon) → the negated forward vector.
    ///
    /// Anything involving **Mars** or another planet, and any pair the geocentric series cannot
    /// compose (e.g. Sun relative to Moon), returns `None`. Dispatch is on [`Body::name`].
    fn relative_position(&self, target: &Body, center: &Body, jd_tdb: f64) -> Option<Vec3> {
        // A body relative to itself is the origin, for any body the series knows or not.
        if target.name == center.name {
            return Some([0.0, 0.0, 0.0]);
        }
        let t = Self::t_tt_jc(jd_tdb);
        match (target.name, center.name) {
            ("Sun", "Earth") => Some(sun_position(t)),
            ("Moon", "Earth") => Some(moon_position(t)),
            ("Earth", "Sun") => Some(neg(sun_position(t))),
            ("Earth", "Moon") => Some(neg(moon_position(t))),
            // Mars and every other planet have no built-in series; geocentric Sun↔Moon pairs
            // are not composed here. `KernelEphemeris` is the path.
            _ => None,
        }
    }
}

/// A kernel-free **whole-solar-system** ephemeris: any body of [`crate::body::SOLAR_SYSTEM`]
/// relative to any other, in the ICRF (equatorial J2000), metres.
///
/// Built from published analytic sources only, composed through heliocentric positions:
///
/// * the planets (and the Earth-Moon barycentre) from the JPL Standish Keplerian elements
///   ([`crate::ephem::standish_state`]), rotated from the J2000 ecliptic to the ICRF;
/// * the Earth and the Moon from the barycentre and the geocentric Moon series
///   ([`crate::ephem::moon_icrf`]), split by the mass ratio `μ_Moon / (μ_Earth + μ_Moon)`;
/// * Phobos, Deimos, the Galilean moons and Titan from [`crate::ephem::satellite_state`],
///   added to their planet.
///
/// The Standish elements give a planetary *system*, which is used for the planet itself: the
/// offset of a planet from its system barycentre (about 200 km for Jupiter) is far below the
/// Standish error at that planet. The Sun is the origin: the Standish positions are heliocentric.
/// The epoch is TDB; the geocentric Moon series takes TT, within 2 ms of it.
#[derive(Debug, Clone, Copy, Default)]
pub struct AnalyticSolarSystem {
    /// The Standish table to use, or `None` for the table whose fit interval contains the epoch
    /// ([`StandishTable::for_epoch`]).
    pub table: Option<StandishTable>,
}

impl AnalyticSolarSystem {
    /// The Standish table used at `jd_tdb`: the forced one, else the one whose fit interval
    /// contains the epoch. `None` outside 3000 BC to 3000 AD.
    pub fn table_at(&self, jd_tdb: f64) -> Option<StandishTable> {
        match self.table {
            Some(t) => Some(t),
            None => StandishTable::for_epoch(Self::t_tdb_jc(jd_tdb)),
        }
    }

    fn t_tdb_jc(jd_tdb: f64) -> f64 {
        (jd_tdb - JD_J2000) / DAYS_PER_JULIAN_CENTURY
    }

    fn planet_icrf(&self, p: Planet, jd_tdb: f64) -> Option<(Vec3, Vec3)> {
        let table = self.table_at(jd_tdb)?;
        let s = standish_state(p, Self::t_tdb_jc(jd_tdb), table)?;
        Some((ecliptic_to_icrf(s.pos_m), ecliptic_to_icrf(s.vel_m_s)))
    }

    /// Heliocentric ICRF position (m) and velocity (m/s) of the body called `name` (as in
    /// [`Body::name`]) at `jd_tdb`. `None` for an unknown name, outside the Standish tables'
    /// interval, or for Pluto outside 1800 AD to 2050 AD.
    pub fn heliocentric_state(&self, name: &str, jd_tdb: f64) -> Option<(Vec3, Vec3)> {
        let planet = |p| self.planet_icrf(p, jd_tdb);
        let with_moon = |sat: Satellite| -> Option<(Vec3, Vec3)> {
            let (pp, pv) = self.planet_icrf(sat.parent(), jd_tdb)?;
            let s = satellite_state(sat, jd_tdb);
            Some((add(pp, s.pos_m), add(pv, s.vel_m_s)))
        };
        match name {
            "Sun" => Some(([0.0; 3], [0.0; 3])),
            "Mercury" => planet(Planet::Mercury),
            "Venus" => planet(Planet::Venus),
            "Earth" | "Moon" => {
                let (bp, bv) = planet(Planet::EarthMoonBarycentre)?;
                let mp = moon_icrf(jd_tdb);
                let mv = moon_icrf_velocity(jd_tdb);
                let mu_e = crate::forces::MU_EARTH;
                let mu_m = crate::forces::MU_MOON;
                // Earth sits μ_M/(μ_E+μ_M) of the way from the barycentre towards minus the Moon.
                let k = if name == "Earth" {
                    -mu_m / (mu_e + mu_m)
                } else {
                    mu_e / (mu_e + mu_m)
                };
                Some((add(bp, scale(mp, k)), add(bv, scale(mv, k))))
            }
            "Mars" => planet(Planet::Mars),
            "Jupiter" => planet(Planet::Jupiter),
            "Saturn" => planet(Planet::Saturn),
            "Uranus" => planet(Planet::Uranus),
            "Neptune" => planet(Planet::Neptune),
            "Pluto" => planet(Planet::Pluto),
            "Phobos" => with_moon(Satellite::Phobos),
            "Deimos" => with_moon(Satellite::Deimos),
            "Io" => with_moon(Satellite::Io),
            "Europa" => with_moon(Satellite::Europa),
            "Ganymede" => with_moon(Satellite::Ganymede),
            "Callisto" => with_moon(Satellite::Callisto),
            "Titan" => with_moon(Satellite::Titan),
            _ => None,
        }
    }

    /// One line naming the model behind the body called `name` at `jd_tdb`.
    pub fn method(&self, name: &str, jd_tdb: f64) -> String {
        let table = self
            .table_at(jd_tdb)
            .map(|t| t.label())
            .unwrap_or("outside every Standish table");
        match name {
            "Sun" => "origin of the heliocentric frame".to_string(),
            "Earth" | "Moon" => format!(
                "Earth-Moon barycentre from {table}, split by the Montenbruck & Gill lunar series \
                 precessed to J2000"
            ),
            "Phobos" | "Deimos" | "Io" | "Europa" | "Ganymede" | "Callisto" => format!(
                "JPL mean elements in the Laplace plane, mean-longitude rate from the IAU \
                 synchronous rotation rate, added to its planet from {table}"
            ),
            "Titan" => format!(
                "IAU synchronous-rotation direction at the JPL mean distance, added to Saturn \
                 from {table}"
            ),
            "Pluto" => "the 1992 Standish Table 1 Pluto row (1800 AD to 2050 AD)".to_string(),
            _ => table.to_string(),
        }
    }
}

impl EphemerisProvider for AnalyticSolarSystem {
    /// `target` relative to `center` in the ICRF (m): the difference of their heliocentric
    /// positions. `None` when either body is unknown or outside the model's interval.
    fn relative_position(&self, target: &Body, center: &Body, jd_tdb: f64) -> Option<Vec3> {
        if target.name == center.name {
            return Some([0.0, 0.0, 0.0]);
        }
        let (t, _) = self.heliocentric_state(target.name, jd_tdb)?;
        let (c, _) = self.heliocentric_state(center.name, jd_tdb)?;
        Some([t[0] - c[0], t[1] - c[1], t[2] - c[2]])
    }
}

/// A **kernel** ephemeris: JPL Development Ephemeris states read from an SPK (Spacecraft and
/// Planet Kernel) file by the engine's own reader ([`crate::naif_kernel::SpkKernel`]).
///
/// Positions are geometric (no light time or aberration), in the J2000 frame of the kernel, in
/// metres; epochs are **TDB** (barycentric dynamical time, which is what SPICE calls ephemeris
/// time) given as a two-part Julian date ([`Jd2`]) so a sub-microsecond epoch survives: one
/// `f64` Julian date near J2000 resolves only about 40 µs, about 4 cm of the Moon's motion.
///
/// Bodies are addressed by NAIF integer code ([`KernelEphemeris::state`]) or, through
/// [`EphemerisProvider`], by [`Body::name`] mapped to its code in [`crate::body::SOLAR_SYSTEM`].
/// A body whose own code is not in the kernel gives `None` (or an error from the coded calls):
/// a planet is never silently replaced by its system barycentre (Jupiter's is about 200 km from
/// Jupiter), so with `de440s.bsp` the providers here supply the Sun, Mercury, Venus, the Earth
/// and the Moon, and the barycentres 0 to 9 by code.
///
/// What its correctness rests on: the reader row (`tests/naif_reader_spice_oracle.rs`), which
/// validates reading, not DE440. The kernel file's SHA-256 is kept so a run can record exactly
/// which kernel produced it ([`KernelEphemeris::kernel_sha256`]).
#[derive(Debug, Clone)]
pub struct KernelEphemeris {
    spk: crate::naif_kernel::SpkKernel,
    sha256: String,
    source: String,
}

/// NAIF codes of the bodies the lunar paths ask a kernel for.
pub mod naif_code {
    /// The Sun.
    pub const SUN: i32 = 10;
    /// The Earth-Moon barycentre.
    pub const EARTH_MOON_BARYCENTRE: i32 = 3;
    /// The Moon.
    pub const MOON: i32 = 301;
    /// The Earth.
    pub const EARTH: i32 = 399;
}

impl KernelEphemeris {
    /// Read an SPK kernel from disk and record its SHA-256.
    pub fn open(path: &std::path::Path) -> Result<Self, String> {
        let bytes = std::fs::read(path).map_err(|e| format!("read {}: {e}", path.display()))?;
        Self::from_bytes(bytes, path.display().to_string())
    }

    /// Wrap the bytes of an SPK kernel; `source` names it in reports (a path, usually).
    pub fn from_bytes(bytes: Vec<u8>, source: String) -> Result<Self, String> {
        use sha2::{Digest, Sha256};
        let sha256 = Sha256::digest(&bytes)
            .iter()
            .map(|b| format!("{b:02x}"))
            .collect();
        let spk = crate::naif_kernel::SpkKernel::from_daf(
            crate::naif_kernel::DafFile::from_bytes(bytes)?,
        )?;
        Ok(Self {
            spk,
            sha256,
            source,
        })
    }

    /// SHA-256 (lower-case hexadecimal) of the kernel file this ephemeris reads.
    pub fn kernel_sha256(&self) -> &str {
        &self.sha256
    }

    /// Where the kernel came from, as given to [`open`](Self::open).
    pub fn source(&self) -> &str {
        &self.source
    }

    /// Ephemeris time (TDB seconds past J2000) of a TDB two-part Julian date, as the
    /// `(hi, lo)` pair the reader takes. The whole-day part is exact, so `hi` carries no rounding.
    pub fn et_of(jd_tdb: Jd2) -> (f64, f64) {
        (
            (jd_tdb.day - JD_J2000) * SECONDS_PER_DAY,
            jd_tdb.frac * SECONDS_PER_DAY,
        )
    }

    /// Geometric state (position m, velocity m/s, acceleration m/s²; J2000) of NAIF body
    /// `target` relative to NAIF body `center` at the TDB epoch `jd_tdb`.
    pub fn state(
        &self,
        target: i32,
        center: i32,
        jd_tdb: Jd2,
    ) -> Result<crate::naif_kernel::State3, String> {
        let (hi, lo) = Self::et_of(jd_tdb);
        self.spk.state(target, center, hi, lo)
    }

    /// Geometric position (m, J2000) of NAIF body `target` relative to `center` at TDB `jd_tdb`.
    pub fn position(&self, target: i32, center: i32, jd_tdb: Jd2) -> Result<Vec3, String> {
        Ok(self.state(target, center, jd_tdb)?[0])
    }

    /// The geocentric Moon (m, J2000) at TDB `jd_tdb`: the kernel counterpart of
    /// [`crate::ephem::moon_position`].
    pub fn moon_geocentric(&self, jd_tdb: Jd2) -> Result<Vec3, String> {
        self.position(naif_code::MOON, naif_code::EARTH, jd_tdb)
    }

    /// The geocentric Sun (m, J2000) at TDB `jd_tdb`: the kernel counterpart of
    /// [`crate::ephem::sun_position`].
    pub fn sun_geocentric(&self, jd_tdb: Jd2) -> Result<Vec3, String> {
        self.position(naif_code::SUN, naif_code::EARTH, jd_tdb)
    }

    /// The geocentric Moon (m, J2000) at a **TT** (terrestrial time) two-part Julian date, the
    /// time argument the analytic series take: TT is carried to TDB with
    /// [`crate::timescales::tdb_minus_tt`] (at most about 1.7 ms, which the series ignore but
    /// which is about 1.7 m of the Moon's motion) before the kernel is read.
    pub fn moon_geocentric_tt(&self, jd_tt: Jd2) -> Result<Vec3, String> {
        self.moon_geocentric(tt_to_tdb_jd2(jd_tt))
    }

    /// The geocentric Sun (m, J2000) at a TT two-part Julian date (see
    /// [`moon_geocentric_tt`](Self::moon_geocentric_tt)).
    pub fn sun_geocentric_tt(&self, jd_tt: Jd2) -> Result<Vec3, String> {
        self.sun_geocentric(tt_to_tdb_jd2(jd_tt))
    }

    /// Position of `target` relative to `center` (m, J2000) at a **UTC** two-part epoch, carried
    /// to TDB by [`utc_to_tdb_jd2`]; `None` as for [`relative_position_jd2`](Self::relative_position_jd2).
    pub fn relative_position_utc(&self, target: &Body, center: &Body, jd_utc: Jd2) -> Option<Vec3> {
        self.relative_position_jd2(target, center, utc_to_tdb_jd2(jd_utc))
    }

    /// The NAIF code of a body named as [`Body::name`] names it.
    pub fn naif_code_of(name: &str) -> Option<i32> {
        crate::body::SOLAR_SYSTEM
            .iter()
            .find(|f| f.name == name)
            .map(|f| f.naif_id)
    }

    /// Position of `target` relative to `center` (m, J2000) at the TDB two-part epoch, or `None`
    /// when either body has no NAIF code or no data in the kernel at that epoch.
    pub fn relative_position_jd2(&self, target: &Body, center: &Body, jd_tdb: Jd2) -> Option<Vec3> {
        if target.name == center.name {
            return Some([0.0, 0.0, 0.0]);
        }
        let t = Self::naif_code_of(target.name)?;
        let c = Self::naif_code_of(center.name)?;
        self.position(t, c, jd_tdb).ok()
    }
}

/// UTC two-part Julian date to TDB: TT = UTC + (TAI − UTC) + 32.184 s, with the leap seconds of
/// [`crate::timescales::tai_minus_utc`], then [`tt_to_tdb_jd2`]. The two-term TDB − TT series
/// is good to a few tens of microseconds; leap seconds are exact from 1972 on.
pub fn utc_to_tdb_jd2(jd_utc: Jd2) -> Jd2 {
    let tt_minus_utc = crate::timescales::tai_minus_utc(jd_utc.total()) + 32.184;
    tt_to_tdb_jd2(jd_utc.add_seconds(tt_minus_utc))
}

/// TT two-part Julian date to TDB, with [`crate::timescales::tdb_minus_tt`].
pub fn tt_to_tdb_jd2(jd_tt: Jd2) -> Jd2 {
    jd_tt.add_seconds(crate::timescales::tdb_minus_tt(jd_tt.total()))
}

impl EphemerisProvider for KernelEphemeris {
    /// The kernel position of `target` relative to `center` (J2000, m). The single-`f64` TDB
    /// epoch is split with [`Jd2::new`]; callers that hold a two-part epoch should use
    /// [`KernelEphemeris::relative_position_jd2`] and keep its precision.
    fn relative_position(&self, target: &Body, center: &Body, jd_tdb: f64) -> Option<Vec3> {
        self.relative_position_jd2(target, center, Jd2::new(jd_tdb))
    }
}

/// Where a lunar or solar geometry path takes the geocentric Moon and Sun from.
///
/// `Analytic` is the default everywhere and is the engine's behaviour before kernels existed,
/// bit for bit: callers keep their own analytic expressions on that branch. `Kernel` reads the
/// given [`KernelEphemeris`].
#[derive(Debug, Clone, Copy, Default)]
pub enum LunisolarSource<'a> {
    /// The analytic Montenbruck & Gill series ([`crate::ephem`]).
    #[default]
    Analytic,
    /// A JPL kernel read by the engine's own reader.
    Kernel(&'a KernelEphemeris),
}

impl LunisolarSource<'_> {
    /// The kernel, when this source is one.
    pub fn kernel(&self) -> Option<&KernelEphemeris> {
        match self {
            LunisolarSource::Analytic => None,
            LunisolarSource::Kernel(k) => Some(k),
        }
    }

    /// `"analytic"` or `"kernel"`, for reports.
    pub fn label(&self) -> &'static str {
        match self {
            LunisolarSource::Analytic => "analytic",
            LunisolarSource::Kernel(_) => "kernel",
        }
    }
}

fn add(a: Vec3, b: Vec3) -> Vec3 {
    [a[0] + b[0], a[1] + b[1], a[2] + b[2]]
}

fn scale(a: Vec3, k: f64) -> Vec3 {
    [a[0] * k, a[1] * k, a[2] * k]
}

/// Negate a 3-vector (the reverse relative-position direction).
fn neg(v: Vec3) -> Vec3 {
    [-v[0], -v[1], -v[2]]
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::ephem::{moon_position, sun_position};

    /// A fixed probe epoch (Julian Date, TDB) inside the analytic series' good span: 2022-001.5.
    const PROBE_JD: f64 = 2_459_580.5;

    fn t_at(jd: f64) -> f64 {
        (jd - 2_451_545.0) / 36_525.0
    }

    /// Sun relative to Earth must be exactly the built-in geocentric `sun_position` at the same
    /// epoch — the seam adds no transform of its own beyond the documented TDB ≈ TT identity.
    #[test]
    fn sun_relative_to_earth_is_the_builtin_sun_series() {
        let p = BuiltinEphemeris;
        let got = p
            .relative_position(&Body::sun(), &Body::earth(), PROBE_JD)
            .expect("builtin must supply Sun relative to Earth");
        let want = sun_position(t_at(PROBE_JD));
        for k in 0..3 {
            assert!(
                (got[k] - want[k]).abs() < 1e-9,
                "Sun/Earth component {k}: got {} want {}",
                got[k],
                want[k]
            );
        }
    }

    /// Moon relative to Earth must be exactly the built-in geocentric `moon_position`.
    #[test]
    fn moon_relative_to_earth_is_the_builtin_moon_series() {
        let p = BuiltinEphemeris;
        let got = p
            .relative_position(&Body::moon(), &Body::earth(), PROBE_JD)
            .expect("builtin must supply Moon relative to Earth");
        let want = moon_position(t_at(PROBE_JD));
        for k in 0..3 {
            assert!(
                (got[k] - want[k]).abs() < 1e-9,
                "Moon/Earth component {k}: got {} want {}",
                got[k],
                want[k]
            );
        }
    }

    /// A body relative to itself is the origin.
    #[test]
    fn body_relative_to_itself_is_the_origin() {
        let p = BuiltinEphemeris;
        for b in [Body::earth(), Body::sun(), Body::moon(), Body::mars()] {
            let got = p
                .relative_position(&b, &b, PROBE_JD)
                .expect("a body relative to itself is always defined");
            assert_eq!(got, [0.0, 0.0, 0.0], "{} relative to itself", b.name);
        }
    }

    /// The reverse pair: Earth relative to the Sun is the negated Sun-relative-to-Earth vector.
    #[test]
    fn earth_relative_to_sun_is_the_negated_sun_vector() {
        let p = BuiltinEphemeris;
        let earth_sun = p
            .relative_position(&Body::earth(), &Body::sun(), PROBE_JD)
            .expect("builtin must supply Earth relative to Sun");
        let sun_earth = sun_position(t_at(PROBE_JD));
        for k in 0..3 {
            assert!(
                (earth_sun[k] + sun_earth[k]).abs() < 1e-9,
                "Earth/Sun component {k} must be the negated Sun/Earth: {} vs {}",
                earth_sun[k],
                -sun_earth[k]
            );
        }
    }

    /// The reverse pair for the Moon: Earth relative to the Moon is the negated Moon vector.
    #[test]
    fn earth_relative_to_moon_is_the_negated_moon_vector() {
        let p = BuiltinEphemeris;
        let earth_moon = p
            .relative_position(&Body::earth(), &Body::moon(), PROBE_JD)
            .expect("builtin must supply Earth relative to Moon");
        let moon_earth = moon_position(t_at(PROBE_JD));
        for k in 0..3 {
            assert!(
                (earth_moon[k] + moon_earth[k]).abs() < 1e-9,
                "Earth/Moon component {k} must be the negated Moon/Earth",
            );
        }
    }

    /// The built-in knows no Mars (or any planet) series, so any Mars pair returns `None` — the
    /// signal a caller uses to reach for a kernel provider ([`KernelEphemeris`]).
    #[test]
    fn mars_relative_to_earth_is_unsupported() {
        let p = BuiltinEphemeris;
        assert!(
            p.relative_position(&Body::mars(), &Body::earth(), PROBE_JD)
                .is_none(),
            "the kernel-free builtin must not invent a Mars position"
        );
        assert!(
            p.relative_position(&Body::earth(), &Body::mars(), PROBE_JD)
                .is_none(),
            "the reverse Mars pair is equally unsupported"
        );
        assert!(
            p.relative_position(&Body::sun(), &Body::mars(), PROBE_JD)
                .is_none(),
            "Sun relative to Mars is unsupported by the geocentric builtin"
        );
    }

    /// A geocentric series cannot compose Sun-relative-to-Moon, so that pair is `None` (it is not
    /// a self-pair, and the builtin only anchors on Earth).
    #[test]
    fn uncomposable_sun_relative_to_moon_is_none() {
        let p = BuiltinEphemeris;
        assert!(
            p.relative_position(&Body::sun(), &Body::moon(), PROBE_JD)
                .is_none(),
            "the builtin does not compose Sun relative to Moon"
        );
    }

    /// `BuiltinEphemeris` is usable behind the trait object the consumers (D0.7/D2/D3) will hold.
    #[test]
    fn usable_as_a_trait_object() {
        let p: &dyn EphemerisProvider = &BuiltinEphemeris;
        assert!(p
            .relative_position(&Body::sun(), &Body::earth(), PROBE_JD)
            .is_some());
    }

    #[test]
    fn analytic_solar_system_is_antisymmetric_and_rejects_unknown_bodies() {
        let e = AnalyticSolarSystem::default();
        let a = e
            .relative_position(&Body::jupiter(), &Body::mars(), PROBE_JD)
            .unwrap();
        let b = e
            .relative_position(&Body::mars(), &Body::jupiter(), PROBE_JD)
            .unwrap();
        for k in 0..3 {
            assert_eq!(a[k], -b[k]);
        }
        assert!(e.heliocentric_state("Vulcan", PROBE_JD).is_none());
        // Outside 3000 BC to 3000 AD there is no table and no position.
        assert!(e
            .heliocentric_state("Mars", 2_451_545.0 + 20.0 * 36_525.0)
            .is_none());
    }

    #[test]
    fn earth_and_moon_balance_about_their_barycentre() {
        // μ_E·r_E + μ_M·r_M = (μ_E + μ_M)·r_EMB: the split must conserve the barycentre.
        let e = AnalyticSolarSystem::default();
        let (re, _) = e.heliocentric_state("Earth", PROBE_JD).unwrap();
        let (rm, _) = e.heliocentric_state("Moon", PROBE_JD).unwrap();
        let t = (PROBE_JD - 2_451_545.0) / 36_525.0;
        let bary = crate::ephem::ecliptic_to_icrf(
            standish_state(Planet::EarthMoonBarycentre, t, StandishTable::Table1)
                .unwrap()
                .pos_m,
        );
        let (me, mm) = (crate::forces::MU_EARTH, crate::forces::MU_MOON);
        for k in 0..3 {
            let w = (me * re[k] + mm * rm[k]) / (me + mm);
            assert!((w - bary[k]).abs() < 1e-3, "component {k}");
        }
        // And the Earth-Moon distance is the lunar series' own.
        let d = [rm[0] - re[0], rm[1] - re[1], rm[2] - re[2]];
        let r = (d[0] * d[0] + d[1] * d[1] + d[2] * d[2]).sqrt();
        assert!((3.5e8..4.1e8).contains(&r));
    }

    /// The cut DE440 kernel committed for the lunar VLBI comparison (2023-12-31..2024-01-02;
    /// segments 3 wrt 0, 301 wrt 3, 399 wrt 3), records bit-identical to NAIF's `de440s.bsp`.
    fn cut_kernel() -> KernelEphemeris {
        KernelEphemeris::open(std::path::Path::new(concat!(
            env!("CARGO_MANIFEST_DIR"),
            "/tests/fixtures/lunar_vlbi_anise_oracle/kernels/de440s_2024-01-01.bsp"
        )))
        .expect("the committed cut kernel opens")
    }

    #[test]
    fn kernel_ephemeris_records_the_kernel_digest() {
        // The digest the fixture's NOTICE.md pins for this file.
        // PIN-SCOPE:    the SHA-256 of the committed cut kernel, equal to the digest its NOTICE.md records.
        // PIN-EXCLUDES: the positions read from the kernel, compared with their oracle elsewhere.
        assert_eq!(
            cut_kernel().kernel_sha256(),
            "8f7986fcc8e2987c9d94efa86e02a2a578631b2cb939d497a395f1c58c3b739e"
        );
    }

    #[test]
    fn kernel_epoch_split_is_exact_on_whole_days() {
        assert_eq!(KernelEphemeris::et_of(Jd2::new(2_451_545.0)), (0.0, 0.0));
        let (hi, lo) = KernelEphemeris::et_of(Jd2::from_parts(2_460_311.0, 0.25));
        assert_eq!(hi, 8_766.0 * 86_400.0);
        assert_eq!(lo, 21_600.0);
    }

    #[test]
    fn kernel_moon_is_the_reader_state_of_301_relative_to_399() {
        let k = cut_kernel();
        let jd = Jd2::from_parts(2_460_311.0, 0.125);
        let (hi, lo) = KernelEphemeris::et_of(jd);
        let spk = crate::naif_kernel::SpkKernel::open(std::path::Path::new(concat!(
            env!("CARGO_MANIFEST_DIR"),
            "/tests/fixtures/lunar_vlbi_anise_oracle/kernels/de440s_2024-01-01.bsp"
        )))
        .unwrap();
        assert_eq!(
            k.moon_geocentric(jd).unwrap(),
            spk.state(301, 399, hi, lo).unwrap()[0]
        );
        // Through the provider trait, by body name, the same vector.
        let via_trait = k
            .relative_position_jd2(&Body::moon(), &Body::earth(), jd)
            .unwrap();
        assert_eq!(via_trait, k.moon_geocentric(jd).unwrap());
        let back = k
            .relative_position_jd2(&Body::earth(), &Body::moon(), jd)
            .unwrap();
        for i in 0..3 {
            assert_eq!(back[i], -via_trait[i]);
        }
    }

    #[test]
    fn kernel_never_substitutes_a_barycentre_or_invents_data() {
        let k = cut_kernel();
        let jd = 2_460_311.0;
        // This cut holds no Mars (499) or Sun (10) data at all, and nothing outside its window.
        assert!(k
            .relative_position(&Body::mars(), &Body::earth(), jd)
            .is_none());
        assert!(k.sun_geocentric(Jd2::new(jd)).is_err());
        assert!(k.moon_geocentric(Jd2::new(2_451_545.0)).is_err());
    }

    #[test]
    fn kernel_and_analytic_moon_agree_within_the_series_stated_bound() {
        // The analytic series states ~0.3 deg; the kernel is DE440. Same epoch, TT to TDB.
        let k = cut_kernel();
        let jd_tt = Jd2::from_parts(2_460_311.0, 0.0);
        let kern = k.moon_geocentric_tt(jd_tt).unwrap();
        let ana = moon_position((jd_tt.total() - 2_451_545.0) / 36_525.0);
        let dot: f64 = (0..3).map(|i| kern[i] * ana[i]).sum();
        let n = |v: Vec3| (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]).sqrt();
        let ang = (dot / (n(kern) * n(ana)))
            .clamp(-1.0, 1.0)
            .acos()
            .to_degrees();
        assert!(ang < 0.3, "{ang} deg");
    }

    #[test]
    fn tt_to_tdb_adds_the_periodic_offset() {
        let jd = Jd2::from_parts(2_460_311.0, 0.0);
        let d = tt_to_tdb_jd2(jd).diff_seconds(jd);
        let want = crate::timescales::tdb_minus_tt(jd.total());
        assert!((d - want).abs() < 1e-9, "{d} vs {want}");
        assert!(d.abs() < 2e-3);
    }
}