diff --git a/.githooks/pre-commit b/.githooks/pre-commit new file mode 100755 index 0000000..0771bd7 --- /dev/null +++ b/.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" diff --git a/examples/two_body.rs b/examples/two_body.rs new file mode 100644 index 0000000..b8273ef --- /dev/null +++ b/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 + HasVelocity>( + _t: T, + y: SVector, + &(mu,): &(T,), +) -> SVector { + let r = L::position(&y); + let mut dy = SVector::::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::, 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 + ); +} diff --git a/readme.md b/readme.md index 32b3080..7df96cf 100644 --- a/readme.md +++ b/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). diff --git a/src/lib.rs b/src/lib.rs index c8ec037..801314c 100644 --- a/src/lib.rs +++ b/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 {} diff --git a/src/state.rs b/src/state.rs new file mode 100644 index 0000000..76d09ca --- /dev/null +++ b/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`. 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 + HasVelocity>( +//! y: &SVector, +//! ) -> Vector3 { +//! 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: 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: Layout { + #[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(y: &SVector) -> Vector3 { + y.fixed_rows::<3>(Self::$index).into_owned() + } + + #[doc = concat!("Writes the ", $what, ", in ", $unit, ".")] + #[inline] + fn $set(y: &mut SVector, value: Vector3) { + 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: Layout { + #[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(y: &SVector) -> T { + y[Self::$index] + } + + #[doc = concat!("Writes the ", $what, ", in ", $unit, ".")] + #[inline] + fn $set(y: &mut SVector, 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!(>::NAME, "OrbitMassCd"); +/// assert_eq!(>::POSITION, 0); +/// assert_eq!(>::VELOCITY, 3); +/// assert_eq!(>::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 ` | 1 | none — an associated const `` | +/// | `vector3 ` | 3 | none — an associated const `` | +/// +/// `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!(>::NAME, "Orbit"); + assert_eq!(>::POSITION, 0); + assert_eq!(>::VELOCITY, 3); + + assert_eq!(>::MASS, 6); + + assert_eq!(>::POSITION, 0); + assert_eq!(>::VELOCITY, 3); + assert_eq!(>::ATTITUDE, 6); + assert_eq!(>::ANGULAR_RATE, 9); + assert_eq!(>::MASS, 12); + assert_eq!(>::ENERGY, 13); + } + + #[test] + fn test_components_describe_the_layout() { + let names: Vec<_> = >::COMPONENTS + .iter() + .map(|c| c.name) + .collect(); + assert_eq!( + names, + [ + "position", + "velocity", + "attitude", + "angular_rate", + "mass", + "energy" + ] + ); + let total: usize = >::COMPONENTS + .iter() + .map(|c| c.width) + .sum(); + assert_eq!(total, 14); + } + + #[test] + fn test_accessors_round_trip() { + let mut y = Vector6::::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>(y: &SVector) -> T { + L::velocity(y).norm() + } + + let mut orbit = Vector6::::zeros(); + Orbit::set_velocity(&mut orbit, Vector3::new(3.0, 4.0, 0.0)); + assert_eq!(speed::(&orbit), 5.0); + + let mut full = SVector::::zeros(); + Full::set_velocity(&mut full, Vector3::new(0.0, 5.0, 12.0)); + assert_eq!(speed::(&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!(>::NAME, "OrbitCdBias"); + assert_eq!(>::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!(>::NAME, "MassOnly"); + assert_eq!(>::MASS, 0); + + let mut y = Vector1::::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::::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); + } +}