astrodyn_gravity 0.2.0

//! Per-source gravity-control configuration ([`GravityControl`] +
//! typed sibling [`GravityControlTyped`]).
//!
//! Ports
//! [`models/environment/gravity/src/gravity_controls.cc`](https://github.com/nasa/jeod/blob/jeod_v5.4.0/models/environment/gravity/src/gravity_controls.cc)
//! and
//! [`spherical_harmonics_gravity_controls.cc`](https://github.com/nasa/jeod/blob/jeod_v5.4.0/models/environment/gravity/src/spherical_harmonics_gravity_controls.cc)
//! from JEOD v5.4.0. A `GravityControl` selects whether a single source
//! contributes point-mass or spherical-harmonics gravity, picks the
//! degree / order, gates the gradient computation, and flags
//! third-body / Battin / relativistic corrections.

use astrodyn_dynamics::GravityAcceleration;
use astrodyn_quantities::aliases::HarmonicDegree;
use astrodyn_quantities::frame_descriptor::FrameUid;
use glam::DMat3;
use glam::DVec3;

use crate::gravity_source::{GravityModel, GravitySource};

/// Self-documenting selector for the gravity-gradient tensor flag at
/// the [`GravityControl`] constructor seam. Replaces the bare `bool`
/// the spherical / non-spherical constructors used to accept, so a
/// call site reads as
/// `GravityControl::new_spherical(earth, GravityGradient::Skip)`
/// rather than the bare-`bool` form, and the reader does not need to
/// remember which boolean polarity meant what.
///
/// This enum *only* gates the [`GravityControl::gradient`] field
/// (compute the gradient tensor in addition to the acceleration
/// vector). The third-body / direct discriminant lives on the
/// separate [`GravityControl::differential`] field, which is set by
/// [`GravityControl::new_third_body`] — *not* by this enum.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum GravityGradient {
    /// Compute the gravity-gradient tensor in addition to the
    /// acceleration vector. Required for gravity-torque interaction
    /// and for the gravity-gradient sensitivity studies in JEOD's
    /// SIM_dyncomp RUN_5B.
    Compute,
    /// Skip the gravity-gradient tensor; only compute the
    /// acceleration vector. The default for a vehicle that does not
    /// model gravity-torque dynamics.
    Skip,
}

impl GravityGradient {
    /// Project to the underlying `bool` storage form
    /// ([`GravityControl::gradient`]). `Compute` maps to `true`,
    /// `Skip` to `false`, matching the bare-bool calling convention
    /// this enum replaces.
    #[inline]
    pub const fn as_bool(self) -> bool {
        matches!(self, Self::Compute)
    }
}

/// Per-source gravity selector — point-mass vs. spherical harmonics,
/// degree / order, gradient flags, and third-body / Battin /
/// relativistic toggles.
///
/// The source is referenced by its inertial-frame [`FrameUid`] — the
/// same value identity in every host (issue #668 collapsed the old
/// host-generic `SourceId` parameter: the runner resolves uid → source
/// index and the Bevy adapter resolves uid → source entity, each at its
/// own boundary, so one control expression serves both backends and a
/// wrong reference fails loudly *by name* instead of silently indexing
/// the wrong source).
#[derive(Debug, Clone, PartialEq)]
pub struct GravityControl {
    /// The source's inertial-frame identity (e.g.
    /// `FrameUid::of::<PlanetInertial<Earth>>()`).
    pub source: FrameUid,
    /// Compute the gravity-gradient tensor in addition to acceleration.
    pub gradient: bool,
    /// If true, use only point-mass (spherical) gravity for this source,
    /// ignoring any spherical harmonics data. Matches JEOD's `spherical` flag
    /// on `GravityControls`.
    pub spherical: bool,
    /// Non-spherical degree to use. Must be <= source degree.
    /// Ignored when `spherical` is true.
    pub degree: usize,
    /// Non-spherical order to use. Must be <= degree and <= source order.
    /// Ignored when `spherical` is true.
    pub order: usize,
    /// If true, exclude point-mass (n=0,1) terms.
    pub perturbing_only: bool,
    /// Degree for gradient computation. Must be <= degree.
    pub gradient_degree: usize,
    /// Order for gradient computation. Must be <= order and <= gradient_degree.
    pub gradient_order: usize,
    /// If true, compute gravity as differential acceleration: the acceleration
    /// of the vehicle toward this source minus the acceleration of the
    /// integration frame origin toward this source. This is the correct
    /// treatment for third-body perturbations (e.g., Sun/Moon when integrating
    /// in an Earth-centered frame).
    ///
    /// Matches JEOD's `GravityIntegFrame::is_third_body` flag. In JEOD, this
    /// is set automatically based on whether the source's inertial frame is a
    /// progeny of the integration frame. Here it is set explicitly per control.
    // JEOD_INV: GV.14 — third-body vs direct gravity classification (set explicitly; JEOD derives from frame tree ancestry)
    pub differential: bool,
    /// If true, use Battin's method for improved numerical accuracy in
    /// third-body (differential) gravity computation. Only meaningful when
    /// `differential` is also true. Off by default in JEOD.
    ///
    /// Battin's method reformulates the differential acceleration to avoid
    /// catastrophic cancellation when the vehicle is close to the integration
    /// frame origin relative to the third-body source distance.
    ///
    /// JEOD ref: `gravity_controls.cc:317-331`.
    pub battin_method: bool,
    /// If true, apply post-Newtonian relativistic correction for this source.
    /// Requires source velocity in `GravitySourceEntry`. Only significant for
    /// Mercury-like orbits near massive bodies.
    pub relativistic: bool,
}

impl GravityControl {
    /// Create a spherical (point-mass only) gravity control.
    ///
    /// Uses only µ/r² acceleration. Any spherical harmonics data on the source
    /// is ignored. For gravity with J2+ harmonics, use [`new_nonspherical`](Self::new_nonspherical).
    ///
    /// `source` accepts `impl Into<FrameUid>` so any future identity
    /// sugar flows through this seam without per-callsite plumbing;
    /// today callers pass a `FrameUid` directly (e.g.
    /// `FrameUid::of::<PlanetInertial<Earth>>()`).
    pub fn new_spherical(source: impl Into<FrameUid>, gradient: GravityGradient) -> Self {
        Self {
            source: source.into(),
            gradient: gradient.as_bool(),
            spherical: true,
            degree: 0,
            order: 0,
            perturbing_only: false,
            gradient_degree: 0,
            gradient_order: 0,
            differential: false,
            battin_method: false,
            relativistic: false,
        }
    }

