Added State Layouts and a few niceties like examples, git hooks, etc

.githooks/pre-commit

@@ -0,0 +1,44 @@
+#!/bin/sh
+# Enable with: git config core.hooksPath .githooks
+# Skip a single commit with: git commit --no-verify
+
+set -e
+
+echo "pre-commit: cargo fmt --check"
+cargo fmt --check
+
+echo "pre-commit: cargo clippy"
+cargo clippy --all-targets --all-features -- -D warnings
+
+echo "pre-commit: cargo test"
+cargo test --quiet
+
+echo "pre-commit: cargo doc"
+RUSTDOCFLAGS="-D warnings" cargo doc --no-deps --quiet
+
+echo "pre-commit: US spelling"
+staged=$(git diff --cached --name-only --diff-filter=ACM -- '*.rs' '*.md' '*.toml')
+if [ -n "$staged" ]; then
+    # shellcheck disable=SC2086
+    if git grep -wniE -e 'metre|metres|kilometre|kilometres' \
+        -e 'licence|defence|programme|catalogue|grey' \
+        -e 'manoeuvre|manoeuvres|manoeuvred|manoeuvring' \
+        -e 'normalise|normalised|normalises|normalising|normalisation' \
+        -e 'stabilise|stabilised|stabilises|stabilising|stabilisation' \
+        -e 'recognise|recognised|recognises|recognising' \
+        -e 'monomorphise|monomorphised|monomorphises|monomorphising' \
+        -e 'serialise|serialised|serialises|serialising|serialisation' \
+        -e 'behaviour|behaviours|colour|colours|centre|centres' \
+        -e 'analyse|analysed|analyses|analysing' \
+        -e 'organise|organised|organisation|optimise|optimised|optimisation' \
+        -e 'initialise|initialised|initialisation|utilise|utilised' \
+        -e 'modelling|modelled|labelled|cancelled|travelled|fuelled' \
+        -e 'favour|favours|honour|honours|neighbour|neighbours' \
+        --cached -- $staged; then
+        echo
+        echo "pre-commit: British spellings found above. This project uses US English."
+        exit 1
+    fi
+fi
+
+echo "pre-commit: ok"

examples/two_body.rs

@@ -0,0 +1,61 @@
+//! Two-body propagation over one orbital period.
+//!
+//! Run with `cargo run --example two_body`.
+
+use astrokit::state::{HasPosition, HasVelocity, Orbit};
+use nalgebra::{SVector, Vector6};
+use ordinary_diffeq::prelude::*;
+use std::f64::consts::PI;
+
+/// 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;
+
+/// The two-body derivative.
+fn two_body<T: Scalar, const D: usize, L: HasPosition<D> + HasVelocity<D>>(
+    _t: T,
+    y: SVector<T, D>,
+    &(mu,): &(T,),
+) -> SVector<T, D> {
+    let r = L::position(&y);
+    let mut dy = SVector::<T, D>::zeros();
+    L::set_position(&mut dy, L::velocity(&y));
+    L::set_velocity(&mut dy, -(r * mu) / r.norm().powi(3));
+    dy
+}
+
+fn main() {
+    let period = 2.0 * PI * (SMA.powi(3) / MU).sqrt();
+
+    let y0 = Vector6::new(
+        4.263868426884883e6,
+        5.146189057155391e6,
+        1.1310208421331816e6,
+        -5923.454461876975,
+        4496.802639690076,
+        1870.3893008991558,
+    );
+
+    let ode = ODE::new(&two_body::<f64, 6, Orbit>, 0.0, period, y0, (MU,));
+    let vern7 = Vern7::new().a_tol(1e-12).r_tol(1e-12);
+    let mut problem = Problem::new(ode, vern7, PIController::default());
+    let solution = problem.solve();
+
+    let end = solution.states.last().unwrap();
+    let position_error = (Orbit::position(end) - Orbit::position(&y0)).norm();
+    let velocity_error = (Orbit::velocity(end) - Orbit::velocity(&y0)).norm();
+
+    println!("period:         {period:.3} s");
+    println!("steps:          {}", solution.states.len() - 1);
+    println!("position error: {position_error:.3e} m");
+    println!("velocity error: {velocity_error:.3e} m/s");
+
+    // Dense output: the state at any time in the arc, not just at a step.
+    let half = solution.interpolate(period / 2.0);
+    println!(
+        "altitude at half a period: {:.3} km",
+        (Orbit::position(&half).norm() - 6.378137e6) / 1e3
+    );
+}

