connor/astrokit
1
0 Fork 0
Code Issues 20 Pull Requests Actions Packages Projects Releases Wiki Activity

Several fixes, plus some conveniences and sane defaults

Browse Source
This commit is contained in:
connor committed 2026-10-03 12:34:49 -04:00
1 parent 9c236d5e97
commit 82074b7c76
13 files changed
+730 -173
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

Cargo.toml
+15 -1
View File
@@ -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"
@@ -48,3 +52,13 @@ harness = false
[[bench]]
name = "propagation"
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
examples/full.rs
+12 -59
View File
@@ -4,37 +4,18 @@
//!
//! 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;
use astrokit::prelude::*;
use astrokit::{attitude, kernels, recipes};
use astrokit::{attitude, kernels};
use nalgebra::Vector3;
/// The truncation the recipe chose, so the header cannot disagree with what ran.
const DEGREE: usize = recipes::LOW_ORBIT_DEGREE;
fn main() -> Result<(), Box<dyn std::error::Error>> {
// `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;
@@ -63,8 +44,7 @@ fn main() -> Result<(), Box<dyn std::error::Error>> {
// 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 = recipes::low_earth_orbit();
let field = GravityModel::Egm96.field(DEGREE, DEGREE)?;
let forces = Regime::LowEarth.forces(Fidelity::Precise, epoch);
let mut propagation = Propagation::new(epoch, initial, &forces)
.frame(frame)
@@ -72,6 +52,11 @@ fn main() -> Result<(), Box<dyn std::error::Error>> {
.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());
@@ -96,6 +81,7 @@ fn main() -> Result<(), Box<dyn std::error::Error>> {
// 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",
@@ -134,27 +120,6 @@ fn main() -> Result<(), Box<dyn std::error::Error>> {
let inertia = &spacecraft.inertia;
let momentum = |rate: Vector3<f64>| (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,
@@ -265,17 +230,5 @@ fn main() -> Result<(), Box<dyn std::error::Error>> {
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!();
// Everything the run was given except the state: the resolved force list, the
// integrator and its tolerances, the vehicle, and whether kernels were attached.
println!("provenance");
println!(
"{}",
serde_json::to_string_pretty(&propagation.provenance())?
);
Ok(())
}
readme.md
+4 -2
View File
@@ -102,8 +102,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
src/atmosphere/mod.rs
+12
View File
@@ -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"),
}
}
}
src/atmosphere/msis.rs
+1 -1
View File
@@ -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;
src/dynamics/mod.rs
+28
View File
@@ -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::*;
src/kernels.rs
+44
View File
@@ -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<PathBuf> {
.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<Arc<Almanac>, KernelError> {
let candidates: Vec<PathBuf> = 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
src/lib.rs
+1
View File
@@ -59,6 +59,7 @@ pub mod prelude {
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,
src/propagation.rs
+13 -1
View File
@@ -523,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."
);
}
}
}
src/provenance.rs
+187 -6
View File
@@ -95,17 +95,114 @@ impl Provenance {
}
}
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<String> = 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;
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::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0),
epoch(),
Orbit
.with_position(Vector3::new(6.978e6, 0.0, 0.0))
.with_velocity(Vector3::new(0.0, 7.56e3, 0.0)),
@@ -118,7 +215,7 @@ mod tests {
/// and read back.
#[test]
fn test_it_round_trips_through_json() {
let forces = recipes::low_earth_orbit();
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)
@@ -131,7 +228,7 @@ mod tests {
/// What it actually carries: the resolved composition rather than a mode name.
#[test]
fn test_it_records_the_resolved_setup() {
let forces = recipes::low_earth_orbit();
let forces = Regime::LowEarth.forces(Fidelity::Precise, epoch());
let record = propagation(&forces).accuracy(Accuracy::High).provenance();
assert_eq!(record.astrokit, VERSION);
@@ -156,8 +253,8 @@ mod tests {
/// Two runs set up the same way record the same thing, and one difference shows.
#[test]
fn test_it_distinguishes_two_setups() {
let low = recipes::low_earth_orbit();
let geo = recipes::geostationary();
let low = Regime::LowEarth.forces(Fidelity::Precise, epoch());
let geo = Regime::Geostationary.forces(Fidelity::Precise, epoch());
assert_eq!(
propagation(&low).provenance(),
propagation(&low).provenance()
@@ -211,3 +308,87 @@ mod tests {
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"));
}
}
src/recipes.rs
+332 -103
View File
@@ -1,19 +1,40 @@
//! Force lists for common orbit regimes.
//! Force lists for common orbit regimes, at a chosen level of detail.
//!
//! Each of these returns a plain `Vec<Dynamics>`, not a mode or a handle. That is the
//! point: a recipe 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.
//! A recipe returns a plain `Vec<Dynamics>`, 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;
//! use astrokit::recipes::{Fidelity, Regime};
//!
//! let mut forces = recipes::low_earth_orbit();
//! 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
@@ -22,124 +43,218 @@
//! 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;
/// Degree and order for a low orbit, where the field is close and rough.
pub const LOW_ORBIT_DEGREE: usize = 70;
/// 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,
}
/// Degree and order for a medium orbit, around half an Earth radius up.
pub const MEDIUM_ORBIT_DEGREE: usize = 16;
/// 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,
}
/// Degree and order for a geostationary orbit, where the field is distant and smooth.
pub const GEOSTATIONARY_DEGREE: usize = 8;
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<Dynamics> {
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.
/// 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<Dynamics> {
vec![TwoBody { body }]
}
/// A low Earth orbit: a close field, the Sun and Moon, sunlight and drag.
///
/// Drag is Harris-Priester, which is always available and reads no space weather. A run
/// that wants the measured atmosphere replaces that one entry; the list is data.
pub fn low_earth_orbit() -> Vec<Dynamics> {
vec![
Gravity {
field: GravityModel::Egm96,
degree: LOW_ORBIT_DEGREE,
order: LOW_ORBIT_DEGREE,
},
ThirdBody { body: Body::MOON },
ThirdBody { body: Body::SUN },
Sunshine,
Drag {
density: DensityModel::harris_priester(),
},
]
}
/// A medium Earth orbit: a navigation constellation's altitude, where there is no
/// atmosphere worth modeling and the field has smoothed out considerably.
pub fn medium_earth_orbit() -> Vec<Dynamics> {
vec![
Gravity {
field: GravityModel::Egm96,
degree: MEDIUM_ORBIT_DEGREE,
order: MEDIUM_ORBIT_DEGREE,
},
ThirdBody { body: Body::MOON },
ThirdBody { body: Body::SUN },
Sunshine,
]
}
/// A geostationary orbit: a distant field, the Sun and Moon, and sunlight.
///
/// No drag. The third bodies matter more here than anywhere lower — they are what drives
/// the inclination a station-keeping budget is mostly spent against.
pub fn geostationary() -> Vec<Dynamics> {
vec![
Gravity {
field: GravityModel::Egm96,
degree: GEOSTATIONARY_DEGREE,
order: GEOSTATIONARY_DEGREE,
},
ThirdBody { body: Body::MOON },
ThirdBody { body: Body::SUN },
Sunshine,
]
}
#[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 (name, forces) in [
("low", low_earth_orbit()),
("medium", medium_earth_orbit()),
("geostationary", geostationary()),
] {
let fields = forces
.iter()
.filter(|force| matches!(force, Gravity { .. }))
.count();
let monopoles = forces
.iter()
.filter(|force| matches!(force, TwoBody { .. }))
.count();
assert_eq!(fields, 1, "{name}");
assert_eq!(monopoles, 0, "{name}: a point mass beside a field");
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 follows the altitude, because a term of degree n falls off as
/// (R/r)^n and there is far less of the field left to resolve further out.
/// 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_falls_with_altitude() {
let degree = |forces: Vec<Dynamics>| match forces[0] {
Gravity { degree, order, .. } => {
assert_eq!(
degree, order,
"a recipe should not truncate the order alone"
);
degree
}
ref other => panic!("gravity should lead the list, found {other:?}"),
};
assert!(degree(low_earth_orbit()) > degree(medium_earth_orbit()));
assert!(degree(medium_earth_orbit()) > degree(geostationary()));
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));
}
/// Drag belongs to the one regime with an atmosphere worth the name.
/// 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<Dynamics>| {
@@ -148,27 +263,141 @@ mod tests {
.filter(|force| matches!(force, Drag { .. }))
.count()
};
assert_eq!(drags(low_earth_orbit()), 1);
assert_eq!(drags(medium_earth_orbit()), 0);
assert_eq!(drags(geostationary()), 0);
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 = low_earth_orbit();
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 text = serde_json::to_string(&low_earth_orbit()).unwrap();
let recipe = Regime::LowEarth.forces(Fidelity::Precise, epoch());
let text = serde_json::to_string(&recipe).unwrap();
assert_eq!(
serde_json::from_str::<Vec<Dynamics>>(&text).unwrap(),
low_earth_orbit()
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::<Vec<_>>()
.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 }]);
src/spacecraft.rs
+16
View File
@@ -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::*;
tests/density_fallback_warning.rs
+65
View File
@@ -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<Mutex<Vec<String>>> = 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]);
}