    /// Create a non-spherical (spherical harmonics) gravity control.
    ///
    /// Evaluates the source's spherical harmonics coefficients up to the given
    /// `degree` and `order`. The source must have a `SphericalHarmonics` model
    /// and the gravity source entry must provide a planet-fixed rotation matrix.
    ///
    /// See [`Self::new_spherical`] for the `impl Into<FrameUid>` rationale.
    pub fn new_nonspherical(
        source: impl Into<FrameUid>,
        degree: usize,
        order: usize,
        gradient: GravityGradient,
    ) -> Self {
        Self {
            source: source.into(),
            gradient: gradient.as_bool(),
            spherical: false,
            degree,
            order,
            perturbing_only: false,
            gradient_degree: 0,
            gradient_order: 0,
            differential: false,
            battin_method: false,
            relativistic: false,
        }
    }

    /// Create a spherical (point-mass) gravity control for a third-body source.
    ///
    /// Third-body sources use differential acceleration: the acceleration of
    /// the vehicle toward this source minus the acceleration of the integration
    /// frame origin toward this source.
    ///
    /// See [`Self::new_spherical`] for the `impl Into<FrameUid>` rationale.
    pub fn new_third_body(source: impl Into<FrameUid>) -> Self {
        Self {
            source: source.into(),
            gradient: false,
            spherical: true,
            degree: 0,
            order: 0,
            perturbing_only: false,
            gradient_degree: 0,
            gradient_order: 0,
            differential: true,
            battin_method: false,
            relativistic: false,
        }
    }

    /// Validate this control against its gravity source.
    ///
    /// Ported from JEOD's
    /// `SphericalHarmonicsGravityControls::check_validity()`, but with
    /// fail-loud semantics: where JEOD logs a non-fatal
    /// `MessageHandler::error` and silently auto-corrects the control,
    /// we panic. A physics simulation has no graceful-degradation mode
    /// (CLAUDE.md "Fail Loudly"); a mission crate that swaps planet
    /// sources and forgets to update the gravity degree, or whose
    /// gradient ordinals fall out of range, must surface the
    /// misconfiguration at construction rather than silently propagate
    /// with a different gravity model than the operator requested.
    ///
    /// # Panics
    /// - `spherical` is false and `degree < 2` (use
    ///   `GravityControl::new_spherical` for point-mass). `degree == 1`
    ///   is also rejected here because Gottlieb returns zero perturbation
    ///   for degree < 2 — accepting it would let `effective_orders`
    ///   silently collapse the request to point-mass and violate Fail
    ///   Loudly.
    /// - `spherical` is false against a `GravityModel::PointMass` source.
    /// - `degree > source.degree` or `order > source.order`.
    /// - `order > degree`.
    /// - `gradient` is true and any of:
    ///   `gradient_degree > degree`, `gradient_degree == 1`,
    ///   `gradient_order > gradient_degree`, `gradient_order > order`.
    // JEOD_INV: GV.03 — check_validity() called on degree/order mutation
    pub fn check_validity(&self, source: &GravitySource) {
        if self.spherical {
            return;
        }

        // JEOD_INV: GV.07 — degree < 2 with spherical=false panics
        // (JEOD `spherical_harmonics_gravity_controls.cc:334-346` logs a
        // non-fatal MessageHandler::error and flips spherical=true; we
        // surface the misconfiguration instead — silently flipping the
        // gravity model under the caller violates Fail Loudly.)
        //
        // `degree == 1` is also rejected: Gottlieb returns zero
        // perturbation for degree < 2 (`calc_nonspherical_with_scratch`),
        // so the per-step `effective_orders` clamp collapses it to 0
        // and the runtime path takes the spherical branch. Accepting
        // `degree == 1` at startup would let that silent fixup propagate
        // to a kernel that produces point-mass acceleration under a
        // control whose configuration *says* non-spherical — exactly the
        // silent-model-change failure mode this gate exists to prevent.
        assert!(
            self.degree >= 2,
            "Non-spherical gravity (spherical=false) requested with degree={} (< 2). \
             Set spherical=true via `GravityControl::new_spherical(...)` \
             for point-mass gravity, or set degree >= 2. degree=1 is meaningless \
             for spherical harmonics (Gottlieb returns zero perturbation for degree < 2).",
            self.degree
        );

        match &source.model {
            GravityModel::SphericalHarmonics(ref data) => {
                // JEOD_INV: GV.04 — degree <= source degree
                // JEOD_INV: GV.19 — source-side degree/order clamp (same check, catalogued separately)
                assert!(
                    self.degree <= data.degree,
                    "Gravity field degree requested ({}) is greater than max gravity field degree ({}).",
                    self.degree, data.degree
                );
                // JEOD_INV: GV.05 — order <= source order
                assert!(
                    self.order <= data.order,
                    "Gravity field order requested ({}) is greater than max gravity field order ({}).",
                    self.order, data.order
                );
            }
            GravityModel::PointMass => {
                panic!(
                    "Non-spherical gravity (spherical=false) is only supported for \
                     SphericalHarmonics gravity models. Set spherical=true for \
                     point-mass gravity sources."
                );
            }
        }

        // JEOD_INV: GV.06 — requested spherical-harmonics order must not exceed requested degree
        assert!(
            self.order <= self.degree,
            "Gravity field order ({}) is greater than gravity field degree ({}).",
            self.order,
            self.degree
        );

        // Gradient validation: JEOD `spherical_harmonics_gravity_controls.cc:395-454`
        // uses MessageHandler::error (non-fatal) and silently auto-corrects
        // invalid ordinals. We treat the same misconfigurations as fatal:
        // an out-of-range gradient ordinal means the caller's expectation
        // of which terms contribute to the gradient tensor diverges from
        // what the kernel will compute, and propagating that divergence
        // produces a wrong (but plausible-looking) torque trajectory.
        if self.gradient {
            // JEOD_INV: GV.08 — gradient_degree <= degree
            assert!(
                self.gradient_degree <= self.degree,
                "Gravity gradient degree ({}) > gravity degree ({}). \
                 Set gradient_degree <= degree, or set gradient_degree=0 to \
                 skip the spherical-harmonics gradient contribution.",
                self.gradient_degree,
                self.degree
            );
            // JEOD_INV: GV.09 — gradient_degree != 1
            // (Gottlieb returns zero perturbation for degree < 2, so a
            // gradient_degree of exactly 1 is meaningless — every n=1 term
            // collapses to point-mass and produces no SH gradient.)
            assert!(
                self.gradient_degree != 1,
                "Gravity gradient degree must not equal 1 (no SH gradient \
                 contribution from degree-1 terms). Set gradient_degree=0 \
                 to skip the SH gradient, or gradient_degree >= 2."
            );
            // JEOD_INV: GV.10 — gradient_order <= gradient_degree
            assert!(
                self.gradient_order <= self.gradient_degree,
                "Gravity gradient order ({}) > gradient degree ({}). \
                 Set gradient_order <= gradient_degree.",
                self.gradient_order,
                self.gradient_degree
            );
            // JEOD_INV: GV.11 — gradient_order <= order
            assert!(
                self.gradient_order <= self.order,
                "Gravity gradient order ({}) > gravity order ({}). \
                 Set gradient_order <= order.",
                self.gradient_order,
                self.order
            );
        }
    }

