Added equinoctial elements. Also a little cleanup

benches/kepler.rs

@@ -8,6 +8,7 @@ use std::hint::black_box;
 
 use astrokit::body::Body;
 use astrokit::cartesian::Cartesian;
+use astrokit::equinoctial::Equinoctial;
 use astrokit::frames::Frame;
 use astrokit::kepler::{self, Degeneracy, Keplerian};
 use criterion::{criterion_group, criterion_main, Criterion};
@@ -121,11 +122,44 @@ fn bench_invariants(c: &mut Criterion) {
     group.finish();
 }
 
+/// Equinoctial against Keplerian on the same states.
+///
+/// The interesting comparison is GEO: Keplerian takes its doubly-degenerate
+/// branch there and equinoctial takes no branch at all.
+fn bench_equinoctial(c: &mut Criterion) {
+    let mut group = c.benchmark_group("equinoctial");
+
+    let leo = leo();
+    let geo = geo();
+
+    group.bench_function("from_cartesian/leo", |b| {
+        b.iter(|| Equinoctial::from_cartesian(&black_box(leo)))
+    });
+    group.bench_function("from_cartesian/geo", |b| {
+        b.iter(|| Equinoctial::from_cartesian(&black_box(geo)))
+    });
+
+    let elements = Equinoctial::from_cartesian(&leo);
+    group.bench_function("to_cartesian", |b| {
+        b.iter(|| black_box(elements).to_cartesian())
+    });
+    group.bench_function("round_trip", |b| {
+        b.iter(|| Equinoctial::from_cartesian(&black_box(leo)).to_cartesian())
+    });
+
+    let keplerian = leo.keplerian();
+    group.bench_function("from_keplerian", |b| {
+        b.iter(|| Equinoctial::from_keplerian(&black_box(keplerian)))
+    });
+    group.finish();
+}
+
 criterion_group!(
     benches,
     bench_scalar_relations,
     bench_anomalies,
     bench_conversions,
-    bench_invariants
+    bench_invariants,
+    bench_equinoctial
 );
 criterion_main!(benches);

benches/propagation.rs

