diff --git a/Cargo.toml b/Cargo.toml index caa769a..3a8fd0b 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -11,7 +11,7 @@ keywords = ["astrodynamics", "orbit", "propagation", "spaceflight"] categories = ["science", "aerospace"] [features] -default = ["autodiff"] +default = ["autodiff", "msis"] autodiff = ["dep:num-dual"] # NRLMSISE-00 as a density model, from its own crate. Brings that crate's `fast` path, # which trades the reference's exact bits for roughly half the time: the largest value it @@ -23,6 +23,9 @@ msis = ["dep:nrlmsise00"] # We handle the kernel fetching, so we don't need anise's default features anise = { version = "0.10.6", default-features = false } hifitime = "4.3" +# A facade only. Nothing is logged unless the program installs a logger, and the one +# thing this crate reports that way is a density model falling back to a fixed profile. +log = "0.4" nalgebra = { version = "0.35", features = ["serde-serialize"] } num-dual = { version = "0.15", optional = true } nrlmsise00 = { version = "0.2", optional = true, features = ["drivers", "fast"] } @@ -33,6 +36,7 @@ thiserror = "2.0" [dev-dependencies] approx = "0.5" +log = { version = "0.4", features = ["std"] } serde_json = "1.0" criterion = "0.7" num-dual = "0.15" @@ -52,3 +56,13 @@ harness = false [[bench]] name = "reference_scenarios" harness = false + +# Measured 2026-10-03, because the unoptimized example took 95 seconds and the readme +# claimed six. The gravity evaluator is this crate's own code, so optimizing only +# dependencies buys little -- that was measured too, at 95 s to 62 s. At opt-level 1 the +# whole example is 6.3 s and `cargo test` goes from 2.8 s to 0.9 s, against 4 s added to +# a clean build (10.0 s to 14.3 s) and no measurable cost to an incremental one. Some +# inlining happens at this level, so a debugger steps less faithfully than at 0; that is +# the only thing given up. +[profile.dev] +opt-level = 1 diff --git a/examples/full.rs b/examples/full.rs index 8724c4d..d9bdb73 100644 --- a/examples/full.rs +++ b/examples/full.rs @@ -4,11 +4,8 @@ //! //! Run with `cargo run --example full`. Needs kernels — `scripts/fetch-kernels.sh`. //! -//! Add `--release` if you run it often. A day of degree-70 gravity takes about six -//! seconds unoptimized and two tenths of a second optimized. - -use std::path::PathBuf; -use std::sync::Arc; +//! About 6 seconds unoptimized and 3 optimized, measured 2026-10-03. It propagates the +//! day seven times use astrokit::dynamics::Dynamics; use astrokit::eclipse; @@ -16,24 +13,9 @@ use astrokit::prelude::*; use astrokit::{attitude, kernels}; use nalgebra::Vector3; -const DEGREE: usize = 70; - fn main() -> Result<(), Box> { - // `kernels::directory` reads ASTROKIT_KERNELS and nothing else, since a path - // relative to a source tree means nothing to a library. A program knows where it - // was run from and can fall back. - let directory = kernels::directory().unwrap_or_else(|| PathBuf::from("data/kernels")); - if !kernels::manifest() - .iter() - .all(|kernel| directory.join(kernel.name).is_file()) - { - println!( - "kernels not found in {}; run scripts/fetch-kernels.sh", - directory.display() - ); - return Ok(()); - } - let almanac = Arc::new(kernels::load(&directory)?); + // ASTROKIT_KERNELS if it is set, else data/kernels. + let almanac = kernels::load_default()?; let epoch = Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0); let frame = Frame::EARTH_J2000; @@ -60,23 +42,9 @@ fn main() -> Result<(), Box> { .with_mass(spacecraft.mass) .with_angular_rate(tumble); - // What a run records of itself. Plain data, so it serializes as the explicit force - // list rather than as a mode flag, and names the coefficient set rather than only - // how far it was truncated. - let forces = vec![ - Gravity { - field: GravityModel::Egm96, - degree: DEGREE, - order: DEGREE, - }, - ThirdBody { body: Body::MOON }, - ThirdBody { body: Body::SUN }, - Sunshine, - Drag { - density: DensityModel::harris_priester(), - }, - ]; - let field = GravityModel::Egm96.field(DEGREE, DEGREE)?; + // A recipe is a starting position, not a mode: it comes back as the plain force list + // this run then records of itself. + let forces = Regime::LowEarth.forces(Fidelity::Precise, epoch); let mut propagation = Propagation::new(epoch, initial, &forces) .frame(frame) @@ -84,6 +52,11 @@ fn main() -> Result<(), Box> { .almanac(almanac.clone()) .duration(arc) .accuracy(Accuracy::Medium); + + // Everything the run was given except the state. + println!("{}", propagation.provenance()); + println!(); + let solution = propagation.run(); let final_state = Full.state(*solution.states.last().unwrap()); @@ -108,6 +81,7 @@ fn main() -> Result<(), Box> { // exception — removing it leaves nothing to orbit — so for that one the harmonics are // swapped for a point mass instead. let labels = [ + // "gravity harmonics", "the moon", "the sun", @@ -146,27 +120,6 @@ fn main() -> Result<(), Box> { let inertia = &spacecraft.inertia; let momentum = |rate: Vector3| (inertia * rate).norm(); - println!("kernels {}", directory.display()); - println!( - "spacecraft {:.0} kg, {} faces, inertia diag({:.1}, {:.1}, {:.1}) kg m^2", - spacecraft.mass, - match &spacecraft.surface { - SurfaceModel::Faceted { faces } => faces.len(), - _ => 1, - }, - inertia[(0, 0)], - inertia[(1, 1)], - inertia[(2, 2)] - ); - println!( - "field EGM96 degree {} order {}, mu {:.9e}, reference radius {} m", - field.degree(), - field.order(), - field.mu(), - field.reference_radius() - ); - println!("epoch {epoch}"); - println!("arc {arc}"); println!( "initial a = {:.3} km, e = {:.6}, i = {:.3} deg, period {}", orbit.semi_major_axis / 1e3, @@ -277,12 +230,5 @@ fn main() -> Result<(), Box> { midway_ground.height / 1e3 ); - // The force list is plain data, so the record of what ran is the composition itself - // rather than the name of a mode. Each body carries the constants it was used with, - // so reading this back does not depend on those constants not having changed. - println!(); - println!("provenance"); - println!("{}", serde_json::to_string_pretty(&forces)?); - Ok(()) } diff --git a/readme.md b/readme.md index 117beda..2345984 100644 --- a/readme.md +++ b/readme.md @@ -201,8 +201,10 @@ cargo run --example attitude # orbit and attitude together, on a tumbling box cargo run --example full # a day of EGM96 gravity, with the ground track ``` -`full` needs kernels. It is worth adding `--release` if you run it often: a day of -degree-70 gravity takes six seconds unoptimized and two tenths of a second optimized. +`full` needs kernels. It takes about 6 seconds unoptimized and 3 with `--release`, +because it propagates the day seven times: once for the trajectory, once for a bare +point mass, and once per force with that force left out, which is how it prices each +of them. ## Development diff --git a/src/atmosphere/mod.rs b/src/atmosphere/mod.rs index 71a424b..5ceaa28 100644 --- a/src/atmosphere/mod.rs +++ b/src/atmosphere/mod.rs @@ -142,3 +142,15 @@ impl DensityModel { } } } + +impl std::fmt::Display for DensityModel { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self { + Self::HarrisPriester { bulge_exponent } => { + write!(f, "Harris-Priester (n = {bulge_exponent})") + } + #[cfg(feature = "msis")] + Self::Nrlmsise00 => write!(f, "NRLMSISE-00"), + } + } +} diff --git a/src/atmosphere/msis.rs b/src/atmosphere/msis.rs index 030b5b5..2a146f4 100644 --- a/src/atmosphere/msis.rs +++ b/src/atmosphere/msis.rs @@ -140,7 +140,7 @@ mod tests { /// An arc reaching past the measured record is refused before integrating, rather /// than extrapolated. No kernels needed: this check runs before the almanac one. #[test] - #[should_panic(expected = "reads space weather")] + #[should_panic(expected = "outside the measured space weather record")] fn test_an_arc_past_the_record_is_refused() { use crate::prelude::*; use crate::spacecraft::Spacecraft; diff --git a/src/dynamics/mod.rs b/src/dynamics/mod.rs index acd015b..14639f9 100644 --- a/src/dynamics/mod.rs +++ b/src/dynamics/mod.rs @@ -436,6 +436,34 @@ where specs.iter().map(Dynamics::build).collect() } +impl std::fmt::Display for Dynamics { + /// One line naming the model and what distinguishes it, for a record someone reads. + /// + /// Deliberately not the serialized form: [`Body`] carries four constants, so a + /// `ThirdBody` renders as a wall of numbers under `Debug` or serde and as the body's + /// name here. + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self { + Self::TwoBody { body } => write!(f, "point mass, {}", body.name()), + Self::Gravity { + field, + degree, + order, + } => { + // `non_exhaustive` binds downstream crates, not this one, so this match + // stays exhaustive and a new field has to be named here. + let name = match field { + GravityModel::Egm96 => "EGM96", + }; + write!(f, "gravity, {name} {degree}x{order}") + } + Self::ThirdBody { body } => write!(f, "third body, {}", body.name()), + Self::Drag { density } => write!(f, "drag, {density}"), + Self::Sunshine => write!(f, "sunlight (pressure, power, eclipse)"), + } + } +} + #[cfg(test)] mod tests { use super::*; diff --git a/src/kernels.rs b/src/kernels.rs index ffbf14a..a29db8b 100644 --- a/src/kernels.rs +++ b/src/kernels.rs @@ -18,6 +18,7 @@ //! `ASTROKIT_KERNELS` for convenience, but nothing here guesses. use std::path::{Path, PathBuf}; +use std::sync::Arc; use anise::almanac::Almanac; use thiserror::Error; @@ -114,6 +115,49 @@ pub(crate) fn test_directory() -> Option { .then_some(candidate) } +/// Loads the kernels from wherever they are, for a program that does not want to care. +/// +/// Looks in `ASTROKIT_KERNELS` if it is set, then in `data/kernels` relative to the +/// working directory. Returns an [`Arc`] because that is what +/// [`Propagation::almanac`](crate::propagation::Propagation::almanac) takes. +/// +/// [`directory`] deliberately does not fall back, because a relative path means nothing +/// to a library deciding on its own. This is the other half of that: a *caller* may +/// reasonably want the convention, and here it is in one call rather than spelled out +/// at every call site. +/// +/// # Errors +/// +/// If no candidate directory holds the whole manifest. The error names the last place +/// it looked and how to fetch them. +pub fn load_default() -> Result, KernelError> { + let candidates: Vec = directory() + .into_iter() + .chain(std::iter::once(PathBuf::from("data/kernels"))) + .collect(); + + let complete = |candidate: &PathBuf| { + manifest() + .iter() + .all(|kernel| candidate.join(kernel.name).is_file()) + }; + let chosen = candidates + .iter() + .find(|candidate| complete(candidate)) + .or_else(|| candidates.last()); + + match chosen { + Some(directory) => load(directory).map(Arc::new), + None => Err(KernelError::Missing { + name: manifest() + .first() + .map(|kernel| kernel.name.to_string()) + .unwrap_or_default(), + directory: PathBuf::from("data/kernels"), + }), + } +} + /// Loads every kernel in the manifest from `directory`. /// /// All or nothing: a missing kernel is an error rather than a quietly smaller diff --git a/src/lib.rs b/src/lib.rs index 09868e4..0d62185 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -35,6 +35,8 @@ pub mod kepler; pub mod kernels; pub mod measurement; pub mod propagation; +pub mod provenance; +pub mod recipes; pub mod spacecraft; pub mod state; pub mod station; @@ -54,8 +56,10 @@ pub mod prelude { pub use crate::geodetic::Geodetic; pub use crate::gravity::{Field, GravityModel}; pub use crate::kepler::{Degeneracy, Keplerian, TranslationalLayout}; - pub use crate::measurement::{Measurement, Observed, Prediction}; + pub use crate::measurement::{Corrections, Measurement, Observed, Prediction}; pub use crate::propagation::{Accuracy, Propagation}; + pub use crate::provenance::Provenance; + pub use crate::recipes::{Fidelity, Regime}; pub use crate::spacecraft::{Battery, Face, Mounting, SolarArray, Spacecraft, SurfaceModel}; pub use crate::state::{ Full, HasAngularRate, HasAttitude, HasEnergy, HasMass, HasPosition, HasVelocity, Layout, diff --git a/src/measurement.rs b/src/measurement.rs index 70acf60..723e4af 100644 --- a/src/measurement.rs +++ b/src/measurement.rs @@ -132,6 +132,18 @@ impl Prediction { } } +/// Which corrections an observation carries. +/// +/// Recorded rather than implied, because the corrections are a choice made per +/// observation and a result is not interpretable without knowing which were applied. +#[derive(Clone, Copy, Debug, Default, PartialEq, Eq, Serialize, Deserialize)] +pub struct Corrections { + /// Moved to where the target was when the light left it. + pub light_time: bool, + /// The line of sight tilted by the observer's own motion. + pub aberration: bool, +} + /// The geometry a measurement is read off, worked out once. /// /// Preparing this costs two frame rotations, so a set of measurements taken at one epoch @@ -145,6 +157,8 @@ pub struct Observed { pub relative: Cartesian, /// The station's own inertial state, which is what aberration is relative to. pub station: Cartesian, + /// Which corrections have been applied to it. + pub corrections: Corrections, /// Inertial to the station's east, north and up. horizon: Matrix3, } @@ -171,6 +185,7 @@ impl Observed { Ok(Self { relative, station: station.state_in::(inertial, epoch, almanac)?, + corrections: Corrections::default(), horizon: station.east_north_up::() * rotation.map(T::lift), }) } @@ -204,6 +219,7 @@ impl Observed { delay = (geometric - target_velocity * delay).norm() / LIGHT_SPEED; } self.relative.position = geometric - target_velocity * delay; + self.corrections.light_time = true; self } @@ -219,6 +235,7 @@ impl Observed { let line = self.relative.position / range; let tilted = line + self.station.velocity / T::lift(LIGHT_SPEED); self.relative.position = tilted.normalize() * range; + self.corrections.aberration = true; self } } @@ -548,6 +565,31 @@ mod tests { assert!((-FRAC_PI_2..=FRAC_PI_2).contains(&elevation), "{elevation}"); } + /// An observation says which corrections it carries, so a recorded measurement is + /// interpretable without knowing how the program that made it was built. + #[test] + fn test_an_observation_records_its_corrections() { + let Some(at) = observed() else { + eprintln!("skipped: kernels not fetched, run scripts/fetch-kernels.sh"); + return; + }; + assert_eq!(at.corrections, Corrections::default()); + assert!(!at.corrections.light_time && !at.corrections.aberration); + + let corrected = at.with_light_time(); + assert!(corrected.corrections.light_time); + assert!(!corrected.corrections.aberration); + + let both = corrected.with_aberration(); + assert!(both.corrections.light_time && both.corrections.aberration); + + let text = serde_json::to_string(&both.corrections).unwrap(); + assert_eq!( + serde_json::from_str::(&text).unwrap(), + both.corrections + ); + } + /// A measurement is part of a run's record, so it serializes. #[test] fn test_a_measurement_round_trips_through_json() { diff --git a/src/propagation.rs b/src/propagation.rs index e6d88d4..08be319 100644 --- a/src/propagation.rs +++ b/src/propagation.rs @@ -13,6 +13,7 @@ use crate::cartesian::Cartesian; use crate::dynamics::{self, Context, Contribution, Dynamics, DynamicsModel}; use crate::events; use crate::frames::Frame; +use crate::provenance::{self, Provenance}; use crate::spacecraft::{Spacecraft, SurfaceModel}; use crate::state::{HasPosition, HasVelocity, Layout, Orbit, State}; @@ -73,6 +74,9 @@ impl Accuracy { } } +/// The tolerance a propagation starts at, before any [`Accuracy`] preset is applied. +const DEFAULT_TOLERANCE: f64 = 1e-10; + /// The magnitude in a component's own units that its absolute tolerance is scaled by. /// /// One for everything whose natural size is around one or larger and which does not pass @@ -245,6 +249,9 @@ pub struct Propagation<'a, T: Scalar + 'static, const D: usize, L: Layout + ' controller: PIController, almanac: Option>, shadow_set_switching: bool, + /// Mirrors what was pushed into `method`, which has no way to report it back. + relative_tolerance: f64, + absolute_tolerances: SVector, callbacks: Vec<&'a mut dyn Callback>>, problem: Option>, } @@ -266,10 +273,15 @@ where spacecraft: None, duration: None, specs: dynamics.to_vec(), - method: MethodName::Vern7.instantiate().a_tol(1e-10).r_tol(1e-10), + method: MethodName::Vern7 + .instantiate() + .a_tol(DEFAULT_TOLERANCE) + .r_tol(DEFAULT_TOLERANCE), controller: PIController::default(), almanac: None, shadow_set_switching: true, + relative_tolerance: DEFAULT_TOLERANCE, + absolute_tolerances: SVector::repeat(DEFAULT_TOLERANCE), callbacks: Vec::new(), problem: None, } @@ -373,28 +385,32 @@ where /// Apply it before [`Propagation::a_tol`], [`Propagation::a_tol_per_component`] or /// [`Propagation::r_tol`] if you want to override one of them. pub fn accuracy(mut self, accuracy: Accuracy) -> Self { - let tol = accuracy.tolerance(); + self.relative_tolerance = accuracy.tolerance(); + self.absolute_tolerances = accuracy.absolute_tolerances::(); self.method = self .method - .a_tol_full(accuracy.absolute_tolerances::()) - .r_tol(tol); + .a_tol_full(self.absolute_tolerances) + .r_tol(self.relative_tolerance); self } /// Sets the absolute tolerance for every component. pub fn a_tol(mut self, a_tol: f64) -> Self { + self.absolute_tolerances = SVector::repeat(a_tol); self.method = self.method.a_tol(a_tol); self } /// Sets the absolute tolerance per component. pub fn a_tol_per_component(mut self, a_tol: SVector) -> Self { + self.absolute_tolerances = a_tol; self.method = self.method.a_tol_full(a_tol); self } /// Sets the relative tolerance. pub fn r_tol(mut self, r_tol: f64) -> Self { + self.relative_tolerance = r_tol; self.method = self.method.r_tol(r_tol); self } @@ -405,11 +421,16 @@ where self } - /// Replaces the integrator, at its default tolerances. + /// Replaces the integrator, keeping the tolerances. /// - /// Tolerances travel with the integrator, so set them after this, not before. + /// Tolerances belong to the integrator in `ordinary-diffeq`, so instantiating a new + /// one would otherwise silently revert them to that integrator's own defaults and + /// make the order of the builder calls matter. pub fn method(mut self, method: MethodName) -> Self { - self.method = method.instantiate(); + self.method = method + .instantiate() + .a_tol_full(self.absolute_tolerances) + .r_tol(self.relative_tolerance); self.problem = None; self } @@ -424,6 +445,29 @@ where &self.specs } + /// What this propagation was set up with, as a record that serializes. + /// + /// Can be taken before or after running; nothing about the setup changes during a + /// run. See [`Provenance`] for what it does and does not promise. + pub fn provenance(&self) -> Provenance { + Provenance { + astrokit: provenance::VERSION.to_string(), + epoch: self.epoch, + duration: self.duration, + frame: self.frame, + layout: L::NAME.to_string(), + dimension: D, + forces: self.specs.clone(), + method: self.method.name(), + relative_tolerance: self.relative_tolerance, + absolute_tolerances: self.absolute_tolerances.iter().copied().collect(), + spacecraft: self.spacecraft.clone(), + kernels_attached: self.almanac.is_some(), + kernels: Provenance::manifest(), + shadow_set_switching: self.shadow_set_switching, + } + } + /// Runs the propagation. /// /// The returned solution borrows this `Propagation`, because @@ -479,7 +523,19 @@ where let covered = nrlmsise00::drivers::Drivers::embedded() .covers_arc(self.epoch.to_mjd_utc_days(), end.to_mjd_utc_days()); if let Err(error) = covered { - panic!("{spec:?} reads space weather: {error}"); + panic!( + "this arc is outside the measured space weather record, so \ + NRLMSISE-00 cannot be evaluated over it: {error}\n\ + The usable record stops 41 days before its last measured \ + day, because the 81-day centered flux average needs that \ + many days after the epoch. An arc at or near the present is \ + therefore already past it, which is expected rather than a \ + fault: extrapolating space weather would produce a density \ + nobody measured.\n\ + For an arc outside the record use \ + `DensityModel::harris_priester()`, a fixed profile that \ + reads no weather and so has no end date." + ); } } } diff --git a/src/provenance.rs b/src/provenance.rs new file mode 100644 index 0000000..b85a48a --- /dev/null +++ b/src/provenance.rs @@ -0,0 +1,394 @@ +//! What produced a run. +//! +//! A record of everything a propagation was given except the state itself: the resolved +//! force list, the integrator and its tolerances, the vehicle, the kernels and the +//! version of this crate. It serializes, so it can sit beside a result and be read back +//! later or compared against another run's. +//! +//! # What it promises, and what it does not +//! +//! It records **what was used**, not that a rerun will reproduce the same bits. Floating +//! point propagation does not offer that, and neither does the serialization: a `f64` +//! written to JSON and read back can differ in the last place, because not every parser +//! is correctly rounded. What the record is for is the question that actually comes up — +//! *were these two runs the same setup?* — and for catching the answer "no" when two +//! results disagree. +//! +//! It also describes the propagation and not the measurements taken from it. Which +//! corrections an observation carried is recorded on the observation, in +//! [`Observed::corrections`](crate::measurement::Observed::corrections). + +use hifitime::{Duration, Epoch}; +use ordinary_diffeq::prelude::MethodName; +use serde::{Deserialize, Serialize}; + +use crate::dynamics::Dynamics; +use crate::frames::Frame; +use crate::kernels; +use crate::spacecraft::Spacecraft; + +/// The version of this crate, as the record writes it. +pub const VERSION: &str = env!("CARGO_PKG_VERSION"); + +/// A kernel by name and expected digest. +#[derive(Clone, Debug, PartialEq, Eq, Serialize, Deserialize)] +pub struct KernelRecord { + /// File name. + pub name: String, + /// Expected SHA-256, lowercase hex. + pub sha256: String, +} + +/// What a propagation was set up with. +/// +/// Produced by [`Propagation::provenance`](crate::propagation::Propagation::provenance), +/// which can be called before or after running. +#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)] +pub struct Provenance { + /// Version of this crate. The versions of its dependencies are not reachable at + /// runtime; this plus a lockfile is what identifies the build. + pub astrokit: String, + /// Epoch the arc started from. + pub epoch: Epoch, + /// How long it was asked to run, if a duration was set. + pub duration: Option, + /// Frame the state was expressed in. + pub frame: Frame, + /// Name of the state layout, and how many elements it has. + pub layout: String, + /// Dimension of the state vector. + pub dimension: usize, + /// The force list as given, in order. This is the resolved composition: a gravity + /// entry names its field and its truncation rather than implying them. + pub forces: Vec, + /// Which integrator ran. + pub method: MethodName, + /// Relative tolerance the step size controller was held to. + pub relative_tolerance: f64, + /// Absolute tolerance per state element, in each element's own units. + pub absolute_tolerances: Vec, + /// The vehicle, if one was attached. + pub spacecraft: Option, + /// Whether kernels were attached at all. + pub kernels_attached: bool, + /// What [`kernels::load`] loads, by name and digest. + /// + /// Describes the kernel set this crate ships a manifest for, which is what + /// `kernels_attached` will have been built from in the ordinary case. An almanac + /// assembled some other way is not described by it, which is why the flag and the + /// list are separate. + pub kernels: Vec, + /// Whether the default shadow-set switch was left on. + pub shadow_set_switching: bool, +} + +impl Provenance { + /// The kernel manifest, as the record holds it. + pub(crate) fn manifest() -> Vec { + kernels::manifest() + .into_iter() + .map(|kernel| KernelRecord { + name: kernel.name.to_string(), + sha256: kernel.sha256.to_string(), + }) + .collect() + } +} + +impl std::fmt::Display for Provenance { + /// The record laid out for a human, which is what printing one is almost always for. + /// + /// Not the serialized form. `Debug` and serde both render a [`Body`](crate::body::Body) + /// as its four constants, so a force list naming the Sun and Moon becomes a wall of + /// numbers; here each force is one line. Use serde when something else has to read + /// it, and this when a person does. + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + writeln!(f, "astrokit {}", self.astrokit)?; + writeln!(f, " epoch {}", self.epoch)?; + match self.duration { + Some(duration) => writeln!(f, " duration {duration}")?, + None => writeln!(f, " duration until an epoch given at run time")?, + } + writeln!( + f, + " frame {} {:?}", + self.frame.center.name(), + self.frame.orientation + )?; + writeln!(f, " state {} ({})", self.layout, self.dimension)?; + + writeln!(f, " integrator {:?}", self.method)?; + writeln!(f, " r_tol {:e}", self.relative_tolerance)?; + // Per component, and usually all the same, so say so rather than printing a row + // of identical numbers. Only energy's floor normally differs — see `Accuracy`. + let first = self.absolute_tolerances.first().copied().unwrap_or(0.0); + if self.absolute_tolerances.iter().all(|t| *t == first) { + writeln!(f, " a_tol {first:e}")?; + } else { + // Grouped into runs, because a wide layout is typically one floor repeated + // with a single component differing — energy, whose natural size is 1e8 J. + // Printing fourteen numbers to show two would bury the one that matters. + let mut runs: Vec<(f64, usize)> = Vec::new(); + for tolerance in &self.absolute_tolerances { + match runs.last_mut() { + Some((value, count)) if value == tolerance => *count += 1, + _ => runs.push((*tolerance, 1)), + } + } + let listed: Vec = runs + .iter() + .map(|(value, count)| match count { + 1 => format!("{value:e}"), + _ => format!("{value:e} x{count}"), + }) + .collect(); + writeln!(f, " a_tol {}", listed.join(", "))?; + } + + for (index, force) in self.forces.iter().enumerate() { + let label = if index == 0 { "forces" } else { "" }; + writeln!(f, " {label:<12} {force}")?; + } + if self.forces.is_empty() { + writeln!(f, " forces none")?; + } + + match &self.spacecraft { + Some(vehicle) => writeln!( + f, + " spacecraft {:.3} kg ({:.3} kg dry), {}", + vehicle.mass, vehicle.dry_mass, vehicle.surface + )?, + None => writeln!(f, " spacecraft none")?, + } + + if self.kernels_attached { + for (index, kernel) in self.kernels.iter().enumerate() { + let label = if index == 0 { "kernels" } else { "" }; + writeln!(f, " {:<12} {}", label, kernel.name)?; + } + if self.kernels.is_empty() { + writeln!(f, " kernels attached, manifest empty")?; + } + } else { + writeln!(f, " kernels none attached")?; + } + + write!( + f, + " shadow set {}", + if self.shadow_set_switching { + "switching on" + } else { + "switching off" + } + ) + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::prelude::*; + use crate::recipes::{self, Fidelity, Regime}; + use crate::spacecraft::Spacecraft; + use nalgebra::Vector3; + + /// Inside the measured space weather record, so a low-orbit recipe resolves to + /// NRLMSISE-00 rather than falling back to the fixed profile. + fn epoch() -> Epoch { + Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0) + } + + fn propagation<'a>(forces: &'a [Dynamics]) -> Propagation<'a, f64, 6, Orbit> { + Propagation::new( + epoch(), + Orbit + .with_position(Vector3::new(6.978e6, 0.0, 0.0)) + .with_velocity(Vector3::new(0.0, 7.56e3, 0.0)), + forces, + ) + .duration(Duration::from_seconds(3600.0)) + } + + /// The record is what sits beside a result, so it has to survive being written out + /// and read back. + #[test] + fn test_it_round_trips_through_json() { + let forces = Regime::LowEarth.forces(Fidelity::Precise, epoch()); + let record = propagation(&forces) + .spacecraft(Spacecraft::cannonball(500.0, 5.0, 2.2, 1.3)) + .accuracy(Accuracy::High) + .provenance(); + + let text = serde_json::to_string_pretty(&record).unwrap(); + assert_eq!(serde_json::from_str::(&text).unwrap(), record); + } + + /// What it actually carries: the resolved composition rather than a mode name. + #[test] + fn test_it_records_the_resolved_setup() { + let forces = Regime::LowEarth.forces(Fidelity::Precise, epoch()); + let record = propagation(&forces).accuracy(Accuracy::High).provenance(); + + assert_eq!(record.astrokit, VERSION); + assert_eq!(record.forces, forces); + assert_eq!(record.layout, "Orbit"); + assert_eq!(record.dimension, 6); + assert_eq!(record.relative_tolerance, Accuracy::High.tolerance()); + assert_eq!(record.absolute_tolerances.len(), 6); + assert!(!record.kernels_attached); + assert!(!record.kernels.is_empty(), "the manifest should be listed"); + + // The gravity entry names its field and its truncation, so the record says which + // coefficients ran rather than implying them. + let gravity = record + .forces + .iter() + .find(|force| matches!(force, Gravity { .. })) + .expect("the recipe has a field"); + assert!(format!("{gravity:?}").contains("Egm96")); + } + + /// Two runs set up the same way record the same thing, and one difference shows. + #[test] + fn test_it_distinguishes_two_setups() { + let low = Regime::LowEarth.forces(Fidelity::Precise, epoch()); + let geo = Regime::Geostationary.forces(Fidelity::Precise, epoch()); + assert_eq!( + propagation(&low).provenance(), + propagation(&low).provenance() + ); + assert_ne!( + propagation(&low).provenance(), + propagation(&geo).provenance() + ); + // Including a difference that is only in the integrator settings. + assert_ne!( + propagation(&low).accuracy(Accuracy::Low).provenance(), + propagation(&low).accuracy(Accuracy::High).provenance() + ); + } + + /// Choosing an integrator keeps the tolerances. They live on the integrator in + /// `ordinary-diffeq`, so instantiating a new one would otherwise revert them and make + /// the order of the builder calls matter. + #[test] + fn test_the_integrator_can_be_chosen_in_any_order() { + let forces = recipes::two_body(Body::EARTH); + let before = propagation(&forces) + .accuracy(Accuracy::High) + .method(MethodName::DormandPrince45) + .provenance(); + let after = propagation(&forces) + .method(MethodName::DormandPrince45) + .accuracy(Accuracy::High) + .provenance(); + + assert_eq!(before, after); + assert_eq!(before.method, MethodName::DormandPrince45); + assert_eq!(before.relative_tolerance, Accuracy::High.tolerance()); + } + + /// Attaching kernels shows in the record, since a run without them is a different run. + #[test] + fn test_kernels_show_as_attached() { + let Some(directory) = crate::kernels::test_directory() else { + eprintln!("skipped: kernels not fetched, run scripts/fetch-kernels.sh"); + return; + }; + let almanac = std::sync::Arc::new(crate::kernels::load(&directory).unwrap()); + let forces = recipes::two_body(Body::EARTH); + assert!( + propagation(&forces) + .almanac(almanac) + .provenance() + .kernels_attached + ); + assert!(!propagation(&forces).provenance().kernels_attached); + } +} + +#[cfg(test)] +mod display_tests { + use super::*; + use crate::body::Body; + use crate::dynamics::Dynamics::{Drag, Gravity, Sunshine, ThirdBody, TwoBody}; + use crate::gravity::GravityModel; + + /// Every force renders on one line, naming itself and what distinguishes it. + /// + /// The reason this exists rather than leaning on `Debug`: a `Body` carries four + /// constants, so the derived form prints a wall of numbers for what a reader wants + /// to see as "Moon". + #[test] + fn test_a_force_renders_as_one_readable_line() { + let moon = ThirdBody { body: Body::MOON }.to_string(); + assert_eq!(moon, "third body, Moon"); + assert!(!moon.contains("naif_id"), "the body's constants leaked in"); + assert!(!moon.contains('\n'), "a force should be one line"); + + assert_eq!( + Gravity { + field: GravityModel::Egm96, + degree: 70, + order: 70, + } + .to_string(), + "gravity, EGM96 70x70" + ); + assert_eq!( + TwoBody { body: Body::EARTH }.to_string(), + "point mass, Earth" + ); + assert!(Sunshine.to_string().starts_with("sunlight")); + assert!(Drag { + density: crate::atmosphere::DensityModel::harris_priester(), + } + .to_string() + .starts_with("drag, Harris-Priester")); + } + + /// The record prints as a block a person can read, with every force named. + #[test] + fn test_the_record_renders_every_force_and_leaks_no_constants() { + let forces = crate::recipes::Regime::LowEarth.forces( + crate::recipes::Fidelity::Precise, + Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0), + ); + let record = Provenance { + astrokit: VERSION.to_string(), + epoch: Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0), + duration: Some(Duration::from_seconds(86400.0)), + frame: Frame::EARTH_J2000, + layout: "Orbit".to_string(), + dimension: 6, + forces: forces.clone(), + method: MethodName::Vern7, + relative_tolerance: 1e-12, + absolute_tolerances: vec![1e-12; 6], + spacecraft: None, + kernels_attached: false, + kernels: Vec::new(), + shadow_set_switching: true, + }; + + let shown = record.to_string(); + for force in &forces { + assert!( + shown.contains(&force.to_string()), + "{force} is missing from the record" + ); + } + assert!(shown.contains(VERSION)); + assert!(shown.contains("Orbit")); + // The thing `Debug` gets wrong, which is the whole reason for this impl. + assert!( + !shown.contains("naif_id") && !shown.contains("equatorial_radius"), + "a body's constants leaked into the human-readable form" + ); + // Identical per-component tolerances collapse to one line rather than six. + assert_eq!(shown.matches("a_tol").count(), 1); + assert!(!shown.contains("1e-12, 1e-12")); + } +} diff --git a/src/recipes.rs b/src/recipes.rs new file mode 100644 index 0000000..68c8ba3 --- /dev/null +++ b/src/recipes.rs @@ -0,0 +1,405 @@ +//! Force lists for common orbit regimes, at a chosen level of detail. +//! +//! A recipe returns a plain `Vec`, not a mode or a handle. That is the point: +//! it is a starting position someone can read, print, edit and record, not a setting +//! whose meaning is hidden inside the library. Swapping the density model or dropping +//! the Moon is a `Vec` operation. +//! +//! ``` +//! use astrokit::prelude::*; +//! use astrokit::recipes::{Fidelity, Regime}; +//! +//! let epoch = Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0); +//! let mut forces = Regime::LowEarth.forces(Fidelity::Balanced, epoch); +//! // The list is ordinary data, so adjusting it needs nothing from this crate. +//! forces.retain(|force| !matches!(force, Drag { .. })); +//! ``` +//! +//! # Why there are two knobs and not one +//! +//! [`Regime`] says where the orbit is and [`Fidelity`] says how much detail to spend on +//! it. They are not interchangeable and neither alone is enough: at a few hundred +//! kilometers the field and the atmosphere dominate, while at geostationary there is no +//! atmosphere worth the name and the Sun and Moon drive almost everything. A single +//! "high fidelity" dial that added forces would be wrong in a way that depends on the +//! regime, which is why it does not exist here. +//! +//! What [`Fidelity`] varies is how much of each force is kept, measured against what +//! that force is actually worth. Over a day at 400 km the hierarchy is the field at +//! about 1030 km of displacement, drag at 35 km, the Sun at 35 m, the Moon at 6 m and +//! radiation pressure at 6 m. So [`Fidelity::Fast`] keeps the first two and drops terms +//! three orders of magnitude below them, while higher levels keep everything and spend +//! the difference on the gravity truncation. +//! +//! This stays explicit the whole way: a recipe returns the resolved list, so the record +//! a run writes names every force and every degree. Nothing about the composition is +//! recoverable only by knowing which level was asked for. +//! +//! # On the gravity degrees +//! +//! A spherical harmonic term of degree `n` falls off as `(R / r)^n`, so what is worth +//! keeping depends entirely on altitude. At a few hundred kilometers `R / r` is about +//! 0.94 and degree 70 still contributes; at geostationary it is 0.15, and degree 8 is +//! already below anything else in the list. Carrying 70 there would cost about twenty +//! times the gravity evaluation for nothing. + +use hifitime::Epoch; + +use crate::atmosphere::DensityModel; +use crate::body::Body; +use crate::dynamics::Dynamics::{self, Drag, Gravity, Sunshine, ThirdBody, TwoBody}; +use crate::gravity::GravityModel; + +/// Where the orbit is, which decides which forces are worth composing at all. +#[non_exhaustive] +#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)] +pub enum Regime { + /// A few hundred kilometers up: a close, rough field and a real atmosphere. + LowEarth, + /// Around half an Earth radius up, where a navigation constellation sits. No + /// atmosphere worth modeling and a noticeably smoother field. + MediumEarth, + /// Circular and equatorial at about 35,786 km, where the Sun and Moon drive the + /// inclination that station keeping is mostly spent against. + Geostationary, +} + +/// How much detail to spend, within a [`Regime`]. +/// +/// This is about the model, not about the integrator. How precisely the equations are +/// solved is [`Accuracy`](crate::propagation::Accuracy), and the two are deliberately +/// separate: converging a solver against a fixed set of forces is exactly how a +/// comparison against another tool is made, and a single combined dial would make that +/// unaskable. +#[non_exhaustive] +#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash, PartialOrd, Ord)] +pub enum Fidelity { + /// The dominant terms only, at a coarse truncation. For a first look, or for an + /// ensemble large enough that the per-run cost is what binds. + Fast, + /// Every force that matters in the regime, at a truncation that is well inside what + /// the other terms contribute. The sensible default. + Balanced, + /// Every force, at a truncation chosen so the field is not the limiting error. + Precise, +} + +impl Regime { + /// The resolved force list for this regime at this level of detail. + /// + /// The epoch decides the density model, and only that — see [`default_density`]. It + /// is the start of the arc; a long arc can still run past the measured record, and + /// [`Propagation::run`](crate::propagation::Propagation::run) says so if it does. + /// + /// Drag does not appear in MEO or GEO because the effect is literally below floating point + /// precision + pub fn forces(self, fidelity: Fidelity, epoch: Epoch) -> Vec { + let mut forces = vec![Gravity { + field: GravityModel::Egm96, + degree: self.degree(fidelity), + order: self.degree(fidelity), + }]; + + // At low altitude the third bodies and radiation pressure are three orders of + // magnitude under drag, so `Fast` spends nothing on them. Higher up they are the + // whole perturbation budget and are kept at every level. + let minor_terms_matter = self != Regime::LowEarth || fidelity != Fidelity::Fast; + if minor_terms_matter { + forces.push(ThirdBody { body: Body::MOON }); + forces.push(ThirdBody { body: Body::SUN }); + forces.push(Sunshine); + } + + if self == Regime::LowEarth { + forces.push(Drag { + density: default_density(epoch), + }); + } + forces + } + + /// The gravity truncation this regime uses at this level of detail. + pub fn degree(self, fidelity: Fidelity) -> usize { + match (self, fidelity) { + (Regime::LowEarth, Fidelity::Fast) => 2, + (Regime::LowEarth, Fidelity::Balanced) => 20, + (Regime::LowEarth, Fidelity::Precise) => 70, + (Regime::MediumEarth, Fidelity::Fast) => 2, + (Regime::MediumEarth, Fidelity::Balanced) => 4, + (Regime::MediumEarth, Fidelity::Precise) => 16, + (Regime::Geostationary, Fidelity::Fast) => 2, + (Regime::Geostationary, Fidelity::Balanced) => 4, + (Regime::Geostationary, Fidelity::Precise) => 8, + } + } +} + +/// The density model a low orbit gets at `epoch`. +/// +/// NRLMSISE-00 when this crate is built with the `msis` feature, which it is by +/// default, because a drag model that ignores solar and geomagnetic conditions is the +/// weakest part of any low-orbit propagation and there is no reason to hand someone the +/// worse one. +/// +/// # Why this depends on the epoch +/// +/// NRLMSISE-00 reads a measured record and will not extrapolate past it. The usable end +/// of that record is 41 days before its last measured day, because the 81-day centered +/// flux average needs that many days after the epoch — so an arc at or near the present +/// is *always* outside it, however fresh the data is. A recipe that returned +/// NRLMSISE-00 regardless would hand almost every forward-looking run a force list that +/// cannot be propagated. +/// +/// So outside the record this falls back to Harris-Priester, which is a fixed profile +/// with no end date, and **logs a warning** saying so. The fallback is a real change of +/// physics and is never silent: it is reported through [`log`] at the time, and the +/// resolved list it produces names Harris-Priester, so a run's provenance records which +/// model actually ran. +/// +/// Choosing explicitly rather than through a recipe avoids the question entirely. +pub fn default_density(epoch: Epoch) -> DensityModel { + #[cfg(feature = "msis")] + { + if nrlmsise00::drivers::Drivers::embedded().covers(epoch.to_mjd_utc_days()) { + return DensityModel::Nrlmsise00; + } + log::warn!( + "no measured space weather covers {epoch}, so this recipe is falling back \ + from NRLMSISE-00 to Harris-Priester, a fixed profile that ignores solar and \ + geomagnetic conditions. The usable record ends 41 days before its last \ + measured day, so an arc near the present is always outside it. Choose the \ + density model explicitly to avoid this." + ); + } + let _ = epoch; + DensityModel::harris_priester() +} + +/// A point mass and nothing else. +/// +/// The closed-form case, for a test that wants an orbit that closes. Not a regime: it +/// describes the model rather than where the orbit is. +pub fn two_body(body: Body) -> Vec { + vec![TwoBody { body }] +} + +#[cfg(test)] +mod tests { + use super::*; + + fn epoch() -> Epoch { + // Inside the measured space weather record, so the low-orbit recipe resolves to + // NRLMSISE-00 rather than falling back. + Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0) + } + + const EVERY: [Regime; 3] = [Regime::LowEarth, Regime::MediumEarth, Regime::Geostationary]; + const LEVELS: [Fidelity; 3] = [Fidelity::Fast, Fidelity::Balanced, Fidelity::Precise]; + + /// One gravity entry and no point mass alongside it. The monopole is the degree-zero + /// term of the same expansion, so composing both would count it twice. + #[test] + fn test_each_recipe_names_gravity_once() { + for regime in EVERY { + for fidelity in LEVELS { + let forces = regime.forces(fidelity, epoch()); + let fields = forces + .iter() + .filter(|force| matches!(force, Gravity { .. })) + .count(); + let monopoles = forces + .iter() + .filter(|force| matches!(force, TwoBody { .. })) + .count(); + assert_eq!(fields, 1, "{regime:?} {fidelity:?}"); + assert_eq!( + monopoles, 0, + "{regime:?} {fidelity:?}: a point mass beside a field" + ); + } + } + } + + /// The truncation never rises with altitude, because a term of degree n falls off as + /// (R/r)^n and there is less of the field left to resolve further out. + /// + /// Non-increasing rather than strictly falling: the coarse end of the dial bottoms + /// out at degree 2 in every regime, since J2 alone is the floor below which there is + /// no point keeping a field at all, and a strict comparison would forbid that. + /// Strictly falling is asserted at [`Fidelity::Precise`], where there is room for it. + #[test] + fn test_the_degree_never_rises_with_altitude() { + for fidelity in LEVELS { + assert!( + Regime::LowEarth.degree(fidelity) >= Regime::MediumEarth.degree(fidelity), + "{fidelity:?}" + ); + assert!( + Regime::MediumEarth.degree(fidelity) >= Regime::Geostationary.degree(fidelity), + "{fidelity:?}" + ); + } + + let precise = Fidelity::Precise; + assert!(Regime::LowEarth.degree(precise) > Regime::MediumEarth.degree(precise)); + assert!(Regime::MediumEarth.degree(precise) > Regime::Geostationary.degree(precise)); + } + + /// And it rises with the level asked for, in every regime. + #[test] + fn test_the_degree_rises_with_fidelity() { + for regime in EVERY { + assert!(regime.degree(Fidelity::Fast) < regime.degree(Fidelity::Balanced)); + assert!(regime.degree(Fidelity::Balanced) < regime.degree(Fidelity::Precise)); + } + } + + /// Drag belongs to the one regime with an atmosphere worth the name, at every level. + #[test] + fn test_only_the_low_orbit_has_drag() { + let drags = |forces: Vec| { + forces + .iter() + .filter(|force| matches!(force, Drag { .. })) + .count() + }; + for fidelity in LEVELS { + assert_eq!( + drags(Regime::LowEarth.forces(fidelity, epoch())), + 1, + "{fidelity:?}" + ); + assert_eq!( + drags(Regime::MediumEarth.forces(fidelity, epoch())), + 0, + "{fidelity:?}" + ); + assert_eq!( + drags(Regime::Geostationary.forces(fidelity, epoch())), + 0, + "{fidelity:?}" + ); + } + } + + /// Raising the level never removes a force. It is a dial on detail, not a different + /// model, and a user comparing two levels should be comparing like with more. + #[test] + fn test_fidelity_only_ever_adds() { + for regime in EVERY { + let fast = regime.forces(Fidelity::Fast, epoch()).len(); + let balanced = regime.forces(Fidelity::Balanced, epoch()).len(); + let precise = regime.forces(Fidelity::Precise, epoch()).len(); + assert!(fast <= balanced, "{regime:?}"); + assert!(balanced <= precise, "{regime:?}"); + } + } + + /// The cheap low-orbit list drops what is three orders of magnitude under drag, and + /// keeps what is not. + #[test] + fn test_a_fast_low_orbit_keeps_the_field_and_the_atmosphere() { + let forces = Regime::LowEarth.forces(Fidelity::Fast, epoch()); + assert!(forces.iter().any(|f| matches!(f, Gravity { .. }))); + assert!(forces.iter().any(|f| matches!(f, Drag { .. }))); + assert!(!forces.iter().any(|f| matches!(f, ThirdBody { .. }))); + assert!(!forces.iter().any(|f| matches!(f, Sunshine))); + } + + /// Further out the third bodies are the budget, so they survive even the cheap list. + #[test] + fn test_a_fast_distant_orbit_keeps_the_third_bodies() { + for regime in [Regime::MediumEarth, Regime::Geostationary] { + let forces = regime.forces(Fidelity::Fast, epoch()); + assert_eq!( + forces + .iter() + .filter(|f| matches!(f, ThirdBody { .. })) + .count(), + 2, + "{regime:?}" + ); + } + } + + /// A recipe is data, not a mode: it can be read, edited and recorded. + #[test] + fn test_a_recipe_is_editable_data() { + let mut forces = Regime::LowEarth.forces(Fidelity::Precise, epoch()); + let before = forces.len(); + forces.retain(|force| !matches!(force, ThirdBody { .. })); + assert_eq!(forces.len(), before - 2); + + // And it survives a round trip, which is what makes it recordable. + let recipe = Regime::LowEarth.forces(Fidelity::Precise, epoch()); + let text = serde_json::to_string(&recipe).unwrap(); + assert_eq!( + serde_json::from_str::>(&text).unwrap(), + recipe + ); + } + + /// The measured atmosphere is what a low orbit gets when it is available, because a + /// density model that ignores space weather is the weakest part of the propagation. + #[test] + fn test_the_low_orbit_uses_the_measured_atmosphere_when_it_is_built_in() { + let forces = Regime::LowEarth.forces(Fidelity::Balanced, epoch()); + let density = forces.iter().find_map(|force| match force { + Drag { density } => Some(*density), + _ => None, + }); + #[cfg(feature = "msis")] + assert_eq!(density, Some(DensityModel::Nrlmsise00)); + #[cfg(not(feature = "msis"))] + assert_eq!(density, Some(DensityModel::harris_priester())); + } + + /// Outside the measured record the low-orbit recipe falls back rather than handing + /// back a force list that cannot be propagated. + /// + /// The usable record ends 41 days before its last measured day, so "now" is always + /// outside it. This test uses a date far enough ahead that it stays outside however + /// fresh the embedded drivers become. + #[cfg(feature = "msis")] + #[test] + fn test_an_epoch_past_the_record_falls_back_to_the_fixed_profile() { + let distant = Epoch::from_gregorian_utc_hms(2100, 1, 1, 0, 0, 0); + assert_eq!(default_density(distant), DensityModel::harris_priester()); + + let forces = Regime::LowEarth.forces(Fidelity::Balanced, distant); + let density = forces.iter().find_map(|force| match force { + Drag { density } => Some(*density), + _ => None, + }); + assert_eq!(density, Some(DensityModel::harris_priester())); + } + + /// And inside it, the measured model is what comes back. + #[cfg(feature = "msis")] + #[test] + fn test_an_epoch_inside_the_record_gets_the_measured_atmosphere() { + assert_eq!(default_density(epoch()), DensityModel::Nrlmsise00); + } + + /// The fallback is a change of physics, so it has to be visible in what the recipe + /// returns rather than only in the log line. A record written from this list names + /// Harris-Priester, which is what actually ran. + #[cfg(feature = "msis")] + #[test] + fn test_the_fallback_is_recorded_in_the_list_itself() { + let distant = Epoch::from_gregorian_utc_hms(2100, 1, 1, 0, 0, 0); + let shown = Regime::LowEarth + .forces(Fidelity::Balanced, distant) + .iter() + .map(|force| force.to_string()) + .collect::>() + .join("; "); + assert!(shown.contains("Harris-Priester"), "{shown}"); + assert!(!shown.contains("NRLMSISE"), "{shown}"); + } + + #[test] + fn test_two_body_is_only_that() { + assert_eq!(two_body(Body::EARTH), vec![TwoBody { body: Body::EARTH }]); + } +} diff --git a/src/spacecraft.rs b/src/spacecraft.rs index 7e34540..664bc7d 100644 --- a/src/spacecraft.rs +++ b/src/spacecraft.rs @@ -391,6 +391,22 @@ impl Spacecraft { } } +impl std::fmt::Display for SurfaceModel { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self { + Self::Cannonball { + area, + drag_coefficient, + srp_coefficient, + } => write!( + f, + "cannonball (area {area} m^2, Cd {drag_coefficient}, Cr {srp_coefficient})" + ), + Self::Faceted { faces } => write!(f, "faceted ({} faces)", faces.len()), + } + } +} + #[cfg(test)] mod tests { use super::*; diff --git a/tests/density_fallback_warning.rs b/tests/density_fallback_warning.rs new file mode 100644 index 0000000..2974c83 --- /dev/null +++ b/tests/density_fallback_warning.rs @@ -0,0 +1,65 @@ +//! The recipe's density fallback is reported, not silent. +//! +//! [`recipes::default_density`] returns NRLMSISE-00 inside the measured space weather +//! record and Harris-Priester outside it. That is a real change of physics, so it is +//! logged at warning level. This lives in its own integration test because installing a +//! logger is process-wide and each test file is its own binary, so nothing else is +//! disturbed by it. +//! +//! What the resolved force list says is covered by unit tests in `recipes`; this covers +//! only that the warning reaches a logger, which those cannot see. + +// There is no fallback to report without NRLMSISE-00 in the build, so the whole file +// compiles away rather than leaving an unused logger behind. +#![cfg(feature = "msis")] + +use astrokit::prelude::*; +use log::{Metadata, Record}; +use std::sync::{Mutex, OnceLock}; + +static LINES: OnceLock>> = OnceLock::new(); + +struct Collect; + +impl log::Log for Collect { + fn enabled(&self, _: &Metadata) -> bool { + true + } + fn log(&self, record: &Record) { + LINES + .get_or_init(Default::default) + .lock() + .expect("the collector is not poisoned") + .push(format!("{} {}", record.level(), record.args())); + } + fn flush(&self) {} +} + +#[test] +fn test_falling_back_to_the_fixed_profile_is_logged() { + log::set_boxed_logger(Box::new(Collect)).expect("no other logger in this binary"); + log::set_max_level(log::LevelFilter::Trace); + + // Inside the record: the measured model, and nothing to report. + let covered = Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0); + let _ = Regime::LowEarth.forces(Fidelity::Balanced, covered); + assert!( + LINES + .get_or_init(Default::default) + .lock() + .unwrap() + .is_empty(), + "an epoch inside the record should say nothing" + ); + + // Far enough ahead to stay outside it however fresh the embedded drivers become. + let distant = Epoch::from_gregorian_utc_hms(2100, 1, 1, 0, 0, 0); + let _ = Regime::LowEarth.forces(Fidelity::Balanced, distant); + + let lines = LINES.get_or_init(Default::default).lock().unwrap().clone(); + assert_eq!(lines.len(), 1, "exactly one warning expected: {lines:?}"); + assert!(lines[0].starts_with("WARN"), "{}", lines[0]); + // It has to name both models, or it does not tell the reader what changed. + assert!(lines[0].contains("NRLMSISE-00"), "{}", lines[0]); + assert!(lines[0].contains("Harris-Priester"), "{}", lines[0]); +}