    /// Returns true if this control's *configuration* selects non-spherical
    /// (spherical-harmonics) computation, i.e. `spherical` is false and
    /// `degree > 0`.
    ///
    /// This is **purely config-based** — it does not consult the source.
    /// For the runtime question "will this control actually compute SH
    /// terms against `source`?", use [`Self::requires_planet_fixed_rotation`].
    /// Examples where `is_nonspherical()` returns true but the runtime
    /// path collapses to spherical:
    /// - source is `GravityModel::PointMass` (no SH data → effective degree = 0)
    /// - the configured degree exceeds the source degree (clamped down to 0)
    /// - the configured degree is 1 (Gottlieb returns zero perturbation for
    ///   degree < 2, so the per-step clamp collapses it to 0)
    pub fn is_nonspherical(&self) -> bool {
        !self.spherical && self.degree > 0
    }

    /// Returns true if evaluating this control against `source` will require
    /// the planet-fixed rotation matrix (i.e., the runtime path is genuinely
    /// non-spherical after clamping). Mirrors the gate inside `evaluate`'s
    /// dispatch.
    ///
    /// Use this for pre-flight checks (e.g., "does this source need a
    /// rotation model wired up?") instead of [`Self::is_nonspherical`],
    /// which is config-only and would over-approximate the requirement.
    pub fn requires_planet_fixed_rotation(&self, source: &GravitySource) -> bool {
        self.effective_orders(source).0 > 0
    }

    /// Compute the effective `(degree, order, gradient_degree, gradient_order)`
    /// quadruple after clamping to the source's bounds. Does not mutate `self`.
    ///
    /// This is the per-step variant of [`Self::check_validity`]: where
    /// `check_validity` is a startup gate that panics on any out-of-range
    /// ordinal (Fail Loudly), `effective_orders` is the runtime path used
    /// by [`Self::evaluate_inner`] on every step. It clamps gracefully
    /// rather than panicking so a control mutated mid-mission, or one
    /// constructed by a test that deliberately bypasses the validation
    /// pipeline, doesn't crash deep inside the spherical-harmonics
    /// kernel. For controls that have passed `check_validity`, the
    /// per-step clamp is a no-op.
    ///
    /// Returns `(0, 0, 0, 0)` for spherical controls, point-mass sources,
    /// or any case where the request collapses to point-mass gravity
    /// (e.g., `spherical=false` against a `GravityModel::PointMass`).
    /// Callers see the spherical branch in [`Self::evaluate_inner`] when
    /// the returned degree is 0.
    // JEOD_INV: GV.04 — degree clamped to source degree (per-step path)
    // JEOD_INV: GV.05 — order clamped to source order
    // JEOD_INV: GV.06 — order clamped to degree
    // JEOD_INV: GV.08 — gradient_degree clamped to degree
    // JEOD_INV: GV.09 — gradient_degree=1 collapses to 0
    // JEOD_INV: GV.10 — gradient_order clamped to gradient_degree
    // JEOD_INV: GV.11 — gradient_order clamped to order
    fn effective_orders(&self, source: &GravitySource) -> (usize, usize, usize, usize) {
        if self.spherical {
            return (0, 0, 0, 0);
        }
        let (src_degree, src_order) = match &source.model {
            GravityModel::SphericalHarmonics(data) => (data.degree, data.order),
            // JEOD_INV: GV.07 — non-spherical against point-mass panics at
            // startup (`check_validity`); here we mirror that with a
            // safety zero quadruple for the rare path where a control
            // bypasses `check_validity`, so `evaluate_inner` takes the
            // spherical branch instead of panicking inside the kernel.
            GravityModel::PointMass => return (0, 0, 0, 0),
        };
        let mut degree = self.degree.min(src_degree);
        // Gottlieb early-returns zero perturbation for degree < 2 (see
        // `calc_nonspherical_with_scratch`), so a configured degree of 1
        // produces no SH contribution. Collapse to 0 here so the runtime
        // predicate `eff_degree > 0` aligns with "the SH path actually does
        // work," and so the planet-fixed rotation matrix is not required
        // for a control whose computation degenerates to point-mass anyway.
        if degree == 1 {
            degree = 0;
        }
        // GV.06: order ≤ degree; GV.05: order ≤ source order
        let order = self.order.min(src_order).min(degree);
        // GV.08: gradient_degree ≤ degree
        let mut gradient_degree = self.gradient_degree.min(degree);
        // GV.09: gradient_degree of exactly 1 is meaningless; collapse to 0
        if gradient_degree == 1 {
            gradient_degree = 0;
        }
        // GV.10: gradient_order ≤ gradient_degree; GV.11: gradient_order ≤ order
        let gradient_order = self.gradient_order.min(gradient_degree).min(order);
        (degree, order, gradient_degree, gradient_order)
    }

    /// Evaluate this gravity control for a single source at the given position.
    ///
    /// Dispatches to spherical (point-mass) or non-spherical (spherical harmonics)
    /// gravity computation based on this control's configuration. For non-spherical
    /// gravity, `t_inertial_pfix` must be `Some` (matching JEOD's requirement that
    /// the planet-fixed frame is subscribed for non-spherical gravity).
    ///
    /// # Arguments
    /// - `source`: the gravity source (mu + model data)
    /// - `position`: body position relative to source center, in inertial frame
    /// - `t_inertial_pfix`: inertial-to-planet-fixed rotation (required for non-spherical)
    ///
    /// # Panics
    /// Panics if non-spherical gravity is requested but `t_inertial_pfix` is `None`.
    // JEOD_INV: GV.13 — gravity source must have inertial frame (planet-fixed rotation required for non-spherical)
    // JEOD_INV: GV.17 — active nonspherical controls subscribe to planet-fixed frame
    pub fn evaluate(
        &self,
        source: &GravitySource,
        position: DVec3,
        t_inertial_pfix: Option<&DMat3>,
        delta_c20: f64,
        has_delta_coeffs: bool,
    ) -> GravityAcceleration {
        self.evaluate_inner(
            source,
            position,
            t_inertial_pfix,
            self.gradient,
            self.gradient_degree,
            self.gradient_order,
            delta_c20,
            has_delta_coeffs,
        )
    }

