Added full cartesian support with (the beginning of) frames

benches/kepler.rs

@@ -6,31 +6,36 @@
 
 use std::hint::black_box;
 
+use astrokit::body::Body;
+use astrokit::cartesian::Cartesian;
+use astrokit::frames::Frame;
 use astrokit::kepler::{self, Degeneracy, Keplerian};
 use criterion::{criterion_group, criterion_main, Criterion};
 use nalgebra::Vector3;
 
-const MU: f64 = 3.986004415e14;
+const MU: f64 = Body::EARTH.mu;
 
 /// Inclined, near-circular LEO. Just outside the degeneracy tolerance, so it
 /// takes the general branch.
-fn leo() -> (Vector3<f64>, Vector3<f64>) {
+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,
     )
 }
 
 /// GEO: circular and equatorial, so it takes the doubly-degenerate branch.
-fn geo() -> (Vector3<f64>, Vector3<f64>) {
+fn geo() -> Cartesian<f64> {
     let radius = 4.2164e7;
-    (
+    Cartesian::new(
         Vector3::new(radius, 0.0, 0.0),
         Vector3::new(0.0, kepler::circular_speed(radius, MU), 0.0),
+        Frame::EARTH_J2000,
     )
 }
 
@@ -73,14 +78,14 @@ fn bench_anomalies(c: &mut Criterion) {
 fn bench_conversions(c: &mut Criterion) {
     let mut group = c.benchmark_group("kepler/conversions");
 
-    let (r, v) = leo();
+    let leo = leo();
     group.bench_function("from_cartesian/leo", |b| {
-        b.iter(|| Keplerian::from_cartesian(black_box(r), black_box(v), black_box(MU)))
+        b.iter(|| Keplerian::from_cartesian(&black_box(leo)))
     });
 
-    let (r_geo, v_geo) = geo();
+    let geo = geo();
     group.bench_function("from_cartesian/geo", |b| {
-        b.iter(|| Keplerian::from_cartesian(black_box(r_geo), black_box(v_geo), black_box(MU)))
+        b.iter(|| Keplerian::from_cartesian(&black_box(geo)))
     });
 
     let elements = Keplerian {
@@ -90,30 +95,28 @@ fn bench_conversions(c: &mut Criterion) {
         raan: 0.5,
         aop: 1.0,
         true_anomaly: 2.0,
+        frame: Frame::EARTH_J2000,
         degeneracy: Degeneracy::None,
     };
     group.bench_function("to_cartesian", |b| {
-        b.iter(|| black_box(elements).to_cartesian(black_box(MU)))
+        b.iter(|| black_box(elements).to_cartesian())
     });
 
     group.bench_function("round_trip", |b| {
-        b.iter(|| {
+        b.iter(|| Keplerian::from_cartesian(&black_box(leo)).to_cartesian())
-            let elements = Keplerian::from_cartesian(black_box(r), black_box(v), MU);
-            black_box(elements.to_cartesian(MU))
-        })
     });
     group.finish();
 }
 
 fn bench_invariants(c: &mut Criterion) {
     let mut group = c.benchmark_group("kepler/invariants");
-    let (r, v) = leo();
+    let leo = leo();
 
     group.bench_function("specific_energy", |b| {
-        b.iter(|| kepler::specific_energy(black_box(r), black_box(v), black_box(MU)))
+        b.iter(|| black_box(leo).specific_energy())
     });
     group.bench_function("specific_angular_momentum", |b| {
-        b.iter(|| kepler::specific_angular_momentum(black_box(r), black_box(v)))
+        b.iter(|| black_box(leo).specific_angular_momentum())
     });
     group.finish();
 }

benches/propagation.rs

@@ -36,7 +36,7 @@ fn velocity() -> Vector3<f64> {
 }
 
 fn two_body() -> Vec<Dynamics> {
-    vec![TwoBody { mu: MU }]
+    vec![TwoBody { body: Body::EARTH }]
 }
 
 /// One derivative evaluation, per layout. The layout should cost nothing.
@@ -45,7 +45,7 @@ fn bench_derivative(c: &mut Criterion) {
 
     let assembly = Assembly::new(
         dynamics::build::<f64, 6, Orbit>(&two_body()),
-        Context::new(epoch()),
+        Context::new(epoch(), Frame::EARTH_J2000),
     );
     let y = Orbit.with_position(position()).with_velocity(velocity());
     group.bench_function("two_body/Orbit", |b| {
@@ -54,7 +54,7 @@ fn bench_derivative(c: &mut Criterion) {
 
     let assembly = Assembly::new(
         dynamics::build::<f64, 7, OrbitMass>(&two_body()),
-        Context::new(epoch()),
+        Context::new(epoch(), Frame::EARTH_J2000),
     );
     let y = OrbitMass
         .with_position(position())
@@ -66,7 +66,7 @@ fn bench_derivative(c: &mut Criterion) {
 
     let assembly = Assembly::new(
         dynamics::build::<f64, 14, Full>(&two_body()),
-        Context::new(epoch()),
+        Context::new(epoch(), Frame::EARTH_J2000),
     );
     let y = Full
         .with_position(position())

examples/two_body.rs

@@ -2,46 +2,38 @@
 //!
 //! Run with `cargo run --example two_body`.
 
-use astrokit::kepler;
 use astrokit::prelude::*;
 use nalgebra::Vector3;
 
-/// Earth's gravitational parameter, m^3/s^2.
-const MU: f64 = 3.986004415e14;
-
-/// Semi-major axis of the test orbit, m.
-const SMA: f64 = 6.7781363e6;
-
-/// Earth's equatorial radius, m. Comes from a body once bodies land.
-const EARTH_RADIUS: f64 = 6.378137e6;
-
 fn main() {
     let epoch = Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0);
-    let period = Duration::from_seconds(kepler::period(SMA, MU));
+    let frame = Frame::EARTH_J2000;
 
-    let y0 = Orbit
+    let initial = Cartesian::new(
-        .with_position(Vector3::new(
+        Vector3::new(
             4.263868426884883e6,
             5.146189057155391e6,
             1.1310208421331816e6,
-        ))
+        ),
-        .with_velocity(Vector3::new(
+        Vector3::new(-5923.454461876975, 4496.802639690076, 1870.3893008991558),
-            -5923.454461876975,
+        frame,
-            4496.802639690076,
+    );
-            1870.3893008991558,
-        ));
 
-    let forces = vec![TwoBody { mu: MU }];
+    let elements = initial.keplerian();
+    let period = Duration::from_seconds(initial.period());
+    let forces = vec![TwoBody { body: Body::EARTH }];
 
-    let mut propagation = Propagation::new(epoch, y0, &forces)
+    let mut propagation = Propagation::from_cartesian(epoch, initial, &forces)
         .duration(period)
         .accuracy(Accuracy::High);
     let solution = propagation.run();
 
-    let end = Orbit.state(*solution.states.last().unwrap());
+    let end = Orbit
+        .state(*solution.states.last().unwrap())
+        .cartesian(frame);
 
-    let elements = y0.keplerian(MU);
     println!("epoch:          {epoch}");
+    println!("body:           {}", initial.frame.center.name());
     println!("period:         {period}");
     println!(
         "elements:       a = {:.3} km, e = {:.6}, i = {:.3} deg ({:?})",
@@ -50,24 +42,30 @@ fn main() {
         elements.inclination.to_degrees(),
         elements.degeneracy
     );
+    println!(
+        "altitude:       {:.3} km",
+        initial.spherical_altitude() / 1e3
+    );
     println!("steps:          {}", solution.states.len() - 1);
     println!(
         "position error: {:.3e} m",
-        (end.position() - y0.position()).norm()
+        (end.position - initial.position).norm()
     );
     println!(
         "velocity error: {:.3e} m/s",
-        (end.velocity() - y0.velocity()).norm()
+        (end.velocity - initial.velocity).norm()
+    );
+    println!(
+        "energy drift:   {:.3e} J/kg",
+        end.specific_energy() - initial.specific_energy()
     );
 
     // Dense output: the state at any time in the arc, not only at a step.
-    let half = Orbit.state(solution.interpolate(period.to_seconds() / 2.0));
+    let half = Orbit
+        .state(solution.interpolate(period.to_seconds() / 2.0))
+        .cartesian(frame);
     println!(
-        "altitude at half a period: {:.3} km",
+        "half a period:  {:.3} km altitude",
-        (half.position().norm() - EARTH_RADIUS) / 1e3
+        half.spherical_altitude() / 1e3
-    );
-    println!(
-        "specific energy drift: {:.3e} J/kg",
-        end.specific_energy(MU) - y0.specific_energy(MU)
     );
 }

src/body.rs

@@ -0,0 +1,128 @@
+//! Celestial bodies.
+
+use serde::{Deserialize, Serialize};
+
+/// A celestial body: the gravitational parameter and shape that force models and
+/// element conversions need.
+///
+/// Carried by a [`Frame`](crate::frames::Frame) and by the dynamics that need a
+/// central body, so that no signature takes a bare `mu`. A propagation about the
+/// Moon is then a matter of saying so rather than of working around an assumption.
+#[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)]
+pub struct Body {
+    /// NAIF integer id, for ANISE ephemeris and orientation lookups.
+    pub naif_id: i32,
+    /// Gravitational parameter, m^3/s^2.
+    pub mu: f64,
+    /// Equatorial radius, m.
+    pub equatorial_radius: f64,
+    /// Flattening, dimensionless. Zero for a sphere.
+    pub flattening: f64,
+}
+
+impl Body {
+    /// Earth.
+    ///
+    /// `mu` is EGM96's. **It must agree with whatever gravity field a
+    /// propagation uses**: the point-mass term and the spherical harmonics are two
+    /// halves of one expansion, and taking `mu` from one model and the harmonics
+    /// from another is an inconsistency. Once a gravity
+    /// field is loaded it should supply its own `mu` and override this.
+    ///
+    /// Radius and flattening are WGS84.
+    pub const EARTH: Self = Self {
+        naif_id: 399,
+        mu: 3.986004415e14,
+        equatorial_radius: 6378137.0,
+        flattening: 1.0 / 298.257223563,
+    };
+
+    /// The Moon. `mu` and radius from DE440; the figure is very nearly spherical.
+    pub const MOON: Self = Self {
+        naif_id: 301,
+        mu: 4.902800118e12,
+        equatorial_radius: 1737400.0,
+        flattening: 0.0012,
+    };
+
+    /// The Sun, treated as spherical.
+    pub const SUN: Self = Self {
+        naif_id: 10,
+        mu: 1.32712440018e20,
+        equatorial_radius: 6.957e8,
+        flattening: 0.0,
+    };
+
+    /// A spherical body, for a test case or a body whose figure does not matter.
+    pub const fn spherical(naif_id: i32, mu: f64, radius: f64) -> Self {
+        Self {
+            naif_id,
+            mu,
+            equatorial_radius: radius,
+            flattening: 0.0,
+        }
+    }
+
+    /// Polar radius, m.
+    pub fn polar_radius(&self) -> f64 {
+        self.equatorial_radius * (1.0 - self.flattening)
+    }
+
+    /// A name for the recognized bodies, or a rendering of the NAIF id.
+    ///
+    /// Diagnostics only. Provenance records the fields, which is the thing that
+    /// actually pins a run down.
+    pub fn name(&self) -> String {
+        match self.naif_id {
+            10 => "Sun".to_string(),
+            199 => "Mercury".to_string(),
+            299 => "Venus".to_string(),
+            301 => "Moon".to_string(),
+            399 => "Earth".to_string(),
+            499 => "Mars".to_string(),
+            599 => "Jupiter".to_string(),
+            699 => "Saturn".to_string(),
+            799 => "Uranus".to_string(),
+            899 => "Neptune".to_string(),
+            id => format!("NAIF {id}"),
+        }
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use approx::assert_relative_eq;
+
+    #[test]
+    fn test_earth_matches_wgs84() {
+        assert_eq!(Body::EARTH.equatorial_radius, 6378137.0);
+        // WGS84's polar radius, to the millimeter.
+        assert_relative_eq!(Body::EARTH.polar_radius(), 6356752.314245, epsilon = 1e-3);
+        assert_eq!(Body::EARTH.name(), "Earth");
+    }
+
+    #[test]
+    fn test_spherical_body_has_equal_radii() {
+        let body = Body::spherical(9999, 1.0e14, 1.0e6);
+        assert_eq!(body.polar_radius(), body.equatorial_radius);
+        assert_eq!(body.name(), "NAIF 9999");
+    }
+
+    /// Bodies go into a run's provenance, so they have to round-trip.
+    #[test]
+    fn test_body_round_trips_through_json() {
+        let json = serde_json::to_string(&Body::EARTH).unwrap();
+        assert_eq!(serde_json::from_str::<Body>(&json).unwrap(), Body::EARTH);
+    }
+
+    /// The gravitational parameters should agree with ANISE's, which come from
+    /// JPL kernels. ANISE works in km^3/s^2.
+    #[test]
+    fn test_gravitational_parameters_are_the_expected_magnitude() {
+        // Earth: EGM96 against the DE440 value, which differ in the 9th digit.
+        assert_relative_eq!(Body::EARTH.mu * 1e-9, 398600.4415, max_relative = 1e-10);
+        assert_relative_eq!(Body::MOON.mu * 1e-9, 4902.800118, max_relative = 1e-10);
+        assert_relative_eq!(Body::SUN.mu * 1e-9, 132712440018.0, max_relative = 1e-10);
+    }
+}

src/cartesian.rs

@@ -0,0 +1,225 @@
+//! A translational state as position, velocity and frame.
+//!
+//! This is the interchange representation. An integrated
+//! [`State`] is layout-shaped and may carry mass, attitude or
+//! estimated parameters; a `Cartesian` is the translational part of it together
+//! with the frame it is expressed in, which is what every element conversion needs
+//! and what no bare pair of vectors can supply.
+
+use nalgebra::Vector3;
+use ordinary_diffeq::scalar::Scalar;
+
+use crate::frames::Frame;
+use crate::kepler::{self, Keplerian};
+use crate::state::{HasPosition, HasVelocity, State};
+
+/// A position and velocity in a reference frame.
+#[derive(Clone, Copy, Debug, PartialEq)]
+pub struct Cartesian<T: Scalar> {
+    /// Position, m.
+    pub position: Vector3<T>,
+    /// Velocity, m/s.
+    pub velocity: Vector3<T>,
+    /// The frame both are expressed in.
+    pub frame: Frame,
+}
+
+impl<T: Scalar> Cartesian<T> {
+    /// A state from a position and velocity in `frame`.
+    pub fn new(position: Vector3<T>, velocity: Vector3<T>, frame: Frame) -> Self {
+        Self {
+            position,
+            velocity,
+            frame,
+        }
+    }
+
+    /// Distance from the frame's origin, m.
+    pub fn radius(&self) -> T {
+        self.position.norm()
+    }
+
+    /// Speed, m/s.
+    pub fn speed(&self) -> T {
+        self.velocity.norm()
+    }
+
+    /// Height above the central body's equatorial radius, m.
+    ///
+    /// Spherical, not geodetic: it subtracts the equatorial radius rather than
+    /// solving for the distance to the reference ellipsoid. For the Earth the two
+    /// differ by up to about 21 km at the poles. Geodetic altitude arrives with
+    /// the frame transformations, which is also where it belongs, since it is only
+    /// meaningful in a body-fixed frame.
+    pub fn spherical_altitude(&self) -> T {
+        self.radius() - self.frame.center.equatorial_radius
+    }
+
+    /// Specific orbital energy, J/kg.
+    pub fn specific_energy(&self) -> T {
+        kepler::specific_energy(self.position, self.velocity, self.frame.mu())
+    }
+
+    /// Specific angular momentum, m^2/s.
+    pub fn specific_angular_momentum(&self) -> Vector3<T> {
+        kepler::specific_angular_momentum(self.position, self.velocity)
+    }
+
+    /// The Keplerian elements of this state.
+    ///
+    /// # Panics
+    ///
+    /// If the frame is not inertial. Orbital elements describe inertial motion, and
+    /// a velocity in a rotating frame is not the inertial velocity — converting one
+    /// anyway would return a plausible-looking wrong answer.
+    pub fn keplerian(&self) -> Keplerian<T> {
+        Keplerian::from_cartesian(self)
+    }
+
+    /// Orbital period, s.
+    pub fn period(&self) -> T {
+        kepler::period(self.keplerian().semi_major_axis, self.frame.mu())
+    }
+
+    /// This state in a layout, with everything the layout also carries left zero.
+    pub fn to_state<const D: usize, L>(&self, layout: L) -> State<T, D, L>
+    where
+        L: HasPosition<D> + HasVelocity<D>,
+    {
+        layout
+            .with_position(self.position)
+            .with_velocity(self.velocity)
+    }
+
+    /// Asserts the frame is inertial, for the operations that require it.
+    pub(crate) fn require_inertial(&self, what: &str) {
+        assert!(
+            self.frame.is_inertial(),
+            "{what} requires an inertial frame, but this state is in {:?} about {}",
+            self.frame.orientation,
+            self.frame.center.name()
+        );
+    }
+}
+
+impl<T: Scalar, const D: usize, L: HasPosition<D> + HasVelocity<D>> State<T, D, L> {
+    /// The translational part of this state, in `frame`.
+    ///
+    /// A layout carries no frame of its own — a propagation holds that — so it has
+    /// to be supplied here.
+    pub fn cartesian(&self, frame: Frame) -> Cartesian<T> {
+        Cartesian::new(self.position(), self.velocity(), frame)
+    }
+}
+
+/// Layouts that can be filled from a [`Cartesian`].
+pub trait FromCartesian<const D: usize>: HasPosition<D> + HasVelocity<D> {
+    /// Starts a [`State`] from a Cartesian state, discarding its frame.
+    fn with_cartesian<T: Scalar>(self, state: &Cartesian<T>) -> State<T, D, Self> {
+        state.to_state(self)
+    }
+}
+
+impl<const D: usize, L: HasPosition<D> + HasVelocity<D>> FromCartesian<D> for L {}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::body::Body;
+    use crate::frames::Orientation;
+    use crate::state::{Full, Orbit, OrbitMass};
+    use approx::assert_relative_eq;
+
+    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,
+        )
+    }
+
+    #[test]
+    fn test_radius_speed_and_altitude() {
+        let state = leo();
+        assert_relative_eq!(state.radius(), 6.778135e6, max_relative = 1e-5);
+        assert_relative_eq!(state.speed(), 7668.5, max_relative = 1e-3);
+        // A 400 km circular orbit, measured off the equatorial radius.
+        assert_relative_eq!(state.spherical_altitude(), 399998.0, max_relative = 1e-3);
+    }
+
+    #[test]
+    fn test_invariants_read_mu_from_the_frame() {
+        let state = leo();
+        assert_relative_eq!(
+            state.specific_energy(),
+            kepler::specific_energy(state.position, state.velocity, Body::EARTH.mu),
+            max_relative = 1e-15
+        );
+        assert_eq!(
+            state.specific_angular_momentum(),
+            state.position.cross(&state.velocity)
+        );
+    }
+
+    /// The same Cartesian state about a different body gives different elements,
+    /// because the frame supplies `mu`. This is the whole reason for the type.
+    #[test]
+    fn test_the_frame_changes_the_elements() {
+        let about_earth = leo();
+        let about_moon = Cartesian::new(
+            about_earth.position,
+            about_earth.velocity,
+            Frame::MOON_J2000,
+        );
+
+        let earth = about_earth.keplerian();
+        let moon = about_moon.keplerian();
+
+        // Same geometry, different central mass, so a different semi-major axis.
+        assert_relative_eq!(earth.inclination, moon.inclination, max_relative = 1e-12);
+        assert!(earth.semi_major_axis != moon.semi_major_axis);
+        // Far too fast for the Moon to hold on to: hyperbolic, so sma is negative.
+        assert!(moon.semi_major_axis < 0.0);
+    }
+
+    #[test]
+    fn test_round_trips_through_a_layout() {
+        let state = leo();
+
+        for_each_layout(&state);
+
+        fn for_each_layout(state: &Cartesian<f64>) {
+            let on_orbit = state.to_state(Orbit);
+            assert_eq!(on_orbit.position(), state.position);
+            assert_eq!(on_orbit.cartesian(state.frame), *state);
+
+            let on_mass = Orbit.with_cartesian(state);
+            assert_eq!(on_mass.velocity(), state.velocity);
+
+            // A wider layout takes the translational part and leaves the rest zero.
+            let on_full = state.to_state(Full);
+            assert_eq!(on_full.position(), state.position);
+            assert_eq!(on_full.mass(), 0.0);
+            assert_eq!(on_full.cartesian(state.frame), *state);
+
+            let on_orbit_mass = OrbitMass.with_cartesian(state).with_mass(500.0);
+            assert_eq!(on_orbit_mass.mass(), 500.0);
+            assert_eq!(on_orbit_mass.cartesian(state.frame), *state);
+        }
+    }
+
+    #[test]
+    #[should_panic(expected = "requires an inertial frame")]
+    fn test_elements_from_a_rotating_frame_are_refused() {
+        let rotating = Cartesian::new(
+            leo().position,
+            leo().velocity,
+            Frame::EARTH_J2000.with_orientation(Orientation::BodyFixed),
+        );
+        let _ = rotating.keplerian();
+    }
+}

src/dynamics/mod.rs

@@ -16,7 +16,7 @@
 //! use astrokit::dynamics::{build, Dynamics};
 //! use astrokit::state::Orbit;
 //!
-//! let spec = vec![Dynamics::TwoBody { mu: 3.986004415e14 }];
+//! let spec = vec![Dynamics::TwoBody { body: astrokit::body::Body::EARTH }];
 //!
 //! // The same spec, at two scalars.
 //! let plain = build::<f64, 6, Orbit>(&spec);
@@ -29,8 +29,9 @@
 //! hidden behind a mode flag.
 //!
 //! A side benefit: because a model is constructed at the propagated scalar, it
-//! stores its constants already converted. [`models::TwoBody`] holds `mu` as a
+//! stores its constants already converted. [`models::TwoBody`] holds its body's
-//! `T`, so `T::lift` happens once at build time and never in the derivative.
+//! `mu` as a `T`, so `T::lift` happens once at build time and never in the
+//! derivative.
 //!
 //! # Force or acceleration
 //!
@@ -47,6 +48,8 @@ use ordinary_diffeq::scalar::Scalar;
 use serde::{Deserialize, Serialize};
 use std::ops::{Add, AddAssign};
 
+use crate::body::Body;
+use crate::frames::Frame;
 use crate::state::{HasPosition, HasVelocity, Layout, State};
 
 pub mod models;
@@ -124,17 +127,20 @@ impl<T: Scalar> Add for Contribution<T> {
 }
 
 /// What a [`DynamicsModel`] needs that is not in the state vector.
-#[derive(Clone, Copy, Debug, PartialEq, Eq)]
+#[derive(Clone, Copy, Debug, PartialEq)]
 pub struct Context {
     /// The epoch the propagation started from. The `t` a model receives is
     /// seconds after this.
     pub epoch: Epoch,
+    /// The frame the state is expressed in, and so the central body a model can
+    /// assume without being told.
+    pub frame: Frame,
 }
 
 impl Context {
-    /// A context anchored at `epoch`.
+    /// A context anchored at `epoch`, in `frame`.
-    pub fn new(epoch: Epoch) -> Self {
+    pub fn new(epoch: Epoch, frame: Frame) -> Self {
-        Self { epoch }
+        Self { epoch, frame }
     }
 
     /// The epoch `t` seconds after the start of the propagation.
@@ -179,10 +185,13 @@ pub trait DynamicsModel<T: Scalar, const D: usize, L: Layout<D>> {
 #[non_exhaustive]
 #[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)]
 pub enum Dynamics {
-    /// Point-mass gravity of a central body.
+    /// Point-mass gravity of a body.
+    ///
+    /// Named explicitly rather than taken from the propagation's frame, because an
+    /// n-body composition needs several and only one of them is the center.
     TwoBody {
-        /// Gravitational parameter, m^3/s^2.
+        /// The attracting body.
-        mu: f64,
+        body: Body,
     },
 }
 
@@ -194,7 +203,7 @@ impl Dynamics {
         L: HasPosition<D> + HasVelocity<D> + 'static,
     {
         match *self {
-            Self::TwoBody { mu } => Box::new(models::TwoBody::new(mu)),
+            Self::TwoBody { body } => Box::new(models::TwoBody::new(body)),
         }
     }
 }
@@ -246,7 +255,7 @@ mod tests {
 
     #[test]
     fn test_spec_round_trips_through_json() {
-        let spec = vec![Dynamics::TwoBody { mu: 3.986004415e14 }];
+        let spec = vec![Dynamics::TwoBody { body: Body::EARTH }];
         let json = serde_json::to_string(&spec).unwrap();
         let back: Vec<Dynamics> = serde_json::from_str(&json).unwrap();
         assert_eq!(spec, back);
@@ -256,7 +265,7 @@ mod tests {
     /// built at two different scalars.
     #[test]
     fn test_one_spec_builds_at_two_scalars() {
-        let spec = vec![Dynamics::TwoBody { mu: 3.986004415e14 }];
+        let spec = vec![Dynamics::TwoBody { body: Body::EARTH }];
 
         let plain = build::<f64, 6, Orbit>(&spec);
         assert_eq!(plain.len(), 1);
@@ -272,7 +281,7 @@ mod tests {
     #[test]
     fn test_context_advances_an_epoch() {
         let epoch = Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0);
-        let context = Context::new(epoch);
+        let context = Context::new(epoch, Frame::EARTH_J2000);
         assert_eq!(context.epoch_at(0.0), epoch);
         assert_eq!(
             context.epoch_at(3600.0),
@@ -283,11 +292,14 @@ mod tests {
     /// A model bounded on a capability works for any layout that has it.
     #[test]
     fn test_two_body_acceleration_is_layout_agnostic() {
-        let context = Context::new(Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0));
+        let context = Context::new(
+            Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0),
+            Frame::EARTH_J2000,
+        );
         let mu = 3.986004415e14;
         let r = Vector3::new(7.0e6, 0.0, 0.0);
 
-        let model = models::TwoBody::<f64>::new(mu);
+        let model = models::TwoBody::<f64>::new(Body::EARTH);
         let on_orbit = Orbit.with_position(r);
         let on_full = Full.with_position(r).with_mass(500.0);
 

src/dynamics/models/two_body.rs

@@ -3,22 +3,31 @@
 use ordinary_diffeq::scalar::Scalar;
 
 use super::super::{Context, Contribution, DynamicsModel};
+use crate::body::Body;
 use crate::state::{HasPosition, State};
 
-/// Point-mass gravity of a single central body, at the origin of the state's
+/// Point-mass gravity of a single body, at the origin of the state's frame.
-/// frame.
 ///
-/// `mu` is stored in the propagated scalar, converted once at construction, so
+/// The body's `mu` is stored in the propagated scalar, converted once at
-/// the derivative does no lifting.
+/// construction, so the derivative does no lifting.
 #[derive(Clone, Copy, Debug, PartialEq)]
 pub struct TwoBody<T: Scalar> {
+    body: Body,
     mu: T,
 }
 
 impl<T: Scalar> TwoBody<T> {
-    /// Point-mass gravity with gravitational parameter `mu`, in m^3/s^2.
+    /// Point-mass gravity of `body`.
-    pub fn new(mu: f64) -> Self {
+    pub fn new(body: Body) -> Self {
-        Self { mu: T::lift(mu) }
+        Self {
+            body,
+            mu: T::lift(body.mu),
+        }
+    }
+
+    /// The attracting body.
+    pub fn body(&self) -> Body {
+        self.body
     }
 
     /// The gravitational parameter, m^3/s^2.
@@ -42,6 +51,7 @@ impl<T: Scalar, const D: usize, L: HasPosition<D>> DynamicsModel<T, D, L> for Tw
 #[cfg(test)]
 mod tests {
     use super::*;
+    use crate::frames::Frame;
     use crate::state::Orbit;
     use approx::assert_relative_eq;
     use hifitime::Epoch;
@@ -51,8 +61,11 @@ mod tests {
 
     #[test]
     fn test_acceleration_points_at_the_origin() {
-        let context = Context::new(Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0));
+        let context = Context::new(
-        let model = TwoBody::<f64>::new(MU);
+            Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0),
+            Frame::EARTH_J2000,
+        );
+        let model = TwoBody::<f64>::new(Body::EARTH);
         let y = Orbit.with_position(Vector3::new(3.0e6, 4.0e6, 0.0));
 
         let c = model.contribution(0.0, &y, &context);

src/frames.rs

@@ -0,0 +1,136 @@
+//! Reference frames: where a state's origin is and how its axes are oriented.
+//!
+//! This module carries the frame *type*. Transformations between orientations —
+//! J2000 to ITRF, topocentric, geodetic site coordinates — are a separate piece of
+//! work and will land here on top of ANISE.
+
+use serde::{Deserialize, Serialize};
+
+use crate::body::Body;
+
+/// How a frame's axes are oriented.
+///
+/// Provisional and deliberately coarse. The distinctions that matter for
+/// sub-arcsecond work — ICRF against EME2000, ITRF against a body-fixed IAU frame
+/// — arrive with the transformations, which is the point at which getting them
+/// wrong would show up. Marked `#[non_exhaustive]` so refining it is not a
+/// breaking change.
+#[non_exhaustive]
+#[derive(Clone, Copy, Debug, PartialEq, Eq, Serialize, Deserialize)]
+pub enum Orientation {
+    /// Inertial, aligned with the J2000 mean equator and equinox.
+    J2000,
+    /// Rotating with the central body.
+    BodyFixed,
+}
+
+/// A reference frame: a central body and an orientation.
+///
+/// Held by every [`Cartesian`](crate::cartesian::Cartesian) and every element set,
+/// which is what lets conversions read `mu` from the state rather than take it as
+/// an argument.
+#[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)]
+pub struct Frame {
+    /// The body at the origin.
+    pub center: Body,
+    /// How the axes are oriented.
+    pub orientation: Orientation,
+}
+
+impl Frame {
+    /// Earth-centered, J2000-oriented. The frame the study's scenarios work in.
+    pub const EARTH_J2000: Self = Self {
+        center: Body::EARTH,
+        orientation: Orientation::J2000,
+    };
+
+    /// Earth-centered, rotating with the Earth.
+    pub const EARTH_FIXED: Self = Self {
+        center: Body::EARTH,
+        orientation: Orientation::BodyFixed,
+    };
+
+    /// Moon-centered, J2000-oriented.
+    pub const MOON_J2000: Self = Self {
+        center: Body::MOON,
+        orientation: Orientation::J2000,
+    };
+
+    /// A J2000-oriented frame centered on `center`.
+    pub const fn inertial(center: Body) -> Self {
+        Self {
+            center,
+            orientation: Orientation::J2000,
+        }
+    }
+
+    /// A frame rotating with `center`.
+    pub const fn body_fixed(center: Body) -> Self {
+        Self {
+            center,
+            orientation: Orientation::BodyFixed,
+        }
+    }
+
+    /// The central body's gravitational parameter, m^3/s^2.
+    pub fn mu(&self) -> f64 {
+        self.center.mu
+    }
+
+    /// Whether orbital motion in this frame is unaccelerated by the frame itself.
+    ///
+    /// Element conversions and the two-body relations are only meaningful in an
+    /// inertial frame; in a rotating one the velocity is not the inertial velocity.
+    pub fn is_inertial(&self) -> bool {
+        matches!(self.orientation, Orientation::J2000)
+    }
+
+    /// The same frame with a different orientation.
+    pub fn with_orientation(&self, orientation: Orientation) -> Self {
+        Self {
+            orientation,
+            ..*self
+        }
+    }
+}
+
+impl Default for Frame {
+    fn default() -> Self {
+        Self::EARTH_J2000
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn test_frame_carries_the_central_bodys_mu() {
+        assert_eq!(Frame::EARTH_J2000.mu(), Body::EARTH.mu);
+        assert_eq!(Frame::MOON_J2000.mu(), Body::MOON.mu);
+    }
+
+    #[test]
+    fn test_only_inertial_frames_report_inertial() {
+        assert!(Frame::EARTH_J2000.is_inertial());
+        assert!(!Frame::EARTH_FIXED.is_inertial());
+        assert!(Frame::inertial(Body::SUN).is_inertial());
+        assert!(!Frame::body_fixed(Body::SUN).is_inertial());
+    }
+
+    #[test]
+    fn test_orientation_can_be_swapped_keeping_the_center() {
+        let fixed = Frame::EARTH_J2000.with_orientation(Orientation::BodyFixed);
+        assert_eq!(fixed, Frame::EARTH_FIXED);
+        assert_eq!(fixed.center, Body::EARTH);
+    }
+
+    #[test]
+    fn test_frame_round_trips_through_json() {
+        let json = serde_json::to_string(&Frame::MOON_J2000).unwrap();
+        assert_eq!(
+            serde_json::from_str::<Frame>(&json).unwrap(),
+            Frame::MOON_J2000
+        );
+    }
+}

src/kepler.rs

@@ -15,6 +15,8 @@ use core::f64::consts::TAU;
 use nalgebra::Vector3;
 use ordinary_diffeq::scalar::Scalar;
 
+use crate::cartesian::Cartesian;
+use crate::frames::Frame;
 use crate::state::{HasPosition, HasVelocity, State};
 
 /// How close to zero an eccentricity or inclination has to be before
@@ -158,6 +160,9 @@ pub struct Keplerian<T: Scalar> {
     pub aop: T,
     /// True anomaly, rad.
     pub true_anomaly: T,
+    /// The frame these describe motion in. Supplies `mu`, so no method here takes
+    /// one.
+    pub frame: Frame,
     /// Which elements, if any, are not what they are named.
     pub degeneracy: Degeneracy,
 }
@@ -176,11 +181,17 @@ pub enum Degeneracy {
 }
 
 impl<T: Scalar> Keplerian<T> {
-    /// Elements from an inertial position and velocity.
+    /// Elements from a Cartesian state.
     ///
-    /// Vallado's RV2COE, with the three singular cases handled as
+    /// Vallado's RV2COE, with the three singular cases handled as [`Keplerian`]
-    /// [`Keplerian`] documents.
+    /// documents and the angles measured by `atan2` rather than `acos`.
-    pub fn from_cartesian(position: Vector3<T>, velocity: Vector3<T>, mu: f64) -> Self {
+    ///
+    /// # Panics
+    ///
+    /// If the state's frame is not inertial.
+    pub fn from_cartesian(state: &Cartesian<T>) -> Self {
+        state.require_inertial("converting to Keplerian elements");
+        let (position, velocity, mu) = (state.position, state.velocity, state.frame.mu());
         let zero = T::lift(0.0);
         let r = position.norm();
 
@@ -265,17 +276,19 @@ impl<T: Scalar> Keplerian<T> {
             raan,
             aop,
             true_anomaly,
+            frame: state.frame,
             degeneracy,
         }
     }
 
-    /// An inertial position and velocity, in meters and meters per second.
+    /// A Cartesian state, in the frame these elements carry.
     ///
-    /// Inverts [`Keplerian::from_cartesian`] for a non-degenerate set. A
+    /// Inverts [`Keplerian::from_cartesian`] for a non-degenerate set. A degenerate
-    /// degenerate set round-trips too, because the substituted angles are zero
+    /// set round-trips too, because the substituted angles are zero where they are
-    /// where they are undefined and the perifocal construction then puts the
+    /// undefined and the perifocal construction then puts the argument of latitude
-    /// argument of latitude or true longitude in the right place.
+    /// or true longitude in the right place.
-    pub fn to_cartesian(&self, mu: f64) -> (Vector3<T>, Vector3<T>) {
+    pub fn to_cartesian(&self) -> Cartesian<T> {
+        let mu = self.frame.mu();
         let one = T::lift(1.0);
         let semi_latus = self.semi_major_axis * (one - self.eccentricity.powi(2));
         let radius = semi_latus / (one + self.eccentricity * self.true_anomaly.cos());
@@ -293,9 +306,10 @@ impl<T: Scalar> Keplerian<T> {
             T::lift(0.0),
         );
 
-        (
+        Cartesian::new(
             self.perifocal_to_inertial(in_plane),
             self.perifocal_to_inertial(velocity_in_plane),
+            self.frame,
         )
     }
 
@@ -319,13 +333,13 @@ impl<T: Scalar> Keplerian<T> {
     }
 
     /// Orbital period, in seconds.
-    pub fn period(&self, mu: f64) -> T {
+    pub fn period(&self) -> T {
-        period(self.semi_major_axis, mu)
+        period(self.semi_major_axis, self.frame.mu())
     }
 
     /// Mean motion, in radians per second.
-    pub fn mean_motion(&self, mu: f64) -> T {
+    pub fn mean_motion(&self) -> T {
-        mean_motion(self.semi_major_axis, mu)
+        mean_motion(self.semi_major_axis, self.frame.mu())
     }
 
     /// Eccentric anomaly, in radians.
@@ -355,31 +369,20 @@ impl<T: Scalar> Keplerian<T> {
 /// these without doing anything.
 pub trait TranslationalLayout<const D: usize>: HasPosition<D> + HasVelocity<D> {
     /// Starts a [`State`] from Keplerian elements.
-    fn with_keplerian<T: Scalar>(self, elements: Keplerian<T>, mu: f64) -> State<T, D, Self> {
+    fn with_keplerian<T: Scalar>(self, elements: Keplerian<T>) -> State<T, D, Self> {
-        let (position, velocity) = elements.to_cartesian(mu);
+        elements.to_cartesian().to_state(self)
-        self.with_position(position).with_velocity(velocity)
     }
 }
 
 impl<const D: usize, L: HasPosition<D> + HasVelocity<D>> TranslationalLayout<D> for L {}
 
 impl<T: Scalar, const D: usize, L: HasPosition<D> + HasVelocity<D>> State<T, D, L> {
-    /// The Keplerian elements of this state, about a body of parameter `mu`.
+    /// The Keplerian elements of this state, in `frame`.
-    pub fn keplerian(&self, mu: f64) -> Keplerian<T> {
+    pub fn keplerian(&self, frame: Frame) -> Keplerian<T> {
-        Keplerian::from_cartesian(self.position(), self.velocity(), mu)
+        self.cartesian(frame).keplerian()
     }
 
-    /// Orbital period, in seconds.
+    /// Specific angular momentum, in m^2/s. Needs no frame.
-    pub fn period(&self, mu: f64) -> T {
-        period(self.keplerian(mu).semi_major_axis, mu)
-    }
-
-    /// Specific orbital energy, in J/kg.
-    pub fn specific_energy(&self, mu: f64) -> T {
-        specific_energy(self.position(), self.velocity(), mu)
-    }
-
-    /// Specific angular momentum, in m^2/s.
     pub fn specific_angular_momentum(&self) -> Vector3<T> {
         specific_angular_momentum(self.position(), self.velocity())
     }
@@ -388,21 +391,33 @@ impl<T: Scalar, const D: usize, L: HasPosition<D> + HasVelocity<D>> State<T, D,
 #[cfg(test)]
 mod tests {
     use super::*;
+    use crate::cartesian::FromCartesian;
     use crate::state::Orbit;
     use approx::assert_relative_eq;
     use core::f64::consts::PI;
 
-    const MU: f64 = 3.986004415e14;
+    /// Earth's gravitational parameter, for the free functions that still take one.
+    const MU: f64 = crate::body::Body::EARTH.mu;
 
     /// The state the propagation tests use, an inclined near-circular LEO.
-    fn leo() -> (Vector3<f64>, Vector3<f64>) {
+    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,
+        )
+    }
+
+    /// A circular orbit of the given radius in the equatorial plane.
+    fn circular_equatorial(radius: f64) -> Cartesian<f64> {
+        Cartesian::new(
+            Vector3::new(radius, 0.0, 0.0),
+            Vector3::new(0.0, circular_speed(radius, MU), 0.0),
+            Frame::EARTH_J2000,
         )
     }
 
@@ -418,23 +433,18 @@ mod tests {
     #[test]
     fn test_circular_speed_gives_a_closed_orbit() {
         let radius = 7.0e6;
-        let speed = circular_speed(radius, MU);
+        let elements = circular_equatorial(radius).keplerian();
-        let elements = Keplerian::from_cartesian(
-            Vector3::new(radius, 0.0, 0.0),
-            Vector3::new(0.0, speed, 0.0),
-            MU,
-        );
         assert_relative_eq!(elements.semi_major_axis, radius, max_relative = 1e-12);
         assert!(elements.eccentricity < 1e-15);
     }
 
     #[test]
     fn test_cartesian_round_trips_through_elements() {
-        let (r, v) = leo();
+        let state = leo();
-        let elements = Keplerian::from_cartesian(r, v, MU);
+        let back = state.keplerian().to_cartesian();
-        let (r2, v2) = elements.to_cartesian(MU);
+        assert_relative_eq!(back.position, state.position, max_relative = 1e-12);
-        assert_relative_eq!(r2, r, max_relative = 1e-12);
+        assert_relative_eq!(back.velocity, state.velocity, max_relative = 1e-12);
-        assert_relative_eq!(v2, v, max_relative = 1e-12);
+        assert_eq!(back.frame, state.frame);
     }
 
     /// An inclined elliptical orbit, where no element is degenerate.
@@ -447,10 +457,10 @@ mod tests {
             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,
         };
-        let (r, v) = elements.to_cartesian(MU);
+        let back = elements.to_cartesian().keplerian();
-        let back = Keplerian::from_cartesian(r, v, MU);
 
         assert_eq!(back.degeneracy, Degeneracy::None);
         assert_relative_eq!(
@@ -478,11 +488,7 @@ mod tests {
     #[test]
     fn test_geo_is_reported_as_circular_equatorial() {
         let radius = 4.2164e7;
-        let elements = Keplerian::from_cartesian(
+        let elements = circular_equatorial(radius).keplerian();
-            Vector3::new(radius, 0.0, 0.0),
-            Vector3::new(0.0, circular_speed(radius, MU), 0.0),
-            MU,
-        );
 
         assert_eq!(elements.degeneracy, Degeneracy::CircularEquatorial);
         assert_relative_eq!(elements.semi_major_axis, radius, max_relative = 1e-12);
@@ -492,11 +498,12 @@ mod tests {
         assert_relative_eq!(elements.true_anomaly, 0.0, epsilon = 1e-12);
 
         // A quarter turn later the true longitude is 90 degrees.
-        let quarter = Keplerian::from_cartesian(
+        let quarter = Cartesian::new(
             Vector3::new(0.0, radius, 0.0),
             Vector3::new(-circular_speed(radius, MU), 0.0, 0.0),
-            MU,
+            Frame::EARTH_J2000,
-        );
+        )
+        .keplerian();
         assert_relative_eq!(quarter.true_anomaly, PI / 2.0, max_relative = 1e-12);
     }
 
@@ -506,11 +513,12 @@ mod tests {
         let radius = 7.0e6;
         let speed = circular_speed(radius, MU);
         // Ascending node along +x, so at +x the argument of latitude is zero.
-        let elements = Keplerian::from_cartesian(
+        let elements = Cartesian::new(
             Vector3::new(radius, 0.0, 0.0),
             Vector3::new(0.0, speed * 0.5, speed * 3.0_f64.sqrt() / 2.0),
-            MU,
+            Frame::EARTH_J2000,
-        );
+        )
+        .keplerian();
 
         assert_eq!(elements.degeneracy, Degeneracy::Circular);
         assert_relative_eq!(
@@ -554,6 +562,7 @@ mod tests {
             raan: 0.0,
             aop: 0.0,
             true_anomaly: 0.0,
+            frame: Frame::EARTH_J2000,
             degeneracy: Degeneracy::None,
         };
         assert_relative_eq!(elements.periapsis(), 1.4e7, max_relative = 1e-12);
@@ -565,21 +574,19 @@ mod tests {
 
     #[test]
     fn test_state_accessors() {
-        let (r, v) = leo();
+        let state = leo();
-        let y = Orbit.with_position(r).with_velocity(v);
+        let y = Orbit.with_cartesian(&state);
 
-        let elements = y.keplerian(MU);
+        let elements = y.keplerian(state.frame);
         assert_relative_eq!(
-            y.period(MU),
+            elements.period(),
             period(elements.semi_major_axis, MU),
             max_relative = 1e-12
         );
-        assert_relative_eq!(
+        assert_eq!(
-            y.specific_energy(MU),
+            y.specific_angular_momentum(),
-            specific_energy(r, v, MU),
+            state.position.cross(&state.velocity)
-            max_relative = 1e-12
         );
-        assert_eq!(y.specific_angular_momentum(), r.cross(&v));
     }
 
     #[test]
@@ -591,10 +598,11 @@ mod tests {
             raan: 0.3,
             aop: 0.6,
             true_anomaly: 1.2,
+            frame: Frame::EARTH_J2000,
             degeneracy: Degeneracy::None,
         };
-        let y = Orbit.with_keplerian(elements, MU);
+        let y = Orbit.with_keplerian(elements);
-        let back = y.keplerian(MU);
+        let back = y.keplerian(Frame::EARTH_J2000);
         assert_relative_eq!(
             back.semi_major_axis,
             elements.semi_major_axis,
@@ -616,8 +624,9 @@ mod tests {
         use anise::prelude::Orbit as AniseOrbit;
         use hifitime::Epoch;
 
-        let (r, v) = leo();
+        let state = leo();
-        let ours = Keplerian::from_cartesian(r, v, MU);
+        let (r, v) = (state.position, state.velocity);
+        let ours = state.keplerian();
 
         let frame = EARTH_J2000.with_mu_km3_s2(MU * 1e-9);
         let theirs = AniseOrbit::new(
@@ -662,7 +671,7 @@ mod tests {
             max_relative = 1e-9
         );
         assert_relative_eq!(
-            period(ours.semi_major_axis, MU),
+            ours.period(),
             theirs.period().unwrap().to_seconds(),
             max_relative = 1e-9
         );
@@ -676,7 +685,8 @@ mod tests {
         use nalgebra::{SVector, Vector6};
         use num_dual::DualSVec64;
 
-        let (r, v) = leo();
+        let leo = leo();
+        let (r, v) = (leo.position, leo.velocity);
         let state = Vector6::new(r[0], r[1], r[2], v[0], v[1], v[2]);
 
         // d(sma)/d(state), by automatic differentiation.
@@ -684,25 +694,24 @@ mod tests {
             |y: Vector6<DualSVec64<6>>| {
                 let position = y.fixed_rows::<3>(0).into_owned();
                 let velocity = y.fixed_rows::<3>(3).into_owned();
-                Keplerian::from_cartesian(position, velocity, MU).semi_major_axis
+                Cartesian::new(position, velocity, Frame::EARTH_J2000)
+                    .keplerian()
+                    .semi_major_axis
             },
             &state,
         );
 
-        assert_relative_eq!(
+        assert_relative_eq!(sma, leo.keplerian().semi_major_axis, max_relative = 1e-12);
-            sma,
-            Keplerian::from_cartesian(r, v, MU).semi_major_axis,
-            max_relative = 1e-12
-        );
 
         // Cross-check the gradient against central differences, which agree only
         // to about the square root of machine epsilon.
         let sma_of = |y: Vector6<f64>| {
-            Keplerian::from_cartesian(
+            Cartesian::new(
                 y.fixed_rows::<3>(0).into_owned(),
                 y.fixed_rows::<3>(3).into_owned(),
-                MU,
+                Frame::EARTH_J2000,
             )
+            .keplerian()
             .semi_major_axis
         };
         for i in 0..6 {

src/lib.rs

@@ -9,7 +9,10 @@
 #![cfg_attr(feature = "autodiff", doc = include_str!("../readme.md"))]
 #![warn(missing_docs)]
 
+pub mod body;
+pub mod cartesian;
 pub mod dynamics;
+pub mod frames;
 pub mod kepler;
 pub mod propagation;
 pub mod state;
@@ -19,7 +22,10 @@ pub mod prelude {
     pub use crate::dynamics::{Context, Contribution, Dynamics, DynamicsModel};
     // The variants, so a script writes `TwoBody { mu }` rather than
     // `Dynamics::TwoBody { mu }`. Note this puts every model's name in scope.
+    pub use crate::body::Body;
+    pub use crate::cartesian::{Cartesian, FromCartesian};
     pub use crate::dynamics::Dynamics::*;
+    pub use crate::frames::{Frame, Orientation};
     pub use crate::kepler::{Degeneracy, Keplerian, TranslationalLayout};
     pub use crate::propagation::{Accuracy, Propagation};
     pub use crate::state::{

src/propagation.rs

@@ -6,8 +6,10 @@ use ordinary_diffeq::integrator::bs3::BS3;
 use ordinary_diffeq::prelude::*;
 use serde::{Deserialize, Serialize};
 
+use crate::cartesian::Cartesian;
 use crate::dynamics::{self, Context, Contribution, Dynamics, DynamicsModel};
-use crate::state::{HasPosition, HasVelocity, Layout, State};
+use crate::frames::Frame;
+use crate::state::{HasPosition, HasVelocity, Layout, Orbit, State};
 
 /// How precisely the equations are solved.
 ///
@@ -189,6 +191,7 @@ pub struct Propagation<
 > {
     epoch: Epoch,
     initial_state: State<T, D, L>,
+    frame: Frame,
     duration: Option<Duration>,
     specs: Vec<Dynamics>,
     integrator: S,
@@ -210,6 +213,7 @@ where
         Self {
             epoch,
             initial_state,
+            frame: Frame::default(),
             duration: None,
             specs: dynamics.to_vec(),
             integrator: Vern7::new().a_tol(1e-10).r_tol(1e-10),
@@ -219,12 +223,40 @@ where
     }
 }
 
+impl<T> Propagation<T, 6, Orbit, Vern7<6>>
+where
+    T: Scalar + 'static,
+    Vern7<6>: Integrator<6>,
+{
+    /// A propagation of a translational state, in that state's own frame.
+    ///
+    /// The shorthand for the common case: a [`Cartesian`] already carries its
+    /// frame, so neither the layout nor the frame has to be named.
+    ///
+    /// Use [`Propagation::new`] with a wider layout when the state also carries
+    /// mass, attitude or estimated parameters.
+    pub fn from_cartesian(epoch: Epoch, state: Cartesian<T>, dynamics: &[Dynamics]) -> Self {
+        let frame = state.frame;
+        Self::new(epoch, state.to_state(Orbit), dynamics).frame(frame)
+    }
+}
+
 impl<T, const D: usize, L, S> Propagation<T, D, L, S>
 where
     T: Scalar + 'static,
     L: HasPosition<D> + HasVelocity<D> + 'static,
     S: Tolerances<D> + Copy,
 {
+    /// The frame the state is expressed in.
+    ///
+    /// Defaults to [`Frame::EARTH_J2000`]. The default lives here rather than in
+    /// the force models, so a propagation about another body is a matter of saying
+    /// so rather than of overriding an assumption in several places.
+    pub fn frame(mut self, frame: Frame) -> Self {
+        self.frame = frame;
+        self
+    }
+
     /// How long to propagate for.
     pub fn duration(mut self, duration: Duration) -> Self {
         self.duration = Some(duration);
@@ -282,6 +314,7 @@ where
         Propagation {
             epoch: self.epoch,
             initial_state: self.initial_state,
+            frame: self.frame,
             duration: self.duration,
             specs: self.specs,
             integrator,
@@ -312,7 +345,7 @@ where
 
         let assembly = Assembly {
             models: dynamics::build::<T, D, L>(&self.specs),
-            context: Context::new(self.epoch),
+            context: Context::new(self.epoch, self.frame),
         };
 
         let ode = ODE::new(
@@ -331,6 +364,7 @@ where
 #[cfg(test)]
 mod tests {
     use super::*;
+    use crate::body::Body;
     use crate::state::{Full, Orbit, OrbitMass};
     use approx::assert_relative_eq;
     use nalgebra::{Vector3, Vector6};
@@ -344,7 +378,7 @@ mod tests {
     }
 
     fn two_body() -> Vec<Dynamics> {
-        vec![Dynamics::TwoBody { mu: MU }]
+        vec![Dynamics::TwoBody { body: Body::EARTH }]
     }
 
     fn leo() -> State<f64, 6, Orbit> {
@@ -565,7 +599,7 @@ mod tests {
             vec![Box::new(ConstantForce {
                 force: Vector3::new(10.0, 0.0, 0.0),
             })],
-            Context::new(epoch()),
+            Context::new(epoch(), Frame::EARTH_J2000),
         );
 
         let y = OrbitMass
@@ -587,9 +621,9 @@ mod tests {
                 Box::new(ConstantForce {
                     force: Vector3::new(10.0, 0.0, 0.0),
                 }),
-                Dynamics::TwoBody { mu: MU }.build(),
+                Dynamics::TwoBody { body: Body::EARTH }.build(),
             ],
-            Context::new(epoch()),
+            Context::new(epoch(), Frame::EARTH_J2000),
         );
 
         let y = OrbitMass
@@ -612,7 +646,7 @@ mod tests {
             vec![Box::new(ConstantForce {
                 force: Vector3::new(10.0, 0.0, 0.0),
             })],
-            Context::new(epoch()),
+            Context::new(epoch(), Frame::EARTH_J2000),
         );
         let _ = state_derivative(0.0, leo().vector, &assembly);
     }