285 lines
11 KiB
Rust
285 lines
11 KiB
Rust
//! This crate's propagator against Orekit's, on the fixtures in `tests/fixtures/orekit`.
|
|
//!
|
|
//! Generated by `scripts/generate-orekit-fixtures.sh`. Each fixture records the setup
|
|
//! Orekit was given, so a disagreement can be attributed to the models rather than left
|
|
//! as a mystery; see `tests/common/mod.rs` for what has to match and why.
|
|
//!
|
|
//! The point-mass fixture comes first and is not about gravity at all: with no field, no
|
|
//! third body, no atmosphere and an inertial frame, there is nothing left to mismatch, so
|
|
//! the figure it produces measures the comparison itself — frames, time scale, units and
|
|
//! `mu`. Every later scenario is read against that floor.
|
|
|
|
mod common;
|
|
|
|
use astrokit::kernels;
|
|
use astrokit::prelude::*;
|
|
use common::ReferenceArc;
|
|
use nalgebra::Vector3;
|
|
use std::path::Path;
|
|
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)
|
|
}
|
|
|
|
/// Propagates the fixture's initial state with this crate and returns the worst position
|
|
/// and velocity disagreement over the arc.
|
|
fn compare(
|
|
arc: &ReferenceArc,
|
|
forces: &[Dynamics],
|
|
almanac: Option<Arc<anise::almanac::Almanac>>,
|
|
) -> (f64, f64) {
|
|
let start = &arc.states[0];
|
|
let epoch = Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0);
|
|
let mut propagation = Propagation::new(
|
|
epoch,
|
|
Orbit
|
|
.with_position(Vector3::from_column_slice(&start.position))
|
|
.with_velocity(Vector3::from_column_slice(&start.velocity)),
|
|
forces,
|
|
)
|
|
.frame(Frame::EARTH_J2000)
|
|
.duration(Duration::from_seconds(arc.duration()))
|
|
.accuracy(Accuracy::High);
|
|
if let Some(almanac) = almanac {
|
|
propagation = propagation.almanac(almanac);
|
|
}
|
|
|
|
let solution = propagation.run();
|
|
let (mut worst_position, mut worst_velocity) = (0.0_f64, 0.0_f64);
|
|
for state in &arc.states {
|
|
let ours = Orbit.state(solution.interpolate(state.elapsed));
|
|
worst_position = worst_position
|
|
.max((ours.position() - Vector3::from_column_slice(&state.position)).norm());
|
|
worst_velocity = worst_velocity
|
|
.max((ours.velocity() - Vector3::from_column_slice(&state.velocity)).norm());
|
|
}
|
|
(worst_position, worst_velocity)
|
|
}
|
|
|
|
#[test]
|
|
fn test_a_point_mass_arc_agrees_with_orekit() {
|
|
let path = Path::new("tests/fixtures/orekit/two-body-low-earth.json");
|
|
let arc = ReferenceArc::load(path).expect("the fixture is committed");
|
|
|
|
assert_eq!(arc.setup.tool, "orekit");
|
|
assert_eq!(arc.duration(), 86400.0);
|
|
assert_eq!(arc.states.len(), 1441);
|
|
|
|
let (position, velocity) = compare(&arc, &[TwoBody { body: Body::EARTH }], None);
|
|
|
|
// Sub-millimeter over a day. Nothing in this case can disagree except the two
|
|
// integrators and the arithmetic, so this is the floor the comparison can measure to:
|
|
// a later fixture that disagrees by more than this is disagreeing about physics.
|
|
assert!(
|
|
position < 1.0e-3,
|
|
"worst position disagreement {position:.3e} m over the day"
|
|
);
|
|
assert!(
|
|
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.
|
|
///
|
|
/// The frames are matched as closely as the two tools allow: GCRF on the inertial side,
|
|
/// which is ICRF-aligned like the frame this crate gets from ANISE, and ITRF93 on the
|
|
/// body-fixed side, which is the realization its SPICE kernel carries. What is left is
|
|
/// about 4 mas between the two rotations, and the centimeter below is mostly that.
|
|
fn gravity(degree: usize, allowed_meters: f64) {
|
|
let Some(almanac) = almanac() else {
|
|
eprintln!("skipped: kernels not fetched, run scripts/fetch-kernels.sh");
|
|
return;
|
|
};
|
|
let path = format!("tests/fixtures/orekit/gravity-{degree}-low-earth.json");
|
|
let arc = ReferenceArc::load(Path::new(&path)).expect("the fixture is committed");
|
|
assert!(arc.setup.forces[0].contains("EGM96"));
|
|
|
|
let (position, velocity) = compare(
|
|
&arc,
|
|
&[Gravity {
|
|
field: GravityModel::Egm96,
|
|
degree,
|
|
order: degree,
|
|
}],
|
|
Some(almanac),
|
|
);
|
|
assert!(
|
|
position < allowed_meters,
|
|
"degree {degree}: worst position disagreement {position:.3e} m over the day"
|
|
);
|
|
println!("gravity {degree}: {position:.3e} m, {velocity:.3e} m/s");
|
|
}
|
|
|
|
#[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.
|
|
gravity(8, 0.03);
|
|
}
|
|
|
|
#[test]
|
|
fn test_a_higher_degree_field_agrees_with_orekit() {
|
|
// Degree 20 rather than 70: the agreement is the same to within a couple of
|
|
// centimeters and the arc costs a fraction of the time in a debug build.
|
|
gravity(20, 0.03);
|
|
}
|
|
|
|
/// The Sun and the Moon over a point mass, from the same DE-440 ephemeris on both sides.
|
|
///
|
|
/// Isolated from the field on purpose: the third-body term is the difference of two
|
|
/// nearly equal attractions, so it is the one most easily got wrong, and a disagreement
|
|
/// here would otherwise hide under the gravity figure.
|
|
#[test]
|
|
fn test_the_sun_and_moon_agree_with_orekit() {
|
|
let Some(almanac) = almanac() else {
|
|
eprintln!("skipped: kernels not fetched, run scripts/fetch-kernels.sh");
|
|
return;
|
|
};
|
|
let path = Path::new("tests/fixtures/orekit/third-body-low-earth.json");
|
|
let arc = ReferenceArc::load(path).expect("the fixture is committed");
|
|
assert!(arc.setup.forces.iter().any(|f| f.contains("Moon")));
|
|
assert!(arc.setup.forces.iter().any(|f| f.contains("Sun")));
|
|
|
|
let (position, velocity) = compare(
|
|
&arc,
|
|
&[
|
|
TwoBody { body: Body::EARTH },
|
|
ThirdBody { body: Body::MOON },
|
|
ThirdBody { body: Body::SUN },
|
|
],
|
|
Some(almanac),
|
|
);
|
|
assert!(
|
|
position < 0.1,
|
|
"worst position disagreement {position:.3e} m over the day"
|
|
);
|
|
println!("third body: {position:.3e} m, {velocity:.3e} m/s");
|
|
}
|
|
|
|
/// Radiation pressure on a cannonball, over a point mass.
|
|
///
|
|
/// The matched knobs are the reference pressure (this crate's 1361 W/m^2 rather than
|
|
/// Orekit's 1367), the cross-section and coefficient, the mass, and a *spherical*
|
|
/// occulting Earth, since this crate's eclipse geometry treats the occulter as a sphere.
|
|
/// What is left unmatched is the solar radius, 6.957e8 here against Orekit's 6.95508e8,
|
|
/// which changes only how wide the penumbra is.
|
|
#[test]
|
|
fn test_radiation_pressure_agrees_with_orekit() {
|
|
let Some(almanac) = almanac() else {
|
|
eprintln!("skipped: kernels not fetched, run scripts/fetch-kernels.sh");
|
|
return;
|
|
};
|
|
let path = Path::new("tests/fixtures/orekit/radiation-low-earth.json");
|
|
let arc = ReferenceArc::load(path).expect("the fixture is committed");
|
|
assert!(arc
|
|
.setup
|
|
.forces
|
|
.iter()
|
|
.any(|f| f.contains("SolarRadiationPressure")));
|
|
|
|
let start = &arc.states[0];
|
|
let mut propagation = Propagation::new(
|
|
Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0),
|
|
Orbit
|
|
.with_position(Vector3::from_column_slice(&start.position))
|
|
.with_velocity(Vector3::from_column_slice(&start.velocity)),
|
|
&[TwoBody { body: Body::EARTH }, Sunshine],
|
|
)
|
|
.frame(Frame::EARTH_J2000)
|
|
.spacecraft(Spacecraft::cannonball(500.0, 5.0, 2.2, 1.3))
|
|
.almanac(almanac)
|
|
.duration(Duration::from_seconds(arc.duration()))
|
|
.accuracy(Accuracy::High);
|
|
|
|
let solution = propagation.run();
|
|
let (mut position, mut velocity) = (0.0_f64, 0.0_f64);
|
|
for state in &arc.states {
|
|
let ours = Orbit.state(solution.interpolate(state.elapsed));
|
|
position =
|
|
position.max((ours.position() - Vector3::from_column_slice(&state.position)).norm());
|
|
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"
|
|
);
|
|
}
|
|
|
|
/// Drag through NRLMSISE-00, which both tools implement from the same reference.
|
|
///
|
|
/// Nothing in the model is intended to differ here, so this is the measurement of how
|
|
/// close two independent implementations get when they are given very nearly the same
|
|
/// weather. The drivers are not quite identical: for this epoch the
|
|
/// daily flux and every three-hourly ap agree exactly, but CSSI publishes the 81-day
|
|
/// average rounded to a tenth, 154.3 against the 154.3333 this crate computes from the
|
|
/// daily values.
|
|
#[cfg(feature = "msis")]
|
|
#[test]
|
|
fn test_msis_drag_agrees_with_orekit() {
|
|
let Some(almanac) = almanac() else {
|
|
eprintln!("skipped: kernels not fetched, run scripts/fetch-kernels.sh");
|
|
return;
|
|
};
|
|
let path = Path::new("tests/fixtures/orekit/msis-drag-low-earth.json");
|
|
let arc = ReferenceArc::load(path).expect("the fixture is committed");
|
|
assert!(arc.setup.forces.iter().any(|f| f.contains("NRLMSISE00")));
|
|
|
|
let start = &arc.states[0];
|
|
let initial = Orbit
|
|
.with_position(Vector3::from_column_slice(&start.position))
|
|
.with_velocity(Vector3::from_column_slice(&start.velocity));
|
|
let epoch = Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0);
|
|
let mut propagation = Propagation::new(
|
|
epoch,
|
|
initial,
|
|
&[
|
|
TwoBody { body: Body::EARTH },
|
|
Drag {
|
|
density: DensityModel::Nrlmsise00,
|
|
},
|
|
],
|
|
)
|
|
.frame(Frame::EARTH_J2000)
|
|
.spacecraft(Spacecraft::cannonball(500.0, 5.0, 2.2, 1.3))
|
|
.almanac(almanac)
|
|
.duration(Duration::from_seconds(arc.duration()))
|
|
.accuracy(Accuracy::High);
|
|
|
|
let solution = propagation.run();
|
|
let (mut position, mut velocity) = (0.0_f64, 0.0_f64);
|
|
for state in &arc.states {
|
|
let ours = Orbit.state(solution.interpolate(state.elapsed));
|
|
position =
|
|
position.max((ours.position() - Vector3::from_column_slice(&state.position)).norm());
|
|
velocity =
|
|
velocity.max((ours.velocity() - Vector3::from_column_slice(&state.velocity)).norm());
|
|
}
|
|
|
|
let drag_free = {
|
|
let mut plain = Propagation::new(epoch, initial, &[TwoBody { body: Body::EARTH }])
|
|
.frame(Frame::EARTH_J2000)
|
|
.duration(Duration::from_seconds(arc.duration()))
|
|
.accuracy(Accuracy::High);
|
|
Orbit.state(*plain.run().states.last().unwrap()).position()
|
|
};
|
|
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",
|
|
position / effect * 100.0
|
|
);
|
|
}
|