    /// Like [`evaluate`](Self::evaluate), but passes `compute_gradient=false`
    /// regardless of this control's `gradient` flag.
    ///
    /// This skips the spherical-harmonics gradient tensor computation (the
    /// expensive part). Point-mass acceleration, potential, and point-mass
    /// gradient are still computed internally by `gravitation()` but the
    /// caller typically reads only `.grav_accel`.
    ///
    /// Use this in hot loops (e.g., RK4 inner stages) where only the
    /// gravitational acceleration vector is needed.
    pub fn evaluate_accel_only(
        &self,
        source: &GravitySource,
        position: DVec3,
        t_inertial_pfix: Option<&DMat3>,
        delta_c20: f64,
        has_delta_coeffs: bool,
    ) -> GravityAcceleration {
        self.evaluate_inner(
            source,
            position,
            t_inertial_pfix,
            false,
            0,
            0,
            delta_c20,
            has_delta_coeffs,
        )
    }

    /// Shared dispatch for [`Self::evaluate`] and [`Self::evaluate_accel_only`].
    ///
    /// All four spherical-harmonic ordinals (`degree`, `order`,
    /// `gradient_degree`, `gradient_order`) are clamped to the source's
    /// bounds via [`Self::effective_orders`] before reaching the
    /// `gravitation` kernel. This makes the per-step path safe against
    /// controls that were constructed outside the validation pipeline
    /// or mutated mid-mission — the kernel never sees an out-of-range
    /// request that would panic deep inside the spherical-harmonics
    /// recurrence. For already-validated controls (the common case),
    /// the clamp is a no-op and Tier 3 baselines stay bit-identical.
    // JEOD_INV: GV.13 — gravity source must have inertial frame (planet-fixed rotation required for non-spherical)
    // JEOD_INV: GV.17 — active nonspherical controls subscribe to planet-fixed frame
    #[allow(clippy::too_many_arguments)]
    fn evaluate_inner(
        &self,
        source: &GravitySource,
        position: DVec3,
        t_inertial_pfix: Option<&DMat3>,
        compute_gradient: bool,
        gradient_degree_request: usize,
        gradient_order_request: usize,
        delta_c20: f64,
        has_delta_coeffs: bool,
    ) -> GravityAcceleration {
        let (eff_degree, eff_order, mut eff_grad_degree, mut eff_grad_order) =
            self.effective_orders(source);
        // The caller may pass an explicit `gradient_degree_request` /
        // `gradient_order_request` — `evaluate_accel_only` passes 0/0 to
        // skip the gradient. Honor the caller's request *capped* by the
        // already-clamped ordinals.
        eff_grad_degree = eff_grad_degree.min(gradient_degree_request);
        eff_grad_order = eff_grad_order
            .min(gradient_order_request)
            .min(eff_grad_degree);

        if eff_degree > 0 {
            // Non-spherical path: requires the planet-fixed rotation.
            let rot = t_inertial_pfix.unwrap_or_else(|| {
                panic!(
                    "Non-spherical gravity (degree={}/order={}) requires planet-fixed \
                     rotation matrix. In JEOD, the planet-fixed frame is always \
                     subscribed for non-spherical gravity.",
                    eff_degree, eff_order
                )
            });
            let kernel_out = crate::gravitation(
                source,
                position,
                rot,
                eff_degree,
                eff_order,
                self.perturbing_only,
                compute_gradient,
                eff_grad_degree,
                eff_grad_order,
                delta_c20,
                has_delta_coeffs,
            );
            // The kernel returns SH in planet-fixed; apply the inverse
            // rotation here. `into_inertial` is gated on
            // `compute_gradient`, so the 9-mul/9-add matrix transform
            // on the gradient tensor is skipped when the caller asked
            // for accel-only (the RK4 inner loop via
            // `evaluate_accel_only`). The accel transform itself still
            // runs once per kernel call — `position` differs per
            // substage, so the resulting inertial-frame vector cannot
            // be hoisted across substages even though `t_parent_this`
            // is.
            kernel_out.into_inertial(rot, compute_gradient)
        } else {
            // Point-mass path: the kernel returns the inertial-frame
            // piece directly (`sh_pfix` is `None`), so the
            // `into_inertial` call below short-circuits the rotation
            // entirely — no matrix-vector ops run on this branch.
            let kernel_out = crate::gravitation(
                source,
                position,
                &DMat3::IDENTITY,
                0,
                0,
                self.perturbing_only,
                compute_gradient,
                eff_grad_degree,
                eff_grad_order,
                0.0,   // point-mass: no SH coefficients to modify
                false, // point-mass: no delta coefficients
            );
            // `sh_pfix` is None on this branch; `into_inertial` short-
            // circuits the rotation and returns the point-mass piece
            // unchanged (or `GravityAcceleration::default()` when
            // `perturbing_only` skips the point-mass term).
            kernel_out.into_inertial(&DMat3::IDENTITY, compute_gradient)
        }
    }
}

/// Typed sibling of [`GravityControl`].
///
/// Field-for-field parity with the untyped form, except the four
/// spherical-harmonic ordinals (`degree`, `order`, `gradient_degree`,
/// `gradient_order`) carry the [`HarmonicDegree`] newtype so the
/// compiler distinguishes them from angular `Angle` or dimensionless
/// `Ratio`.
///
/// Cross-field invariants like `degree <= source.degree` (JEOD
/// `GV.03`–`GV.11`) remain runtime-checked via
/// [`GravityControlTyped::check_validity`] (which delegates to the
/// untyped [`GravityControl::check_validity`]) — the type system
/// can prove ordinals are distinct kinds, not that one specific
/// ordinal is bounded by another's runtime value.
#[derive(Debug, Clone)]
pub struct GravityControlTyped {
    /// The source's inertial-frame identity (see [`GravityControl::source`]).
    pub source: FrameUid,
    /// Compute the gravity-gradient tensor in addition to acceleration.
    pub gradient: bool,
    /// Use only point-mass gravity for this source.
    pub spherical: bool,
    /// Spherical-harmonics degree to use (must be ≤ source degree).
    pub degree: HarmonicDegree,
    /// Spherical-harmonics order to use (must be ≤ source order, ≤ degree).
    pub order: HarmonicDegree,
    /// Exclude `n=0,1` (point-mass) terms from the SH evaluation.
    pub perturbing_only: bool,
    /// Degree for gradient computation (must be ≤ degree).
    pub gradient_degree: HarmonicDegree,
    /// Order for gradient computation (must be ≤ order, ≤ gradient_degree).
    pub gradient_order: HarmonicDegree,
    /// Treat this source as a third-body and use differential acceleration.
    pub differential: bool,
    /// Use Battin's method for differential gravity (third-body only).
    pub battin_method: bool,
    /// Apply post-Newtonian relativistic correction.
    pub relativistic: bool,
}