@@ -117,9 +117,9 @@ fn bench_one_period(c: &mut Criterion) {
 
 /// The cost of a state transition matrix relative to a plain propagation.
 ///
-/// Stage 1 measured 10.6x for a 6x6 on a two-body arc. Central differences would
-/// need thirteen propagations and lose about half the digits, so anything near
-/// that ratio is a win; this tracks it as the force models get heavier.
+/// Central differences would need thirteen propagations for a 6x6 and lose about
+/// half the available digits. This tracks the ratio as the force models get
+/// heavier.
 #[cfg(feature = "autodiff")]
 fn bench_state_transition_matrix(c: &mut Criterion) {
     use nalgebra::{Matrix6, Vector6};

readme.md

@@ -10,9 +10,9 @@ also differentiate.
 
 The crate is still being built to a milestone plan:
 
-- **M1** — the study-critical subset: two-body and spherical-harmonic gravity,
-  drag, third-body, cannonball SRP, GCRF/ITRF frames, ground stations and the
-  four measurement types, with a validation suite against Orekit reference
+- **M1** — the core subset: two-body and spherical-harmonic gravity, drag,
+  third-body, cannonball SRP, GCRF/ITRF frames, ground stations and the four
+  measurement types, with a validation suite against Orekit reference
   trajectories.
 - **M2** — general-purpose completion: higher-degree gravity fields, more
   density models, box-wing SRP and eclipses, maneuvers, covariance propagation,

src/equinoctial.rs

@@ -0,0 +1,603 @@
+//! Modified equinoctial elements.
+//!
+//! Unlike [`Keplerian`], these have no degenerate case on a circular or equatorial
+//! orbit. The one singularity is exactly retrograde motion, `i = 180`, which this
+//! module refuses rather than returning a diverging `tan(i/2)`.
+//!
+//! Conversions to and from Cartesian go directly, through the equinoctial basis
+//! vectors, not via [`Keplerian`] — routing through the degenerate representation
+//! would defeat the purpose.
+//!
+//! Longitudes are *true*, matching [`Keplerian`]'s use of true anomaly.
+//!
+//! Walker, Ireland and Owens, "A set of modified equinoctial orbit elements",
+//! Celestial Mechanics 36 (1985).
+
+use nalgebra::Vector3;
+use ordinary_diffeq::scalar::Scalar;
+
+use crate::cartesian::Cartesian;
+use crate::frames::Frame;
+use crate::kepler::{self, wrap_two_pi, Degeneracy, Keplerian};
+use crate::state::{HasPosition, HasVelocity, State};
+
+/// How close to exactly retrograde an orbit may be before the prograde
+/// convention is refused.
+///
+/// `h` and `k` carry `tan(i/2)`, which diverges at `i = 180`. At this tolerance
+/// they reach about 1e7, which is representable but already meaningless.
+pub const RETROGRADE_TOLERANCE: f64 = 1e-7;
+
+/// Modified equinoctial orbital elements.
+#[derive(Clone, Copy, Debug, PartialEq)]
+pub struct Equinoctial<T: Scalar> {
+    /// Semi-latus rectum `a(1 - e^2)`, m.
+    pub semi_latus_rectum: T,
+    /// `e cos(aop + raan)`.
+    pub f: T,
+    /// `e sin(aop + raan)`.
+    pub g: T,
+    /// `tan(i/2) cos(raan)`.
+    pub h: T,
+    /// `tan(i/2) sin(raan)`.
+    pub k: T,
+    /// True longitude `raan + aop + true anomaly`, rad.
+    pub true_longitude: T,
+    /// The frame these describe motion in. Supplies `mu`.
+    pub frame: Frame,
+}
+
+impl<T: Scalar> Equinoctial<T> {
+    /// Elements from a Cartesian state, without going through [`Keplerian`].
+    ///
+    /// # Panics
+    ///
+    /// If the frame is not inertial, or if the orbit is within
+    /// [`RETROGRADE_TOLERANCE`] of exactly retrograde.
+    pub fn from_cartesian(state: &Cartesian<T>) -> Self {
+        state.require_inertial("converting to equinoctial elements");
+        let (position, velocity) = (state.position, state.velocity);
+        let mu = state.frame.mu();
+
+        let momentum = position.cross(&velocity);
+        let normal = momentum.normalize();
+
+        // h and k come straight from the orbit normal: for a prograde orbit the
+        // normal is [sin i sin raan, -sin i cos raan, cos i], and
+        // tan(i/2) = sin i / (1 + cos i).
+        let denominator = T::lift(1.0) + normal[2];
+        assert!(
+            denominator.re() > RETROGRADE_TOLERANCE,
+            "equinoctial elements in the prograde convention are singular at \
+             i = 180 degrees; this orbit is within {RETROGRADE_TOLERANCE:e} of it"
+        );
+        let h = -normal[1] / denominator;
+        let k = normal[0] / denominator;
+
+        // The semi-latus rectum needs no eccentricity: p = |h|^2 / mu.
+        let semi_latus_rectum = momentum.norm_squared() / mu;
+
+        // f and g are the eccentricity vector's components in the equinoctial
+        // basis, and the true longitude is the position's angle in it.
+        let (f_hat, g_hat) = Self::basis(h, k);
+        let eccentricity_vector = (position
+            * (velocity.norm_squared() - T::lift(mu) / position.norm())
+            - velocity * position.dot(&velocity))
+            / T::lift(mu);
+
+        Self {
+            semi_latus_rectum,
+            f: eccentricity_vector.dot(&f_hat),
+            g: eccentricity_vector.dot(&g_hat),
+            h,
+            k,
+            true_longitude: wrap_two_pi(position.dot(&g_hat).atan2(position.dot(&f_hat))),
+            frame: state.frame,
+        }
+    }
+
+    /// The equinoctial basis vectors, in the frame's axes.
+    ///
+    /// `f_hat` is the direction true longitude is measured from and `g_hat`
+    /// completes the orbital plane. Both are unit vectors and their cross product
+    /// is the orbit normal.
+    fn basis(h: T, k: T) -> (Vector3<T>, Vector3<T>) {
+        let one = T::lift(1.0);
+        let scale = one / (one + h.powi(2) + k.powi(2));
+        (
+            Vector3::new(
+                (one - k.powi(2) + h.powi(2)) * scale,
+                h * k * 2.0 * scale,
+                -k * 2.0 * scale,
+            ),
+            Vector3::new(
+                h * k * 2.0 * scale,
+                (one + k.powi(2) - h.powi(2)) * scale,
+                h * 2.0 * scale,
+            ),
+        )
+    }
+
+    /// Checks that the basis is orthonormal and right-handed about the orbit
+    /// normal implied by `h` and `k`. Used by the tests.
+    #[cfg(test)]
+    fn basis_normal(h: T, k: T) -> Vector3<T> {
+        let (f_hat, g_hat) = Self::basis(h, k);
+        f_hat.cross(&g_hat)
+    }
+
+    /// A Cartesian state, in the frame these elements carry.
+    pub fn to_cartesian(&self) -> Cartesian<T> {
+        let one = T::lift(1.0);
+        let mu = self.frame.mu();
+        let (sin_l, cos_l) = (self.true_longitude.sin(), self.true_longitude.cos());
+
+        let alpha_squared = self.h.powi(2) - self.k.powi(2);
+        let s_squared = one + self.h.powi(2) + self.k.powi(2);
+        let w = one + self.f * cos_l + self.g * sin_l;
+        let radius = self.semi_latus_rectum / w;
+        let speed_scale = (T::lift(mu) / self.semi_latus_rectum).sqrt();
+        let hk2 = self.h * self.k * 2.0;
+
+        let position = Vector3::new(
+            cos_l + alpha_squared * cos_l + hk2 * sin_l,
+            sin_l - alpha_squared * sin_l + hk2 * cos_l,
+            (self.h * sin_l - self.k * cos_l) * 2.0,
+        ) * (radius / s_squared);
+
+        let velocity = Vector3::new(
+            -(sin_l + alpha_squared * sin_l - hk2 * cos_l + self.g - hk2 * self.f
+                + alpha_squared * self.g),
+            -(-cos_l + alpha_squared * cos_l + hk2 * sin_l - self.f
+                + hk2 * self.g
+                + alpha_squared * self.f),
+            (self.h * cos_l + self.k * sin_l + self.f * self.h + self.g * self.k) * 2.0,
+        ) * (speed_scale / s_squared);
+
+        Cartesian::new(position, velocity, self.frame)
+    }
+
+    /// Elements from a Keplerian set, by the defining relations.
+    ///
+    /// # Panics
+    ///
+    /// If the orbit is within [`RETROGRADE_TOLERANCE`] of exactly retrograde.
+    pub fn from_keplerian(elements: &Keplerian<T>) -> Self {
+        assert!(
+            (core::f64::consts::PI - elements.inclination.re()).abs() > RETROGRADE_TOLERANCE,
+            "equinoctial elements in the prograde convention are singular at \
+             i = 180 degrees"
+        );
+        let longitude_of_periapsis = elements.aop + elements.raan;
+        let tangent = (elements.inclination * 0.5).tan();
+        Self {
+            semi_latus_rectum: elements.semi_major_axis
+                * (T::lift(1.0) - elements.eccentricity.powi(2)),
+            f: elements.eccentricity * longitude_of_periapsis.cos(),
+            g: elements.eccentricity * longitude_of_periapsis.sin(),
+            h: tangent * elements.raan.cos(),
+            k: tangent * elements.raan.sin(),
+            true_longitude: wrap_two_pi(longitude_of_periapsis + elements.true_anomaly),
+            frame: elements.frame,
+        }
+    }
+
+    /// The Keplerian set these describe.
+    ///
+    /// Where the Keplerian set is degenerate the usual substitutions apply and
+    /// [`Keplerian::degeneracy`] says so, exactly as
+    /// [`Keplerian::from_cartesian`] would report it.
+    pub fn to_keplerian(&self) -> Keplerian<T> {
+        self.to_cartesian().keplerian()
+    }
+
+    /// Semi-major axis, m. Negative on a hyperbolic orbit.
+    pub fn semi_major_axis(&self) -> T {
+        self.semi_latus_rectum / (T::lift(1.0) - self.eccentricity().powi(2))
+    }
+
+    /// Eccentricity, dimensionless.
+    pub fn eccentricity(&self) -> T {
+        (self.f.powi(2) + self.g.powi(2)).sqrt()
+    }
+
+    /// Inclination, rad.
+    pub fn inclination(&self) -> T {
+        (self.h.powi(2) + self.k.powi(2)).sqrt().atan() * 2.0
+    }
+
+    /// Right ascension of the ascending node, rad. Undefined for an equatorial
+    /// orbit, where `h` and `k` are both zero.
+    pub fn raan(&self) -> T {
+        self.k.atan2(self.h)
+    }
+
+    /// Orbital period, s.
+    pub fn period(&self) -> T {
+        kepler::period(self.semi_major_axis(), self.frame.mu())
+    }
+
+    /// Mean motion, rad/s.
+    pub fn mean_motion(&self) -> T {
+        kepler::mean_motion(self.semi_major_axis(), self.frame.mu())
+    }
+}
+
+impl<T: Scalar> Keplerian<T> {
+    /// These elements as a non-singular set.
+    pub fn equinoctial(&self) -> Equinoctial<T> {
+        Equinoctial::from_keplerian(self)
+    }
+}
+
+impl<T: Scalar> Cartesian<T> {
+    /// The equinoctial elements of this state.
+    pub fn equinoctial(&self) -> Equinoctial<T> {
+        Equinoctial::from_cartesian(self)
+    }
+}
+
+/// Starts a [`State`] from equinoctial elements, for any layout with a position
+/// and a velocity.
+pub trait FromEquinoctial<const D: usize>: HasPosition<D> + HasVelocity<D> {
+    /// Starts a [`State`] from equinoctial elements.
+    fn with_equinoctial<T: Scalar>(self, elements: Equinoctial<T>) -> State<T, D, Self> {
+        elements.to_cartesian().to_state(self)
+    }
+}
+
+impl<const D: usize, L: HasPosition<D> + HasVelocity<D>> FromEquinoctial<D> for L {}
+
+/// Whether a Keplerian set's degeneracy makes its angles untrustworthy.
+///
+/// A convenience for choosing a representation: if this is true, the equinoctial
+/// set is the one to work in.
+pub fn prefers_equinoctial(degeneracy: Degeneracy) -> bool {
+    degeneracy != Degeneracy::None
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::body::Body;
+    use crate::frames::Orientation;
+    use crate::state::Orbit;
+    use approx::assert_relative_eq;
+
+    const MU: f64 = Body::EARTH.mu;
+
+    fn inclined_elliptical() -> Cartesian<f64> {
+        Keplerian {
+            semi_major_axis: 2.0e7,
+            eccentricity: 0.3,
+            inclination: 45.0_f64.to_radians(),
+            raan: 30.0_f64.to_radians(),
+            aop: 60.0_f64.to_radians(),
+            true_anomaly: 120.0_f64.to_radians(),
+            frame: Frame::EARTH_J2000,
+            degeneracy: Degeneracy::None,
+        }
+        .to_cartesian()
+    }
+
+    fn circular_equatorial(radius: f64) -> Cartesian<f64> {
+        Cartesian::new(
+            Vector3::new(radius, 0.0, 0.0),
+            Vector3::new(0.0, kepler::circular_speed(radius, MU), 0.0),
+            Frame::EARTH_J2000,
+        )
+    }
+
+    fn leo() -> Cartesian<f64> {
+        Cartesian::new(
+            Vector3::new(
+                4.263868426884883e6,
+                5.146189057155391e6,
+                1.1310208421331816e6,
+            ),
+            Vector3::new(-5923.454461876975, 4496.802639690076, 1870.3893008991558),
+            Frame::EARTH_J2000,
+        )
+    }
+
+    /// The invariant that was wrong: the basis vectors' cross product has to be
+    /// the orbit normal that `h` and `k` describe, not some other unit vector.
+    #[test]
+    fn test_basis_is_right_handed_about_the_orbit_normal() {
+        for (inclination, raan) in [
+            (45.0_f64, 30.0_f64),
+            (0.0, 0.0),
+            (98.0, 200.0),
+            (170.0, 45.0),
+        ] {
+            let (inclination, raan) = (inclination.to_radians(), raan.to_radians());
+            let tangent: f64 = (inclination / 2.0).tan();
+            let (h, k) = (tangent * raan.cos(), tangent * raan.sin());
+
+            let expected = Vector3::new(
+                raan.sin() * inclination.sin(),
+                -raan.cos() * inclination.sin(),
+                inclination.cos(),
+            );
+            assert_relative_eq!(
+                Equinoctial::<f64>::basis_normal(h, k),
+                expected,
+                max_relative = 1e-12,
+                epsilon = 1e-14
+            );
+        }
+    }
+
+    #[test]
+    fn test_cartesian_round_trips() {
+        for state in [inclined_elliptical(), leo(), circular_equatorial(4.2164e7)] {
+            let back = state.equinoctial().to_cartesian();
+            // An equatorial orbit's z components are both about 1e-9 m, so a purely
+            // relative comparison on them is meaningless.
+            assert_relative_eq!(
+                back.position,
+                state.position,
+                max_relative = 1e-11,
+                epsilon = 1e-6
+            );
+            assert_relative_eq!(
+                back.velocity,
+                state.velocity,
+                max_relative = 1e-11,
+                epsilon = 1e-9
+            );
+            assert_eq!(back.frame, state.frame);
+        }
+    }
+
+    /// The direct Cartesian conversion has to agree with the route through
+    /// Keplerian, where Keplerian is trustworthy.
+    #[test]
+    fn test_direct_conversion_agrees_with_the_keplerian_route() {
+        let state = inclined_elliptical();
+        let direct = state.equinoctial();
+        let via_keplerian = state.keplerian().equinoctial();
+
+        assert_relative_eq!(
+            direct.semi_latus_rectum,
+            via_keplerian.semi_latus_rectum,
+            max_relative = 1e-11
+        );
+        assert_relative_eq!(
+            direct.f,
+            via_keplerian.f,
+            max_relative = 1e-11,
+            epsilon = 1e-12
+        );
+        assert_relative_eq!(
+            direct.g,
+            via_keplerian.g,
+            max_relative = 1e-11,
+            epsilon = 1e-12
+        );
+        assert_relative_eq!(
+            direct.h,
+            via_keplerian.h,
+            max_relative = 1e-11,
+            epsilon = 1e-12
+        );
+        assert_relative_eq!(
+            direct.k,
+            via_keplerian.k,
+            max_relative = 1e-11,
+            epsilon = 1e-12
+        );
+        assert_relative_eq!(
+            direct.true_longitude,
+            via_keplerian.true_longitude,
+            max_relative = 1e-11
+        );
+    }
+
+    /// The point of the module: a circular equatorial orbit is a degenerate
+    /// Keplerian set and an ordinary equinoctial one.
+    #[test]
+    fn test_geo_is_degenerate_in_keplerian_and_not_in_equinoctial() {
+        let state = circular_equatorial(4.2164e7);
+
+        assert_eq!(state.keplerian().degeneracy, Degeneracy::CircularEquatorial);
+        assert!(prefers_equinoctial(state.keplerian().degeneracy));
+
+        // Circular and equatorial is just the origin of the f-g-h-k space.
+        let elements = state.equinoctial();
+        assert_relative_eq!(elements.f, 0.0, epsilon = 1e-15);
+        assert_relative_eq!(elements.g, 0.0, epsilon = 1e-15);
+        assert_relative_eq!(elements.h, 0.0, epsilon = 1e-15);
+        assert_relative_eq!(elements.k, 0.0, epsilon = 1e-15);
+        assert_relative_eq!(elements.semi_latus_rectum, 4.2164e7, max_relative = 1e-12);
+        assert_relative_eq!(elements.true_longitude, 0.0, epsilon = 1e-12);
+    }
+
+    #[test]
+    fn test_derived_quantities_match_the_keplerian_set() {
+        let state = inclined_elliptical();
+        let keplerian = state.keplerian();
+        let equinoctial = state.equinoctial();
+
+        assert_relative_eq!(
+            equinoctial.semi_major_axis(),
+            keplerian.semi_major_axis,
+            max_relative = 1e-11
+        );
+        assert_relative_eq!(
+            equinoctial.eccentricity(),
+            keplerian.eccentricity,
+            max_relative = 1e-11
+        );
+        assert_relative_eq!(
+            equinoctial.inclination(),
+            keplerian.inclination,
+            max_relative = 1e-11
+        );
+        assert_relative_eq!(
+            equinoctial.raan(),
+            keplerian.raan,
+            max_relative = 1e-11,
+            epsilon = 1e-12
+        );
+        assert_relative_eq!(
+            equinoctial.period(),
+            keplerian.period(),
+            max_relative = 1e-11
+        );
+        assert_relative_eq!(
+            equinoctial.mean_motion(),
+            keplerian.mean_motion(),
+            max_relative = 1e-11
+        );
+    }
+
+    #[test]
+    fn test_keplerian_round_trips_through_equinoctial() {
+        let original = inclined_elliptical().keplerian();
+        let back = original.equinoctial().to_keplerian();
+
+        assert_relative_eq!(
+            back.semi_major_axis,
+            original.semi_major_axis,
+            max_relative = 1e-11
+        );
+        assert_relative_eq!(
+            back.eccentricity,
+            original.eccentricity,
+            max_relative = 1e-11
+        );
+        assert_relative_eq!(back.inclination, original.inclination, max_relative = 1e-11);
+        assert_relative_eq!(back.raan, original.raan, max_relative = 1e-11);
+        assert_relative_eq!(
+            back.aop,
+            original.aop,
+            max_relative = 1e-11,
+            epsilon = 1e-12
+        );
+        assert_relative_eq!(
+            back.true_anomaly,
+            original.true_anomaly,
+            max_relative = 1e-11
+        );
+    }
+
+    #[test]
+    fn test_state_from_equinoctial() {
+        let elements = inclined_elliptical().equinoctial();
+        let y = Orbit.with_equinoctial(elements);
+        let back = y.cartesian(Frame::EARTH_J2000).equinoctial();
+        assert_relative_eq!(back.f, elements.f, max_relative = 1e-11, epsilon = 1e-12);
+        assert_relative_eq!(
+            back.true_longitude,
+            elements.true_longitude,
+            max_relative = 1e-11
+        );
+    }
+
+    #[test]
+    #[should_panic(expected = "requires an inertial frame")]
+    fn test_a_rotating_frame_is_refused() {
+        let rotating = Cartesian::new(
+            leo().position,
+            leo().velocity,
+            Frame::EARTH_J2000.with_orientation(Orientation::BodyFixed),
+        );
+        let _ = rotating.equinoctial();
+    }
+
+    #[test]
+    #[should_panic(expected = "singular at i = 180")]
+    fn test_exactly_retrograde_is_refused() {
+        let radius = 7.0e6;
+        // Equatorial, going the wrong way round: i = 180 exactly.
+        let retrograde = Cartesian::new(
+            Vector3::new(radius, 0.0, 0.0),
+            Vector3::new(0.0, -kepler::circular_speed(radius, MU), 0.0),
+            Frame::EARTH_J2000,
+        );
+        let _ = retrograde.equinoctial();
+    }
+
+    /// A slightly retrograde orbit is fine — only the pole itself is singular.
+    #[test]
+    fn test_near_retrograde_is_fine() {
+        let radius = 7.0e6;
+        let speed = kepler::circular_speed(radius, MU);
+        // i = 170 degrees.
+        let inclination = 170.0_f64.to_radians();
+        let state = Cartesian::new(
+            Vector3::new(radius, 0.0, 0.0),
+            Vector3::new(0.0, speed * inclination.cos(), speed * inclination.sin()),
+            Frame::EARTH_J2000,
+        );
+
+        let elements = state.equinoctial();
+        assert_relative_eq!(elements.inclination(), inclination, max_relative = 1e-10);
+        let back = elements.to_cartesian();
+        assert_relative_eq!(back.position, state.position, max_relative = 1e-10);
+    }
+
+    /// Equinoctial elements are differentiable where Keplerian angles are not:
+    /// on a circular orbit, `d(aop)/d(state)` is meaningless, but `d(f)/d(state)`
+    /// is perfectly well behaved.
+    #[cfg(feature = "autodiff")]
+    #[test]
+    fn test_elements_differentiate_on_a_circular_orbit() {
+        use nalgebra::{SVector, Vector6};
+        use num_dual::DualSVec64;
+
+        let state = circular_equatorial(4.2164e7);
+        let packed = Vector6::new(
+            state.position[0],
+            state.position[1],
+            state.position[2],
+            state.velocity[0],
+            state.velocity[1],
+            state.velocity[2],
+        );
+
+        let (value, gradient): (f64, SVector<f64, 6>) = num_dual::gradient(
+            |y: Vector6<DualSVec64<6>>| {
+                Cartesian::new(
+                    y.fixed_rows::<3>(0).into_owned(),
+                    y.fixed_rows::<3>(3).into_owned(),
+                    Frame::EARTH_J2000,
+                )
+                .equinoctial()
+                .f
+            },
+            &packed,
+        );
+
+        assert_relative_eq!(value, 0.0, epsilon = 1e-15);
+        // Every partial is finite, which is the property that matters.
+        for partial in gradient.iter() {
+            assert!(partial.is_finite(), "gradient has a non-finite entry");
+        }
+        // And it agrees with central differences.
+        let f_of = |y: Vector6<f64>| {
+            Cartesian::new(
+                y.fixed_rows::<3>(0).into_owned(),
+                y.fixed_rows::<3>(3).into_owned(),
+                Frame::EARTH_J2000,
+            )
+            .equinoctial()
+            .f
+        };
+        for i in 0..6 {
+            let h = 1e-4 * packed[i].abs().max(1.0);
+            let mut plus = packed;
+            let mut minus = packed;
+            plus[i] += h;
+            minus[i] -= h;
+            let difference = (f_of(plus) - f_of(minus)) / (2.0 * h);
+            assert_relative_eq!(
+                gradient[i],
+                difference,
+                max_relative = 1e-5,
+                epsilon = 1e-16
+            );
+        }
+    }
+}

src/frames.rs

@@ -38,7 +38,7 @@ pub struct Frame {
 }
 
 impl Frame {
-    /// Earth-centered, J2000-oriented. The frame the study's scenarios work in.
+    /// Earth-centered, J2000-oriented.
     pub const EARTH_J2000: Self = Self {
         center: Body::EARTH,
         orientation: Orientation::J2000,

src/kepler.rs

@@ -123,7 +123,7 @@ fn angle_about<T: Scalar>(from: Vector3<T>, to: Vector3<T>, unit_normal: Vector3
 }
 
 /// Wraps an angle into `[0, 2pi)`.
-fn wrap_two_pi<T: Scalar>(angle: T) -> T {
+pub(crate) fn wrap_two_pi<T: Scalar>(angle: T) -> T {
     let turns = (angle.re() / TAU).floor();
     angle - TAU * turns
 }
@@ -144,8 +144,8 @@ fn wrap_two_pi<T: Scalar>(angle: T) -> T {
 /// - **circular and equatorial** — both are undefined, so `raan` and `aop` are
 ///   zero and `true_anomaly` carries the true longitude
 ///
-/// GEO hits the third case and near-circular LEO the first, so these are the
-/// common cases here rather than corner cases.
+/// A geostationary orbit hits the third case and a near-circular one the first,
+/// so these are common rather than corner cases.
 #[derive(Clone, Copy, Debug, PartialEq)]
 pub struct Keplerian<T: Scalar> {
     /// Semi-major axis, m.
@@ -483,8 +483,8 @@ mod tests {
         );
     }
 
-    /// GEO is circular and equatorial, the doubly-singular case, and it is one of
-    /// the study's scenarios rather than a corner case.
+    /// A geostationary orbit is circular and equatorial: the doubly-singular case,
+    /// and a thoroughly ordinary orbit.
     #[test]
     fn test_geo_is_reported_as_circular_equatorial() {
         let radius = 4.2164e7;

src/lib.rs

@@ -12,6 +12,7 @@
 pub mod body;
 pub mod cartesian;
 pub mod dynamics;
+pub mod equinoctial;
 pub mod frames;
 pub mod kepler;
 pub mod propagation;
@@ -25,6 +26,7 @@ pub mod prelude {
     pub use crate::body::Body;
     pub use crate::cartesian::{Cartesian, FromCartesian};
     pub use crate::dynamics::Dynamics::*;
+    pub use crate::equinoctial::{Equinoctial, FromEquinoctial};
     pub use crate::frames::{Frame, Orientation};
     pub use crate::kepler::{Degeneracy, Keplerian, TranslationalLayout};
     pub use crate::propagation::{Accuracy, Propagation};

src/state.rs

@@ -756,7 +756,8 @@ mod tests {
         assert_eq!(speed::<f64, 14, Full>(&full), 13.0);
     }
 
-    /// The `Orbit + k estimated parameters` case the filters will need.
+    /// The `Orbit + k estimated parameters` case, for a filter that estimates
+    /// more than the trajectory.
     #[test]
     fn test_user_defined_layout_with_estimated_parameters() {
         define_state! {