This commit is contained in:
connor committed 2026-10-03 12:49:20 -04:00
1 parent 3331531d4a
commit eade33bd6c
5 files changed
+42 -37

No files matched your search

+1 -4
View File
@@ -17,7 +17,6 @@
//! Orekit does and this does not, and it is arguably the better treatment, so the
//! radiation row is not purely a speed difference.
use std::path::Path;
use std::sync::Arc;
use astrokit::dynamics::{self, Context};
@@ -31,9 +30,7 @@ use std::hint::black_box;
const ARC: f64 = 86400.0;
fn almanac() -> Option<Arc<anise::almanac::Almanac>> {
let directory = kernels::directory()
.unwrap_or_else(|| Path::new(env!("CARGO_MANIFEST_DIR")).join("data/kernels"));
kernels::load(&directory).ok().map(Arc::new)
kernels::load_default().ok()
}
fn bench(c: &mut Criterion) {
-5
View File
@@ -230,11 +230,6 @@ mod tests {
// approximating.
assert_relative_eq!(node_rate, expected_node, max_relative = 1e-2);
assert_relative_eq!(periapsis_rate, expected_periapsis, max_relative = 1e-2);
println!(
"J2 secular rates: node {:.2e}, periapsis {:.2e} relative",
(node_rate - expected_node).abs() / expected_node.abs(),
(periapsis_rate - expected_periapsis).abs() / expected_periapsis.abs()
);
// And the signs, which first-order theory gets exactly right.
assert!(node_rate < 0.0, "a prograde node should regress");
+3 -8
View File
@@ -653,12 +653,12 @@ mod tests {
let position_error = (end.position() - leo().position()).norm();
let velocity_error = (end.velocity() - leo().velocity()).norm();
// Measured 2026-10-03: 9.8e-7 m and 1.1e-9 m/s, an order inside both bounds.
assert!(position_error < 1e-5, "position error {position_error:e} m");
assert!(
velocity_error < 1e-8,
"velocity error {velocity_error:e} m/s"
);
println!("one-period closure: {position_error:.2e} m, {velocity_error:.2e} m/s");
}
/// Two-body motion conserves specific energy and angular momentum. Both are
@@ -675,20 +675,15 @@ mod tests {
|y: &State<f64, 6, Orbit>| y.velocity().norm_squared() * 0.5 - MU / y.position().norm();
let momentum = |y: &State<f64, 6, Orbit>| y.position().cross(&y.velocity());
// Measured 2026-10-03: energy drifts 1.1e-13 relative over five periods and
// angular momentum 5.3e-14, both an order inside the bound below.
let e0 = energy(&leo());
let h0 = momentum(&leo());
let (mut worst_energy, mut worst_momentum) = (0.0_f64, 0.0_f64);
for state in &solution.states {
let y = Orbit.state(*state);
assert_relative_eq!(energy(&y), e0, max_relative = 1e-12);
assert_relative_eq!(momentum(&y).norm(), h0.norm(), max_relative = 1e-12);
worst_energy = worst_energy.max((energy(&y) - e0).abs() / e0.abs());
worst_momentum =
worst_momentum.max((momentum(&y).norm() - h0.norm()).abs() / h0.norm());
}
println!(
"five-period drift: energy {worst_energy:.2e}, angular momentum {worst_momentum:.2e} relative"
);
}
/// The same spec and the same dynamics, propagated on three layouts of three
+25 -7
View File
@@ -29,9 +29,7 @@ use std::sync::Arc;
/// The kernels, or `None` if they have not been fetched.
fn almanac() -> Option<Arc<anise::almanac::Almanac>> {
let directory = kernels::directory()
.unwrap_or_else(|| Path::new(env!("CARGO_MANIFEST_DIR")).join("data/kernels"));
kernels::load(&directory).ok().map(Arc::new)
kernels::load_default().ok()
}
/// Loads a fixture and checks it is the one it claims to be.
@@ -97,6 +95,7 @@ fn test_a_point_mass_arc_agrees_with_nyx() {
let (position, velocity) = compare(&arc, &[TwoBody { body: Body::EARTH }], None, None);
// Measured 2026-10-03: 2.18e-4 m and 2.36e-7 m/s.
assert!(
position < 1.0e-3,
"worst position disagreement {position:.3e} m over the day"
@@ -105,7 +104,6 @@ fn test_a_point_mass_arc_agrees_with_nyx() {
velocity < 1.0e-6,
"worst velocity disagreement {velocity:.3e} m/s over the day"
);
println!("point mass: {position:.3e} m, {velocity:.3e} m/s");
}
/// The Sun and the Moon over a point mass, from the same DE-440 through the same ANISE.
@@ -129,11 +127,15 @@ fn test_the_sun_and_moon_agree_with_nyx() {
Some(almanac),
);
// Measured 2026-10-03: 1.97e-4 m and 2.13e-7 m/s.
assert!(
position < 1.0e-3,
"worst position disagreement {position:.3e} m over the day"
);
println!("third body: {position:.3e} m, {velocity:.3e} m/s");
assert!(
velocity < 1.0e-6,
"worst velocity disagreement {velocity:.3e} m/s over the day"
);
}
/// EGM96, read by both sides from the same GMAT coefficient file, and turned with the
@@ -164,16 +166,22 @@ fn gravity(degree: usize, allowed_meters: f64) {
position < allowed_meters,
"degree {degree}: worst position disagreement {position:.3e} m over the day"
);
println!("gravity {degree}: {position:.3e} m, {velocity:.3e} m/s");
assert!(
velocity < 1.0e-6,
"degree {degree}: worst velocity disagreement {velocity:.3e} m/s over the day"
);
}
#[test]
fn test_a_low_degree_field_agrees_with_nyx() {
// Measured 2026-10-03: 1.02e-4 m, against the one meter a day this is aiming at.
gravity(8, 0.03);
}
#[test]
fn test_a_higher_degree_field_agrees_with_nyx() {
// Measured 2026-10-03: 2.06e-5 m. Tighter than degree 8 rather than looser, which
// is the step placement differing, not the field.
gravity(20, 0.03);
}
@@ -194,5 +202,15 @@ fn test_radiation_pressure_agrees_with_nyx() {
Some(almanac),
);
println!("radiation: {position:.3e} m, {velocity:.3e} m/s");
// Measured 2026-10-03: 4.91e-4 m and 5.32e-7 m/s. This had no assertion at all
// when it was written, so it could not have failed; the figures below are the
// same order as the other matched scenarios.
assert!(
position < 1.0e-2,
"worst position disagreement {position:.3e} m over the day"
);
assert!(
velocity < 1.0e-5,
"worst velocity disagreement {velocity:.3e} m/s over the day"
);
}
+13 -13
View File
@@ -20,9 +20,7 @@ use std::sync::Arc;
/// The kernels, or `None` if they have not been fetched.
fn almanac() -> Option<Arc<anise::almanac::Almanac>> {
let directory = kernels::directory()
.unwrap_or_else(|| Path::new(env!("CARGO_MANIFEST_DIR")).join("data/kernels"));
kernels::load(&directory).ok().map(Arc::new)
kernels::load_default().ok()
}
/// Propagates the fixture's initial state with this crate and returns the worst position
@@ -82,7 +80,6 @@ fn test_a_point_mass_arc_agrees_with_orekit() {
velocity < 1.0e-6,
"worst velocity disagreement {velocity:.3e} m/s over the day"
);
println!("point mass: {position:.3e} m, {velocity:.3e} m/s");
}
/// EGM96, read by both sides from the same coefficient file.
@@ -113,12 +110,15 @@ fn gravity(degree: usize, allowed_meters: f64) {
position < allowed_meters,
"degree {degree}: worst position disagreement {position:.3e} m over the day"
);
println!("gravity {degree}: {position:.3e} m, {velocity:.3e} m/s");
assert!(
velocity < 1.0e-4,
"degree {degree}: worst velocity disagreement {velocity:.3e} m/s over the day"
);
}
#[test]
fn test_a_low_degree_field_agrees_with_orekit() {
// Measured at 1.1 cm, against the one meter a day this is aiming at.
// Measured at 1.089 cm, against the one meter a day this is aiming at.
gravity(8, 0.03);
}
@@ -154,11 +154,16 @@ fn test_the_sun_and_moon_agree_with_orekit() {
],
Some(almanac),
);
// Measured 2026-10-03: 8.36e-4 m, so adding the Sun and Moon did not move the
// disagreement off the point-mass floor at all.
assert!(
position < 0.1,
position < 1.0e-3,
"worst position disagreement {position:.3e} m over the day"
);
println!("third body: {position:.3e} m, {velocity:.3e} m/s");
assert!(
velocity < 1.0e-6,
"worst velocity disagreement {velocity:.3e} m/s over the day"
);
}
/// Radiation pressure on a cannonball, over a point mass.
@@ -205,7 +210,6 @@ fn test_radiation_pressure_agrees_with_orekit() {
velocity =
velocity.max((ours.velocity() - Vector3::from_column_slice(&state.velocity)).norm());
}
println!("radiation: {position:.3e} m, {velocity:.3e} m/s");
assert!(
position < 0.01,
"worst position disagreement {position:.3e} m over the day"
@@ -272,10 +276,6 @@ fn test_msis_drag_agrees_with_orekit() {
let last = arc.states.last().expect("checked on load");
let effect = (Vector3::from_column_slice(&last.position) - drag_free).norm();
println!(
"msis drag: {position:.3e} m, {velocity:.3e} m/s, {:.4}% of a {effect:.0} m effect",
position / effect * 100.0
);
assert!(
position < 1.0e-5 * effect,
"worst disagreement {position:.3e} m is {:.4}% of the {effect:.0} m drag moves it",