impl GravityControlTyped {
    /// Spherical (point-mass) typed control. See
    /// [`GravityControl::new_spherical`] for the `impl Into<FrameUid>`
    /// rationale.
    pub fn new_spherical(source: impl Into<FrameUid>, gradient: GravityGradient) -> Self {
        Self {
            source: source.into(),
            gradient: gradient.as_bool(),
            spherical: true,
            degree: HarmonicDegree::default(),
            order: HarmonicDegree::default(),
            perturbing_only: false,
            gradient_degree: HarmonicDegree::default(),
            gradient_order: HarmonicDegree::default(),
            differential: false,
            battin_method: false,
            relativistic: false,
        }
    }

    /// Non-spherical (spherical-harmonics) typed control. See
    /// [`GravityControl::new_spherical`] for the `impl Into<FrameUid>`
    /// rationale.
    pub fn new_nonspherical(
        source: impl Into<FrameUid>,
        degree: HarmonicDegree,
        order: HarmonicDegree,
        gradient: GravityGradient,
    ) -> Self {
        Self {
            source: source.into(),
            gradient: gradient.as_bool(),
            spherical: false,
            degree,
            order,
            perturbing_only: false,
            gradient_degree: HarmonicDegree::default(),
            gradient_order: HarmonicDegree::default(),
            differential: false,
            battin_method: false,
            relativistic: false,
        }
    }

    /// Third-body (point-mass + differential) typed control. See
    /// [`GravityControl::new_spherical`] for the `impl Into<FrameUid>`
    /// rationale.
    pub fn new_third_body(source: impl Into<FrameUid>) -> Self {
        Self {
            source: source.into(),
            gradient: false,
            spherical: true,
            degree: HarmonicDegree::default(),
            order: HarmonicDegree::default(),
            perturbing_only: false,
            gradient_degree: HarmonicDegree::default(),
            gradient_order: HarmonicDegree::default(),
            differential: true,
            battin_method: false,
            relativistic: false,
        }
    }
}

impl GravityControlTyped {
    /// Validate this typed control against its gravity source.
    ///
    /// Delegates to [`GravityControl::check_validity`] on the untyped
    /// projection — runtime-checked invariants (`GV.03`–`GV.11`)
    /// stay in the canonical f64 path. The validator panics on any
    /// misconfiguration (see [`GravityControl::check_validity`] for
    /// the exhaustive list); on success the typed control is
    /// unchanged.
    // JEOD_INV: GV.03 — check_validity() called on degree/order mutation
    pub fn check_validity(&self, source: &GravitySource) {
        let untyped = self.to_untyped();
        untyped.check_validity(source);
    }

    /// Drop the [`HarmonicDegree`] newtypes and emit the untyped
    /// storage form. Cross-field invariants (GV.03–GV.11) remain
    /// runtime-checked via the resulting
    /// [`GravityControl::check_validity`].
    pub fn to_untyped(&self) -> GravityControl {
        GravityControl {
            source: self.source.clone(),
            gradient: self.gradient,
            spherical: self.spherical,
            degree: self.degree.get(),
            order: self.order.get(),
            perturbing_only: self.perturbing_only,
            gradient_degree: self.gradient_degree.get(),
            gradient_order: self.gradient_order.get(),
            differential: self.differential,
            battin_method: self.battin_method,
            relativistic: self.relativistic,
        }
    }

    /// Wrap an untyped [`GravityControl`] as typed. Lossless conversion.
    pub fn from_untyped_unchecked(c: &GravityControl) -> Self {
        Self {
            source: c.source.clone(),
            gradient: c.gradient,
            spherical: c.spherical,
            degree: HarmonicDegree::from(c.degree),
            order: HarmonicDegree::from(c.order),
            perturbing_only: c.perturbing_only,
            gradient_degree: HarmonicDegree::from(c.gradient_degree),
            gradient_order: HarmonicDegree::from(c.gradient_order),
            differential: c.differential,
            battin_method: c.battin_method,
            relativistic: c.relativistic,
        }
    }
}

// No `Default` impls: a control's source identity must be supplied
// explicitly — there is no meaningful "default source", and minting an
// identity by accident is exactly the failure mode `FrameUid` keying
// exists to kill (issue #668; `FrameUid` itself has no `Default` for
// the same reason).

#[cfg(test)]
mod tests {
    use super::*;
    use crate::spherical_harmonics_gravity_source::SphericalHarmonicsData;
    use astrodyn_quantities::frame_descriptor::{FrameClass, FrameRole, Namespace, Tag};

    /// A fixed source identity for tests that only need *a* uid — the
    /// clamping/validation logic under test never resolves it.
    fn test_uid() -> FrameUid {
        FrameUid::external(
            Namespace(2),
            FrameClass::PlanetInertial,
            FrameRole::Primary,
            Tag::Named("test-source".into()),
        )
    }

    /// Build an in-memory spherical-harmonics source with all zero
    /// coefficients (degenerate, but sufficient for exercising the
    /// degree/order clamping logic — the actual numeric values don't
    /// matter; we're checking that the kernel runs without panicking
    /// and returns a finite acceleration). The coefficient arrays
    /// follow the triangular `cnm[n].len() == n + 1` shape that
    /// `SphericalHarmonicsData::new` expects.
    fn dummy_sh_source(degree: usize, order: usize) -> GravitySource {
        let mu = 3.986_004_415e14;
        let radius = 6_378_137.0;
        let cnm: Vec<Vec<f64>> = (0..=degree).map(|n| vec![0.0_f64; n + 1]).collect();
        let snm: Vec<Vec<f64>> = (0..=degree).map(|n| vec![0.0_f64; n + 1]).collect();
        let data = SphericalHarmonicsData::new(degree, order, radius, mu, cnm, snm, true, 0.0);
        GravitySource {
            mu,
            model: GravityModel::SphericalHarmonics(Box::new(data)),
        }
    }

    /// `effective_orders` for a spherical control returns all zeros
    /// regardless of the source's degree.
    #[test]
    fn effective_orders_spherical_returns_zeros() {
        let src = dummy_sh_source(8, 8);
        let ctrl = GravityControl::new_spherical(test_uid(), GravityGradient::Skip);
        assert_eq!(ctrl.effective_orders(&src), (0, 0, 0, 0));
    }

    /// `effective_orders` against a `PointMass` source collapses any
    /// non-spherical request to zeros (mirrors GV.07 startup auto-correct).
    #[test]
    fn effective_orders_against_point_mass_collapses_to_zero() {
        let src = GravitySource {
            mu: 3.986_004_415e14,
            model: GravityModel::PointMass,
        };
        let ctrl = GravityControl {
            spherical: false,
            degree: 8,
            order: 8,
            ..GravityControl::new_spherical(test_uid(), GravityGradient::Skip)
        };
        assert_eq!(ctrl.effective_orders(&src), (0, 0, 0, 0));
    }