readme.md

@@ -15,7 +15,7 @@ The crate is still being built to a milestone plan:
   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, manoeuvres, covariance propagation,
+  density models, box-wing SRP and eclipses, maneuvers, covariance propagation,
   relative motion, TLE and CCSDS I/O.
 
 The tables under [Validation](#validation) and [Performance](#performance) are
@@ -25,12 +25,13 @@ the acceptance criteria.
 
 - Everything is SI
 - Time is typed, using `hifitime`
+- US English spelling
 
 ## Foundations
 
 These are depended on, not reimplemented.
 
-| Crate | Licence | Provides |
+| Crate | License | Provides |
 |---|---|---|
 | [`hifitime`](https://crates.io/crates/hifitime) 4.3 | MPL-2.0 | TAI/TT/ET/TDB/UTC/GPS time scales, leap seconds |
 | [`anise`](https://crates.io/crates/anise) 0.10 | MPL-2.0 | Reference frames, planetary ephemerides, SPICE-format kernels |
@@ -75,6 +76,18 @@ figures, so the references are measured here rather than cited.
 | Full-force LEO, 1 day | not yet measured | — | — | — |
 | Ensemble of N states, scaling vs. cores | not yet measured | — | — | — |
 
-## Licence
+## Development
+
+Enable the pre-commit hook once per clone:
+
+```sh
+git config core.hooksPath .githooks
+```
+
+It runs `cargo fmt --check`, `cargo clippy --all-targets --all-features -- -D
+warnings`, `cargo test`, `cargo doc` with warnings denied, and a US English
+spelling check. `git commit --no-verify` skips it.
+
+## License
 
 MIT. See [LICENSE](./LICENSE).

src/lib.rs

@@ -2,12 +2,14 @@
 //!
 //! # Conventions
 //!
-//! - Everything is SI: metres, seconds, kilograms, radians.
+//! - Everything is SI: meters, seconds, kilograms, radians.
 //! - Times are always typed using hifitime
 
 // Needed to build docs
 #![cfg_attr(feature = "autodiff", doc = include_str!("../readme.md"))]
 #![warn(missing_docs)]
 
+pub mod state;
+
 /// The types you want in scope to set up and run a propagation.
 pub mod prelude {}

src/state.rs

@@ -0,0 +1,494 @@
+//! State layouts: which quantities an integrated state vector holds, and where.
+//!
+//! The state dimension has to be a compile-time constant, because `num-dual`'s
+//! vector duals are only efficient at a `Const<N>`. Summing layer widths in type
+//! position would need `generic_const_exprs`, which is not going to stabilize,
+//! so [`define_state!`](crate::define_state) does the sum at macro-expansion time instead.
+//!
+//! Dynamics should be written against the capability traits rather than a
+//! concrete layout, so that they work for any state carrying what they need:
+//!
+//! ```
+//! use astrokit::state::{HasPosition, HasVelocity};
+//! use nalgebra::{SVector, Vector3};
+//! use ordinary_diffeq::scalar::Scalar;
+//!
+//! /// Specific angular momentum, in m^2/s.
+//! fn angular_momentum<T: Scalar, const D: usize, L: HasPosition<D> + HasVelocity<D>>(
+//!     y: &SVector<T, D>,
+//! ) -> Vector3<T> {
+//!     L::position(y).cross(&L::velocity(y))
+//! }
+//! ```
+
+use core::fmt::Debug;
+
+use nalgebra::{SVector, Vector3};
+use ordinary_diffeq::scalar::Scalar;
+
+/// One named span of a [`Layout`]. A component's offset is the sum of the widths
+/// preceding it.
+#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
+pub struct Component {
+    /// The component's name, as written at the layout's definition site.
+    pub name: &'static str,
+    /// How many state vector elements the component occupies.
+    pub width: usize,
+}
+
+/// A state layout of dimension `D`.
+///
+/// Layouts are zero-sized marker types, so the accessors monomorphize to a
+/// direct offset. Define one with [`define_state!`](crate::define_state) rather than by hand.
+pub trait Layout<const D: usize>: Copy + Clone + Debug + Default {
+    /// The layout's name, as written at its definition site
+    const NAME: &'static str;
+
+    /// The components, in index order.
+    const COMPONENTS: &'static [Component];
+}
+
+/// Defines a capability trait for a three-element component.
+macro_rules! vector_capability {
+    ($trait:ident, $index:ident, $get:ident, $set:ident, $what:expr, $unit:expr) => {
+        #[doc = concat!("Layouts carrying ", $what, ".")]
+        pub trait $trait<const D: usize>: Layout<D> {
+            #[doc = concat!("Index of the first of the three ", $what, " elements.")]
+            const $index: usize;
+
+            #[doc = concat!("The ", $what, ", in ", $unit, ".")]
+            ///
+            /// A derivative vector shares the state's layout, so this reads the
+            /// corresponding derivative when applied to one.
+            #[inline]
+            fn $get<T: Scalar>(y: &SVector<T, D>) -> Vector3<T> {
+                y.fixed_rows::<3>(Self::$index).into_owned()
+            }
+
+            #[doc = concat!("Writes the ", $what, ", in ", $unit, ".")]
+            #[inline]
+            fn $set<T: Scalar>(y: &mut SVector<T, D>, value: Vector3<T>) {
+                y.fixed_rows_mut::<3>(Self::$index).copy_from(&value);
+            }
+        }
+    };
+}
+
+/// Defines a capability trait for a one-element component.
+macro_rules! scalar_capability {
+    ($trait:ident, $index:ident, $get:ident, $set:ident, $what:expr, $unit:expr) => {
+        #[doc = concat!("Layouts carrying ", $what, ".")]
+        pub trait $trait<const D: usize>: Layout<D> {
+            #[doc = concat!("Index of the ", $what, " element.")]
+            const $index: usize;
+
+            #[doc = concat!("The ", $what, ", in ", $unit, ".")]
+            ///
+            /// A derivative vector shares the state's layout, so this reads the
+            /// corresponding derivative when applied to one.
+            #[inline]
+            fn $get<T: Scalar>(y: &SVector<T, D>) -> T {
+                y[Self::$index]
+            }
+
+            #[doc = concat!("Writes the ", $what, ", in ", $unit, ".")]
+            #[inline]
+            fn $set<T: Scalar>(y: &mut SVector<T, D>, value: T) {
+                y[Self::$index] = value;
+            }
+        }
+    };
+}
+
+vector_capability!(
+    HasPosition,
+    POSITION,
+    position,
+    set_position,
+    "an inertial position",
+    "meters"
+);
+vector_capability!(
+    HasVelocity,
+    VELOCITY,
+    velocity,
+    set_velocity,
+    "an inertial velocity",
+    "meters per second"
+);
+vector_capability!(
+    HasAttitude,
+    ATTITUDE,
+    attitude,
+    set_attitude,
+    "an attitude, as Modified Rodrigues Parameters",
+    "dimensionless units"
+);
+vector_capability!(
+    HasAngularRate,
+    ANGULAR_RATE,
+    angular_rate,
+    set_angular_rate,
+    "a body angular rate",
+    "radians per second"
+);
+scalar_capability!(HasMass, MASS, mass, set_mass, "a mass", "kilograms");
+scalar_capability!(
+    HasEnergy,
+    ENERGY,
+    energy,
+    set_energy,
+    "a power-system stored energy",
+    "joules"
+);
+
+/// Defines a state layout: a zero-sized type implementing [`Layout`] plus one
+/// capability trait per recognized component.
+///
+/// ```
+/// use astrokit::define_state;
+/// use astrokit::state::{HasMass, HasPosition, HasVelocity, Layout};
+///
+/// define_state! {
+///     /// Translational state plus mass and an estimated drag coefficient.
+///     pub struct OrbitMassCd {
+///         position,
+///         velocity,
+///         mass,
+///         scalar drag_coefficient,
+///     }
+/// }
+///
+/// assert_eq!(<OrbitMassCd as Layout<8>>::NAME, "OrbitMassCd");
+/// assert_eq!(<OrbitMassCd as HasPosition<8>>::POSITION, 0);
+/// assert_eq!(<OrbitMassCd as HasVelocity<8>>::VELOCITY, 3);
+/// assert_eq!(<OrbitMassCd as HasMass<8>>::MASS, 6);
+/// assert_eq!(OrbitMassCd::drag_coefficient, 7);
+/// ```
+///
+/// # Recognized components
+///
+/// | Component | Width | Capability |
+/// |---|---|---|
+/// | `position` | 3 | [`HasPosition`] |
+/// | `velocity` | 3 | [`HasVelocity`] |
+/// | `attitude` | 3 | [`HasAttitude`] |
+/// | `angular_rate` | 3 | [`HasAngularRate`] |
+/// | `mass` | 1 | [`HasMass`] |
+/// | `energy` | 1 | [`HasEnergy`] |
+/// | `scalar <name>` | 1 | none — an associated const `<name>` |
+/// | `vector3 <name>` | 3 | none — an associated const `<name>` |
+///
+/// `scalar` and `vector3` cover quantities this crate has no opinion about,
+/// which in practice means a filter's estimated parameters. They get an index
+/// constant and no capability trait.
+///
+/// Components are laid out in the order written; the dimension is their total
+/// width.
+#[macro_export]
+macro_rules! define_state {
+    (
+        $(#[$meta:meta])*
+        $vis:vis struct $name:ident {
+            $( $kind:tt $($label:ident)? ),+ $(,)?
+        }
+    ) => {
+        $(#[$meta])*
+        #[derive(Clone, Copy, Debug, Default, PartialEq, Eq, Hash)]
+        $vis struct $name;
+
+        impl $crate::state::Layout<{ 0 $( + $crate::__state_width!($kind) )+ }> for $name {
+            const NAME: &'static str = ::core::stringify!($name);
+            const COMPONENTS: &'static [$crate::state::Component] = &[
+                $( $crate::state::Component {
+                    name: $crate::__state_label!($kind $($label)?),
+                    width: $crate::__state_width!($kind),
+                } ),+
+            ];
+        }
+
+        $crate::__state_capabilities! {
+            @layout $name
+            @dim { 0 $( + $crate::__state_width!($kind) )+ }
+            @offset { 0 }
+            @rest [ $( $kind $($label)? ),+ ]
+        }
+    };
+}
+
+/// The width of one recognized component. Not public API.
+#[doc(hidden)]
+#[macro_export]
+macro_rules! __state_width {
+    (position) => {
+        3
+    };
+    (velocity) => {
+        3
+    };
+    (attitude) => {
+        3
+    };
+    (angular_rate) => {
+        3
+    };
+    (mass) => {
+        1
+    };
+    (energy) => {
+        1
+    };
+    (scalar) => {
+        1
+    };
+    (vector3) => {
+        3
+    };
+}
+
+/// The name of one component. Not public API.
+#[doc(hidden)]
+#[macro_export]
+macro_rules! __state_label {
+    (scalar $label:ident) => {
+        ::core::stringify!($label)
+    };
+    (vector3 $label:ident) => {
+        ::core::stringify!($label)
+    };
+    ($kind:tt) => {
+        ::core::stringify!($kind)
+    };
+}
+
+/// Walks the component list carrying the running offset, emitting one capability
+/// impl per component. Not public API.
+#[doc(hidden)]
+#[macro_export]
+macro_rules! __state_capabilities {
+    (@layout $name:ident @dim $dim:tt @offset $offset:tt @rest []) => {};
+
+    (
+        @layout $name:ident
+        @dim { $($dim:tt)* }
+        @offset { $($offset:tt)* }
+        @rest [ $kind:tt $($label:ident)? $(, $($rest:tt)*)? ]
+    ) => {
+        $crate::__state_one! { $name, { $($dim)* }, { $($offset)* }, $kind $($label)? }
+
+        $crate::__state_capabilities! {
+            @layout $name
+            @dim { $($dim)* }
+            @offset { $($offset)* + $crate::__state_width!($kind) }
+            @rest [ $( $($rest)* )? ]
+        }
+    };
+}
+
+/// Emits the capability impl for a single component. Not public API.
+#[doc(hidden)]
+#[macro_export]
+macro_rules! __state_one {
+    ($name:ident, { $($dim:tt)* }, { $($offset:tt)* }, position) => {
+        impl $crate::state::HasPosition<{ $($dim)* }> for $name {
+            const POSITION: usize = $($offset)*;
+        }
+    };
+    ($name:ident, { $($dim:tt)* }, { $($offset:tt)* }, velocity) => {
+        impl $crate::state::HasVelocity<{ $($dim)* }> for $name {
+            const VELOCITY: usize = $($offset)*;
+        }
+    };
+    ($name:ident, { $($dim:tt)* }, { $($offset:tt)* }, attitude) => {
+        impl $crate::state::HasAttitude<{ $($dim)* }> for $name {
+            const ATTITUDE: usize = $($offset)*;
+        }
+    };
+    ($name:ident, { $($dim:tt)* }, { $($offset:tt)* }, angular_rate) => {
+        impl $crate::state::HasAngularRate<{ $($dim)* }> for $name {
+            const ANGULAR_RATE: usize = $($offset)*;
+        }
+    };
+    ($name:ident, { $($dim:tt)* }, { $($offset:tt)* }, mass) => {
+        impl $crate::state::HasMass<{ $($dim)* }> for $name {
+            const MASS: usize = $($offset)*;
+        }
+    };
+    ($name:ident, { $($dim:tt)* }, { $($offset:tt)* }, energy) => {
+        impl $crate::state::HasEnergy<{ $($dim)* }> for $name {
+            const ENERGY: usize = $($offset)*;
+        }
+    };
+    ($name:ident, { $($dim:tt)* }, { $($offset:tt)* }, scalar $label:ident) => {
+        impl $name {
+            #[doc = concat!("Index of the `", ::core::stringify!($label), "` element.")]
+            #[allow(non_upper_case_globals)]
+            pub const $label: usize = $($offset)*;
+        }
+    };
+    ($name:ident, { $($dim:tt)* }, { $($offset:tt)* }, vector3 $label:ident) => {
+        impl $name {
+            #[doc = concat!(
+                "Index of the first of the three `", ::core::stringify!($label), "` elements."
+            )]
+            #[allow(non_upper_case_globals)]
+            pub const $label: usize = $($offset)*;
+        }
+    };
+}
+
+define_state! {
+    /// Position and velocity. All that orbit propagation without mass or
+    /// attitude needs.
+    pub struct Orbit { position, velocity }
+}
+
+define_state! {
+    /// Translational state plus mass, for a maneuver that expends propellant.
+    pub struct OrbitMass { position, velocity, mass }
+}
+
+define_state! {
+    /// Position, velocity, attitude, angular rate, mass and stored energy.
+    ///
+    /// The layout low-thrust modeling needs: thrust direction depends on the
+    /// attitude, the mass changes enough to matter, and available power bounds
+    /// the thrust.
+    pub struct Full {
+        position,
+        velocity,
+        attitude,
+        angular_rate,
+        mass,
+        energy,
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use nalgebra::{Vector1, Vector6};
+
+    #[test]
+    fn test_builtin_layout_dimensions_and_offsets() {
+        assert_eq!(<Orbit as Layout<6>>::NAME, "Orbit");
+        assert_eq!(<Orbit as HasPosition<6>>::POSITION, 0);
+        assert_eq!(<Orbit as HasVelocity<6>>::VELOCITY, 3);
+
+        assert_eq!(<OrbitMass as HasMass<7>>::MASS, 6);
+
+        assert_eq!(<Full as HasPosition<14>>::POSITION, 0);
+        assert_eq!(<Full as HasVelocity<14>>::VELOCITY, 3);
+        assert_eq!(<Full as HasAttitude<14>>::ATTITUDE, 6);
+        assert_eq!(<Full as HasAngularRate<14>>::ANGULAR_RATE, 9);
+        assert_eq!(<Full as HasMass<14>>::MASS, 12);
+        assert_eq!(<Full as HasEnergy<14>>::ENERGY, 13);
+    }
+
+    #[test]
+    fn test_components_describe_the_layout() {
+        let names: Vec<_> = <Full as Layout<14>>::COMPONENTS
+            .iter()
+            .map(|c| c.name)
+            .collect();
+        assert_eq!(
+            names,
+            [
+                "position",
+                "velocity",
+                "attitude",
+                "angular_rate",
+                "mass",
+                "energy"
+            ]
+        );
+        let total: usize = <Full as Layout<14>>::COMPONENTS
+            .iter()
+            .map(|c| c.width)
+            .sum();
+        assert_eq!(total, 14);
+    }
+
+    #[test]
+    fn test_accessors_round_trip() {
+        let mut y = Vector6::<f64>::zeros();
+        Orbit::set_position(&mut y, Vector3::new(1.0, 2.0, 3.0));
+        Orbit::set_velocity(&mut y, Vector3::new(4.0, 5.0, 6.0));
+
+        assert_eq!(y, Vector6::new(1.0, 2.0, 3.0, 4.0, 5.0, 6.0));
+        assert_eq!(Orbit::position(&y), Vector3::new(1.0, 2.0, 3.0));
+        assert_eq!(Orbit::velocity(&y), Vector3::new(4.0, 5.0, 6.0));
+    }
+
+    /// A function bounded on capabilities compiles once and works for every
+    /// layout that has them.
+    #[test]
+    fn test_capability_bounds_are_layout_agnostic() {
+        fn speed<T: Scalar, const D: usize, L: HasVelocity<D>>(y: &SVector<T, D>) -> T {
+            L::velocity(y).norm()
+        }
+
+        let mut orbit = Vector6::<f64>::zeros();
+        Orbit::set_velocity(&mut orbit, Vector3::new(3.0, 4.0, 0.0));
+        assert_eq!(speed::<f64, 6, Orbit>(&orbit), 5.0);
+
+        let mut full = SVector::<f64, 14>::zeros();
+        Full::set_velocity(&mut full, Vector3::new(0.0, 5.0, 12.0));
+        assert_eq!(speed::<f64, 14, Full>(&full), 13.0);
+    }
+
+    /// The `Orbit + k estimated parameters` case the filters will need.
+    #[test]
+    fn test_user_defined_layout_with_estimated_parameters() {
+        define_state! {
+            /// Orbit plus an estimated drag coefficient and a range bias.
+            struct OrbitCdBias {
+                position,
+                velocity,
+                scalar drag_coefficient,
+                scalar range_bias,
+            }
+        }
+
+        assert_eq!(<OrbitCdBias as Layout<8>>::NAME, "OrbitCdBias");
+        assert_eq!(<OrbitCdBias as HasPosition<8>>::POSITION, 0);
+        assert_eq!(OrbitCdBias::drag_coefficient, 6);
+        assert_eq!(OrbitCdBias::range_bias, 7);
+    }
+
+    #[test]
+    fn test_layout_of_a_single_component() {
+        define_state! {
+            /// Degenerate but legal: one scalar.
+            struct MassOnly { mass }
+        }
+
+        assert_eq!(<MassOnly as Layout<1>>::NAME, "MassOnly");
+        assert_eq!(<MassOnly as HasMass<1>>::MASS, 0);
+
+        let mut y = Vector1::<f64>::zeros();
+        MassOnly::set_mass(&mut y, 42.0);
+        assert_eq!(MassOnly::mass(&y), 42.0);
+    }
+
+    /// Accessors have to work at a dual scalar, or the layout is unusable in the
+    /// differentiated path.
+    #[cfg(feature = "autodiff")]
+    #[test]
+    fn test_accessors_at_a_dual_scalar() {
+        use approx::assert_relative_eq;
+        use num_dual::Dual64;
+
+        let mut y = Vector6::<Dual64>::zeros();
+        let seeded = Dual64::from_re(3.0).derivative();
+        Orbit::set_position(
+            &mut y,
+            Vector3::new(seeded, Dual64::lift(4.0), Dual64::lift(0.0)),
+        );
+
+        // d|r|/dx = x/|r| = 3/5
+        let range = Orbit::position(&y).norm();
+        assert_relative_eq!(range.re, 5.0, max_relative = 1e-15);
+        assert_relative_eq!(range.eps, 0.6, max_relative = 1e-15);
+    }
+}