    /// Out-of-range degree/order are clamped down to the source's bounds
    /// and to each other (GV.04, GV.05, GV.06).
    #[test]
    fn effective_orders_clamps_degree_order_to_source() {
        let src = dummy_sh_source(8, 8);
        let ctrl = GravityControl {
            spherical: false,
            degree: 100,
            order: 100,
            ..GravityControl::new_spherical(test_uid(), GravityGradient::Skip)
        };
        assert_eq!(ctrl.effective_orders(&src), (8, 8, 0, 0));
    }

    /// Order > degree clamps order down to degree (GV.06).
    #[test]
    fn effective_orders_clamps_order_to_degree() {
        let src = dummy_sh_source(8, 8);
        let ctrl = GravityControl {
            spherical: false,
            degree: 4,
            order: 8,
            ..GravityControl::new_spherical(test_uid(), GravityGradient::Skip)
        };
        assert_eq!(ctrl.effective_orders(&src), (4, 4, 0, 0));
    }

    /// Gradient ordinals are clamped: `gradient_degree=1` collapses to
    /// 0 (GV.09), `gradient_order ≤ gradient_degree` (GV.10), and
    /// `gradient_order ≤ order` (GV.11).
    #[test]
    fn effective_orders_clamps_gradient_ordinals() {
        let src = dummy_sh_source(8, 6);
        let ctrl = GravityControl {
            spherical: false,
            degree: 8,
            order: 4,
            gradient: true,
            gradient_degree: 1, // → collapses to 0
            gradient_order: 5,  // > gradient_degree, > order → clamped
            ..GravityControl::new_spherical(test_uid(), GravityGradient::Skip)
        };
        // After clamping: degree=8, order=4 (≤ src.order=6), gradient_degree=0, gradient_order=0.
        assert_eq!(ctrl.effective_orders(&src), (8, 4, 0, 0));
    }

    /// Regression test for H4: a control with out-of-range
    /// `gradient_degree` constructed without going through `check_validity`
    /// previously panicked deep inside `gravitation`. After the fix,
    /// `evaluate_inner` clamps gracefully and returns a finite acceleration.
    #[test]
    fn evaluate_does_not_panic_on_out_of_range_gradient_degree() {
        let src = dummy_sh_source(4, 4);
        let ctrl = GravityControl {
            spherical: false,
            degree: 4,
            order: 4,
            gradient: true,
            gradient_degree: 100, // wildly out of range
            gradient_order: 100,
            ..GravityControl::new_spherical(test_uid(), GravityGradient::Skip)
        };
        let pos = DVec3::new(7_000_000.0, 0.0, 0.0);
        let rot = DMat3::IDENTITY;
        // No panic; the kernel sees clamped (gradient_degree, gradient_order)=(4, 4).
        let result = ctrl.evaluate(&src, pos, Some(&rot), 0.0, false);
        assert!(result.grav_accel.is_finite());
    }

    /// Regression test for H4: a control with out-of-range `degree` /
    /// `order` against a real spherical-harmonics source no longer
    /// panics deep in `calc_nonspherical`.
    #[test]
    fn evaluate_does_not_panic_on_out_of_range_degree() {
        let src = dummy_sh_source(4, 4);
        let ctrl = GravityControl {
            spherical: false,
            degree: 100,
            order: 100,
            ..GravityControl::new_spherical(test_uid(), GravityGradient::Skip)
        };
        let pos = DVec3::new(7_000_000.0, 0.0, 0.0);
        let rot = DMat3::IDENTITY;
        let result = ctrl.evaluate(&src, pos, Some(&rot), 0.0, false);
        assert!(result.grav_accel.is_finite());
    }

    /// `degree == 1` is meaningless for spherical harmonics: Gottlieb
    /// returns zero perturbation for `degree < 2`. `effective_orders`
    /// collapses such a control to all-zeros so the runtime path takes
    /// the spherical branch and does not require the planet-fixed
    /// rotation matrix.
    #[test]
    fn effective_orders_collapses_degree_one_to_zero() {
        let src = dummy_sh_source(8, 8);
        let ctrl = GravityControl {
            spherical: false,
            degree: 1,
            order: 1,
            ..GravityControl::new_spherical(test_uid(), GravityGradient::Skip)
        };
        assert_eq!(ctrl.effective_orders(&src), (0, 0, 0, 0));
    }

    /// `requires_planet_fixed_rotation` returns false when the runtime
    /// path collapses to spherical (`PointMass` source, `degree=0`,
    /// `degree=1`), even if `is_nonspherical()` is true. Pre-flight
    /// checks must use this method, not `is_nonspherical()`.
    #[test]
    fn requires_rotation_false_when_runtime_collapses() {
        let sh = dummy_sh_source(8, 8);
        let pm = GravitySource {
            mu: 3.986_004_415e14,
            model: GravityModel::PointMass,
        };
        let degree_one = GravityControl {
            spherical: false,
            degree: 1,
            order: 1,
            ..GravityControl::new_spherical(test_uid(), GravityGradient::Skip)
        };
        let against_pm = GravityControl {
            spherical: false,
            degree: 8,
            order: 8,
            ..GravityControl::new_spherical(test_uid(), GravityGradient::Skip)
        };
        assert!(degree_one.is_nonspherical()); // config says yes
        assert!(!degree_one.requires_planet_fixed_rotation(&sh)); // runtime says no
        assert!(against_pm.is_nonspherical()); // config says yes
        assert!(!against_pm.requires_planet_fixed_rotation(&pm)); // runtime says no

        let real_sh = GravityControl {
            spherical: false,
            degree: 4,
            order: 4,
            ..GravityControl::new_spherical(test_uid(), GravityGradient::Skip)
        };
        assert!(real_sh.requires_planet_fixed_rotation(&sh));
    }

    /// Regression test for the degree-1 panic (review feedback on PR #182):
    /// a non-spherical control with `degree=1` against a SH source no
    /// longer panics when `t_inertial_pfix` is `None`, because
    /// `effective_orders` collapses degree=1 to 0 (Gottlieb produces zero
    /// perturbation for degree < 2 anyway).
    #[test]
    fn evaluate_degree_one_does_not_require_rotation() {
        let src = dummy_sh_source(8, 8);
        let ctrl = GravityControl {
            spherical: false,
            degree: 1,
            order: 1,
            ..GravityControl::new_spherical(test_uid(), GravityGradient::Skip)
        };
        let pos = DVec3::new(7_000_000.0, 0.0, 0.0);
        // No rotation matrix supplied; would have panicked previously.
        let result = ctrl.evaluate(&src, pos, None, 0.0, false);
        // Should match point-mass gravity from the source (the SH branch
        // contributes zero for degree < 2 and `perturbing_only` is false).
        assert!(result.grav_accel.is_finite());
        assert!(result.grav_accel.length() > 0.0);
    }

    // ---- check_validity fail-loud sites ------------------------------
    //
    // `check_validity` panics on every misconfiguration that JEOD's
    // `MessageHandler::error` would have silently auto-corrected. Each
    // test below pins one panic class so a future regression that
    // re-introduces silent auto-correction trips immediately.

    /// `spherical=false` with `degree=0` is a misconfiguration: the
    /// caller meant point-mass gravity but did not flip `spherical`.
    /// JEOD silently flips it; we panic so the operator sees the
    /// inconsistency.
    // JEOD_INV: GV.07 — negative test for degree < 2 with spherical=false
    #[test]
    #[should_panic(
        expected = "Non-spherical gravity (spherical=false) requested with degree=0 (< 2)"
    )]
    fn check_validity_panics_on_zero_degree_with_spherical_false() {
        let src = dummy_sh_source(8, 8);
        let ctrl = GravityControl {
            spherical: false,
            degree: 0,
            order: 0,
            ..GravityControl::new_spherical(test_uid(), GravityGradient::Skip)
        };
        ctrl.check_validity(&src);
    }

    /// `spherical=false` with `degree=1` is a misconfiguration: the per-step
    /// `effective_orders` clamp collapses degree=1 to 0 (Gottlieb returns
    /// zero perturbation for degree < 2), which would silently change the
    /// gravity model under the operator. The startup gate rejects it so the
    /// inconsistency between the configured `spherical=false` and the
    /// effectively-point-mass runtime path surfaces immediately.
    // JEOD_INV: GV.07 — negative test for degree == 1 with spherical=false
    #[test]
    #[should_panic(
        expected = "Non-spherical gravity (spherical=false) requested with degree=1 (< 2)"
    )]
    fn check_validity_panics_on_degree_one_with_spherical_false() {
        let src = dummy_sh_source(8, 8);
        let ctrl = GravityControl {
            spherical: false,
            degree: 1,
            order: 1,
            ..GravityControl::new_spherical(test_uid(), GravityGradient::Skip)
        };
        ctrl.check_validity(&src);
    }

    /// `gradient_degree > degree`: JEOD clamps; we panic — clamping
    /// silently changes which SH terms contribute to the gradient
    /// tensor and produces a torque trajectory that differs from what
    /// the operator requested.
    // JEOD_INV: GV.08 — negative test for gradient_degree > degree
    #[test]
    #[should_panic(expected = "Gravity gradient degree (12) > gravity degree (4)")]
    fn check_validity_panics_on_gradient_degree_above_degree() {
        let src = dummy_sh_source(8, 8);
        let ctrl = GravityControl {
            spherical: false,
            degree: 4,
            order: 4,
            gradient: true,
            gradient_degree: 12,
            gradient_order: 0,
            ..GravityControl::new_spherical(test_uid(), GravityGradient::Skip)
        };
        ctrl.check_validity(&src);
    }

    /// `gradient_degree == 1`: Gottlieb returns zero perturbation for
    /// degree < 2, so a gradient_degree of exactly 1 is meaningless.
    /// JEOD resets to 0; we panic so the misconfiguration surfaces.
    // JEOD_INV: GV.09 — negative test for gradient_degree == 1
    #[test]
    #[should_panic(expected = "Gravity gradient degree must not equal 1")]
    fn check_validity_panics_on_gradient_degree_equal_one() {
        let src = dummy_sh_source(8, 8);
        let ctrl = GravityControl {
            spherical: false,
            degree: 4,
            order: 4,
            gradient: true,
            gradient_degree: 1,
            gradient_order: 0,
            ..GravityControl::new_spherical(test_uid(), GravityGradient::Skip)
        };
        ctrl.check_validity(&src);
    }

    /// `gradient_order > gradient_degree`: JEOD clamps; we panic.
    // JEOD_INV: GV.10 — negative test for gradient_order > gradient_degree
    #[test]
    #[should_panic(expected = "Gravity gradient order (5) > gradient degree (2)")]
    fn check_validity_panics_on_gradient_order_above_gradient_degree() {
        let src = dummy_sh_source(8, 8);
        let ctrl = GravityControl {
            spherical: false,
            degree: 4,
            order: 4,
            gradient: true,
            gradient_degree: 2,
            gradient_order: 5,
            ..GravityControl::new_spherical(test_uid(), GravityGradient::Skip)
        };
        ctrl.check_validity(&src);
    }

    /// `gradient_order > order`: JEOD clamps; we panic. The
    /// `gradient_order` check fires before the `> order` check when
    /// both fail, so this test sets `gradient_order` *below*
    /// `gradient_degree` and *above* `order` to isolate GV.11.
    // JEOD_INV: GV.11 — negative test for gradient_order > order
    #[test]
    #[should_panic(expected = "Gravity gradient order (4) > gravity order (2)")]
    fn check_validity_panics_on_gradient_order_above_order() {
        let src = dummy_sh_source(8, 8);
        let ctrl = GravityControl {
            spherical: false,
            degree: 8,
            order: 2,
            gradient: true,
            gradient_degree: 4,
            gradient_order: 4,
            ..GravityControl::new_spherical(test_uid(), GravityGradient::Skip)
        };
        ctrl.check_validity(&src);
    }

    /// A control that satisfies every invariant must not panic. Pinned
    /// alongside the panic tests above so a future overzealous tightening
    /// of `check_validity` is caught immediately rather than at the next
    /// downstream test run.
    #[test]
    fn check_validity_accepts_well_formed_control() {
        let src = dummy_sh_source(8, 8);
        let ctrl = GravityControl {
            spherical: false,
            degree: 4,
            order: 4,
            gradient: true,
            gradient_degree: 4,
            gradient_order: 4,
            ..GravityControl::new_spherical(test_uid(), GravityGradient::Skip)
        };
        ctrl.check_validity(&src); // does not panic
    }

    /// The typed sibling delegates to the untyped validator; the same
    /// misconfiguration must panic on the typed surface too.
    // JEOD_INV: GV.07 — negative test for the typed surface
    #[test]
    #[should_panic(expected = "Non-spherical gravity (spherical=false) requested with degree=0")]
    fn typed_check_validity_panics_on_zero_degree_with_spherical_false() {
        use astrodyn_quantities::aliases::HarmonicDegree;
        let src = dummy_sh_source(8, 8);
        let ctrl = GravityControlTyped {
            spherical: false,
            degree: HarmonicDegree::from(0_usize),
            order: HarmonicDegree::from(0_usize),
            ..GravityControlTyped::new_spherical(test_uid(), GravityGradient::Skip)
        };
        ctrl.check_validity(&src);
    }

    /// `degree > source.degree`: the requested expansion exceeds the
    /// stored coefficient table. JEOD clamps via
    /// `MessageHandler::error`; we panic so a control configured for
    /// (say) 70x70 against an 8x8 source does not silently degenerate
    /// to 8x8 and produce a quietly-wrong trajectory.
    // JEOD_INV: GV.04 — negative test: degree > source.degree
    #[test]
    #[should_panic(
        expected = "Gravity field degree requested (12) is greater than max gravity field degree (8)."
    )]
    fn gv_04_panics_on_degree_above_source_degree() {
        let src = dummy_sh_source(8, 8);
        let ctrl = GravityControl {
            spherical: false,
            degree: 12,
            order: 8,
            ..GravityControl::new_spherical(test_uid(), GravityGradient::Skip)
        };
        ctrl.check_validity(&src);
    }

    /// `order > source.order`: the requested order exceeds the stored
    /// coefficient table. Source is built with `order < degree` so the
    /// GV.05 check fires before the GV.06 `order <= degree` check.
    // JEOD_INV: GV.05 — negative test: order > source.order
    #[test]
    #[should_panic(
        expected = "Gravity field order requested (6) is greater than max gravity field order (4)."
    )]
    fn gv_05_panics_on_order_above_source_order() {
        // Source: degree=8, order=4 (sectorial truncation).
        let src = dummy_sh_source(8, 4);
        let ctrl = GravityControl {
            spherical: false,
            degree: 6,
            order: 6,
            ..GravityControl::new_spherical(test_uid(), GravityGradient::Skip)
        };
        ctrl.check_validity(&src);
    }

    /// `order > degree`: order exceeds the requested degree even though
    /// both are within the source's bounds. JEOD clamps; we panic — a
    /// silently-clamped order would change which sectorial harmonics
    /// contribute to acceleration.
    // JEOD_INV: GV.06 — negative test: order > degree
    #[test]
    #[should_panic(expected = "Gravity field order (5) is greater than gravity field degree (4).")]
    fn gv_06_panics_on_order_above_degree() {
        let src = dummy_sh_source(8, 8);
        let ctrl = GravityControl {
            spherical: false,
            degree: 4,
            order: 5,
            ..GravityControl::new_spherical(test_uid(), GravityGradient::Skip)
        };
        ctrl.check_validity(&src);
    }

    /// Source-side degree clamp at initialization: a request whose
    /// degree exceeds the stored coefficient table fails before the
    /// kernel runs. JEOD's
    /// `SphericalHarmonicsGravitySource::initialize_source` clamps
    /// `degree`/`order` to `max_degree` at startup; we enforce the same
    /// bound at `check_validity` and fail loudly rather than clamping,
    /// so an out-of-table request surfaces as a configuration error.
    /// Shares the panic site with GV.04 — same `assert!` — and is
    /// catalogued separately to mirror JEOD's source-side line.
    // JEOD_INV: GV.19 — negative test: source-side degree clamp
    #[test]
    #[should_panic(
        expected = "Gravity field degree requested (50) is greater than max gravity field degree (8)."
    )]
    fn gv_19_panics_on_degree_above_source_max() {
        let src = dummy_sh_source(8, 8);
        let ctrl = GravityControl {
            spherical: false,
            degree: 50,
            order: 8,
            ..GravityControl::new_spherical(test_uid(), GravityGradient::Skip)
        };
        ctrl.check_validity(&src);
    }

    /// Active nonspherical control must subscribe to the planet-fixed
    /// rotation: `evaluate_inner` panics when `eff_degree > 0` and the
    /// caller passes `None` for the rotation matrix. JEOD subscribes
    /// to the planet-fixed frame unconditionally for non-spherical
    /// gravity (`gravity_controls.cc::initialize_control` registers the
    /// subscription as a structural pre-step). We surface the missing
    /// rotation as a fatal error rather than producing zero / NaN
    /// acceleration silently. Shares the runtime enforcement site
    /// with GV.13 but the catalog tracks the subscription invariant
    /// separately.
    // JEOD_INV: GV.17 — negative test: active nonspherical control without rotation
    #[test]
    #[should_panic(
        expected = "Non-spherical gravity (degree=4/order=4) requires planet-fixed rotation matrix."
    )]
    fn gv_17_panics_on_nonspherical_without_rotation() {
        let src = dummy_sh_source(8, 8);
        let ctrl = GravityControl {
            spherical: false,
            degree: 4,
            order: 4,
            ..GravityControl::new_spherical(test_uid(), GravityGradient::Skip)
        };
        // `evaluate` (not `accumulate_gravity`) drives the kernel-level
        // assert in `evaluate_inner` directly; this test pins the
        // panic site closest to where the rotation matrix is consumed.
        let pos = DVec3::new(7_000_000.0, 0.0, 0.0);
        let _ = ctrl.evaluate(&src, pos, None, 0.0, false);
    }

    // ---- proptest round-trips (#398) ----------------------------------

    use proptest::prelude::*;

    fn arb_gravity_control() -> impl Strategy<Value = GravityControl> {
        (
            // Identity entropy: arbitrary tag names in an external
            // namespace, so the round-trips cover the uid field too
            // (issue #668 — the source reference is a FrameUid value).
            "[a-z]{1,8}",
            any::<bool>(),
            any::<bool>(),
            0usize..=64,
            0usize..=64,
            any::<bool>(),
            0usize..=64,
            0usize..=64,
            any::<bool>(),
            any::<bool>(),
            any::<bool>(),
        )
            .prop_map(
                |(
                    source_name,
                    gradient,
                    spherical,
                    degree,
                    order,
                    perturbing_only,
                    gradient_degree,
                    gradient_order,
                    differential,
                    battin_method,
                    relativistic,
                )| GravityControl {
                    source: FrameUid::external(
                        Namespace(2),
                        FrameClass::PlanetInertial,
                        FrameRole::Primary,
                        Tag::Named(source_name.into()),
                    ),
                    gradient,
                    spherical,
                    degree,
                    order,
                    perturbing_only,
                    gradient_degree,
                    gradient_order,
                    differential,
                    battin_method,
                    relativistic,
                },
            )
    }

    proptest! {
        #[test]
        fn round_trip_gravity_control_untyped_typed_untyped(orig in arb_gravity_control()) {
            let typed = GravityControlTyped::from_untyped_unchecked(&orig);
            prop_assert_eq!(typed.to_untyped(), orig);
        }

        #[test]
        fn round_trip_gravity_control_typed_untyped_typed(orig in arb_gravity_control()) {
            let typed = GravityControlTyped::from_untyped_unchecked(&orig);
            let lifted = GravityControlTyped::from_untyped_unchecked(&typed.to_untyped());
            prop_assert_eq!(lifted.to_untyped(), typed.to_untyped());
        }
    }
}