Added the MSIS scenario to nyx comparison
This commit is contained in:
connor committed
1 parent
eade33bd6c
commit
92a906a693
4 files changed
+217
-5
No files matched your search
@@ -109,11 +109,27 @@ is its accuracy floor. On the point-mass arc its stored states agree with Kepler
|
||||
| EGM96 8x8 | 1 m, 1 mm/s | 0.102 mm, 1.1e-7 m/s | `low_degree_field` |
|
||||
| EGM96 20x20 | 1 m, 1 mm/s | 0.021 mm, 2.2e-8 m/s | `higher_degree_field` |
|
||||
| Cannonball radiation pressure | 1 m, 1 mm/s | 0.491 mm, 5.3e-7 m/s | `radiation_pressure` |
|
||||
| Drag, NRLMSISE-00 | 1% of the drag effect | 8.55 m, 9.7e-3 m/s — 0.029% | `msis_drag` |
|
||||
|
||||
Tests are in `tests/nyx_comparison.rs`, fixtures from
|
||||
`scripts/generate-nyx-fixtures.sh`. nyx is AGPL-3.0-or-later, so it is run offline to
|
||||
produce data and is not a dependency of anything here.
|
||||
|
||||
The drag row is the only one where the two sides disagree by more than a millimeter, and
|
||||
the reason is known. Both tools are given numerically identical space weather — the
|
||||
fixture's drivers are written out of the same record this crate reads, and nyx takes the
|
||||
81-day centered average from a file column, so it carries this crate's own value rather
|
||||
than a separately published one. With the inputs identical the densities still differ by
|
||||
2.6e-4. Run against case 0 of the seventeen standard NRLMSISE-00 cases, which pin the
|
||||
NRL reference to sixteen digits, this crate reproduces the published total mass density
|
||||
exactly and nyx is 1.08e-3 from it. So that 0.029% is nyx departing from the reference.
|
||||
|
||||
One element is unmatched and it is this crate's side: nyx measures geodetic altitude
|
||||
against the planetary constants kernel's Earth and this crate against WGS84, putting the
|
||||
same spacecraft 0.400 m apart. That is worth about 7e-6 of the density, two orders below
|
||||
the model difference, and it cannot be removed from the nyx side — ANISE replaces a
|
||||
frame's ellipsoid from the loaded kernel.
|
||||
|
||||
### Analytic invariants
|
||||
|
||||
| Case | Target | Achieved | Test |
|
||||
@@ -146,13 +162,13 @@ the minimum of thirty runs is reported rather than the mean.
|
||||
| EGM96 20x20 | 2473 steps, 54.5 ms | 2158 steps, 193.6 ms | 1940 steps, 52.9 ms |
|
||||
| Sun and Moon | 1794 steps, 35.9 ms | 2152 steps, 79.7 ms | 1932 steps, 36.1 ms |
|
||||
| Radiation pressure | 2067 steps, 28.0 ms | 2570 steps, 166.9 ms | 2035 steps, 133.6 ms |
|
||||
| Drag, NRLMSISE-00 | 1859 steps, 52.0 ms | not yet measured | 2092 steps, 357.2 ms |
|
||||
| Drag, NRLMSISE-00 | 1859 steps, 52.0 ms | 2452 steps, 385.0 ms | 2092 steps, 357.2 ms |
|
||||
| Ensemble of N states, scaling vs. cores | not yet measured | — | — |
|
||||
|
||||
Step counts are reported alongside the times because the tolerances are not literally
|
||||
equal: Orekit is asked for a nanometer of position accuracy and this runs at
|
||||
`Accuracy::High`. The two do comparable work — within a quarter of each other on every
|
||||
row — so a ratio of times is close to a ratio of per-step cost. A future scenario whose
|
||||
`Accuracy::High`. The three do comparable work — within about a third of each other on
|
||||
every row — so a ratio of times is close to a ratio of per-step cost. A future scenario whose
|
||||
step counts diverge needs them brought back together before its ratio means anything.
|
||||
|
||||
The radiation row is not like for like however the steps come out. Orekit's radiation
|
||||
|
||||
@@ -79,6 +79,9 @@ hifitime = "4.3"
|
||||
nalgebra = "0.35"
|
||||
serde = { version = "1.0", features = ["derive"] }
|
||||
serde_json = "1.0"
|
||||
# Our own crate, MIT. The generator writes nyx's space weather out of the same embedded
|
||||
# record astrokit reads, so both tools see identical drivers.
|
||||
nrlmsise00 = { version = "0.2", features = ["drivers"] }
|
||||
TOML
|
||||
|
||||
cat > "$work/gen/src/main.rs" <<'RUSTEOF'
|
||||
@@ -89,9 +92,12 @@ use std::sync::Arc;
|
||||
use anise::constants::celestial_objects::{MOON, SUN};
|
||||
use anise::constants::frames::{EARTH_ITRF93, EARTH_J2000};
|
||||
use anise::prelude::{Almanac, Orbit};
|
||||
use anise::structure::spacecraft::{Mass, SRPData};
|
||||
use anise::structure::spacecraft::{DragData, Mass, SRPData};
|
||||
use hifitime::{Duration, Epoch};
|
||||
use nyx_space::dynamics::drag::nrlmsise00::{GeomagneticMode, Nrlmsise00Flags};
|
||||
use nyx_space::dynamics::drag::{AtmDensity, Drag};
|
||||
use nyx_space::dynamics::{GravityField, OrbitalDynamics, SolarPressure, SpacecraftDynamics};
|
||||
use nyx_space::io::space_weather::{SpaceWeatherData, StaticSpaceWeather};
|
||||
use nyx_space::io::gravity::GravityFieldData;
|
||||
use nyx_space::propagators::{IntegratorOptions, Propagator};
|
||||
use nyx_space::Spacecraft;
|
||||
@@ -145,6 +151,52 @@ rather than the integration: on the point-mass arc at this tolerance the stored
|
||||
agree with Kepler to 2.7e-6 m while reading the same trajectory on a 60 s grid is 8.1e-4 m. \
|
||||
Sampling the steps themselves removes that term from the comparison.";
|
||||
|
||||
|
||||
/// Writes the space weather nyx reads, from the record astrokit reads.
|
||||
///
|
||||
/// nyx wants a CelesTrak-shaped CSV and takes the 81-day centered average straight out
|
||||
/// of a column rather than computing it. Both of those are useful here: emitting the
|
||||
/// file from 's own embedded drivers means the two tools see numerically
|
||||
/// identical weather, where the Orekit comparison had to reconcile two published
|
||||
/// records and still lost 0.028% to CSSI rounding its average to a tenth.
|
||||
///
|
||||
/// Covers the arc plus a margin either side, because nyx reaches three days back to
|
||||
/// build the 57-hour ap history.
|
||||
fn write_space_weather(path: &str, epoch: Epoch) {
|
||||
use std::fmt::Write;
|
||||
|
||||
let drivers = nrlmsise00::drivers::Drivers::embedded();
|
||||
let start = epoch.to_mjd_utc_days().floor() as i64 - MARGIN_DAYS;
|
||||
let end = (epoch.to_mjd_utc_days() + ARC_SECONDS / 86400.0).floor() as i64 + MARGIN_DAYS;
|
||||
|
||||
let mut out = String::from(
|
||||
"DATE,BSRN,ND,KP1,KP2,KP3,KP4,KP5,KP6,KP7,KP8,KP_SUM, AP1,AP2,AP3,AP4,AP5,AP6,AP7,AP8,AP_AVG,CP,C9,ISN, F10.7_OBS,F10.7_ADJ,F10.7_DATA_TYPE,F10.7_OBS_CENTER81,F10.7_OBS_LAST81, F10.7_ADJ_CENTER81,F10.7_ADJ_LAST81
|
||||
",
|
||||
);
|
||||
for mjd in start..=end {
|
||||
let day = drivers.day(mjd as f64).expect("inside the record");
|
||||
let centered = drivers.centered_flux(mjd as f64).expect("inside the record");
|
||||
let date = Epoch::from_mjd_utc(mjd as f64).to_gregorian_str(hifitime::TimeScale::UTC);
|
||||
let date = &date[..10];
|
||||
write!(out, "{date},0,0,,,,,,,,,,").expect("string");
|
||||
for value in day.ap {
|
||||
write!(out, "{value},").expect("string");
|
||||
}
|
||||
write!(out, "{},,,,", day.ap_daily).expect("string");
|
||||
writeln!(
|
||||
out,
|
||||
"{},{},OBS,{centered},{centered},{centered},{centered}",
|
||||
day.flux, day.flux
|
||||
)
|
||||
.expect("string");
|
||||
}
|
||||
std::fs::write(path, out).expect("write space weather");
|
||||
}
|
||||
|
||||
/// nyx reaches three days back for the ap history; four is margin on that.
|
||||
const MARGIN_DAYS: i64 = 4;
|
||||
|
||||
|
||||
fn main() {
|
||||
let args: Vec<String> = std::env::args().collect();
|
||||
let (kernels, out_dir, scenario) = (&args[1], &args[2], args[3].as_str());
|
||||
@@ -166,6 +218,8 @@ fn main() {
|
||||
let epoch = Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0);
|
||||
// The same initial conditions the Orekit fixtures use.
|
||||
let high = Orbit::new(6978.0, 0.0, 0.0, 0.0, 7.560, 0.0, epoch, earth);
|
||||
// Drag is measured lower down, where there is an atmosphere worth integrating.
|
||||
let low = Orbit::new(6778.0, 0.0, 0.0, 0.0, 7.668, 0.0, epoch, earth);
|
||||
|
||||
let (orbit, dynamics, forces, mut notes) = match scenario {
|
||||
"two-body" => (
|
||||
@@ -249,6 +303,68 @@ fn main() {
|
||||
],
|
||||
)
|
||||
}
|
||||
"drag-msis" => {
|
||||
let fixed = almanac.frame_info(EARTH_ITRF93).expect("itrf93");
|
||||
// A temporary file rather than one beside the fixtures: nyx reads its
|
||||
// weather from a path, and `--time` has no output directory to put it in.
|
||||
let csv = std::env::temp_dir()
|
||||
.join("astrokit-nyx-space-weather.csv")
|
||||
.to_string_lossy()
|
||||
.into_owned();
|
||||
write_space_weather(&csv, epoch);
|
||||
let weather = SpaceWeatherData::from_csv_file(&csv, StaticSpaceWeather::default())
|
||||
.expect("space weather");
|
||||
std::fs::remove_file(&csv).ok();
|
||||
let drag = Arc::new(Drag {
|
||||
density: AtmDensity::NRLMSISE00 {
|
||||
weather,
|
||||
flags: Some(Nrlmsise00Flags {
|
||||
// astrokit passes the seven-element ap history the model also
|
||||
// accepts; nyx's default is the single daily Ap. This is the
|
||||
// same switch that was worth 1.1% against Orekit.
|
||||
geomagnetic: GeomagneticMode::ExtendedHistory57h,
|
||||
..Default::default()
|
||||
}),
|
||||
},
|
||||
frame: fixed,
|
||||
estimate: false,
|
||||
});
|
||||
(
|
||||
low,
|
||||
SpacecraftDynamics::from_model(OrbitalDynamics::two_body(), drag),
|
||||
vec![
|
||||
format!("Earth point mass, mu = {EGM96_MU_KM3_S2} km^3/s^2 from the frame"),
|
||||
"Drag(AtmDensity::NRLMSISE00, ExtendedHistory57h, mean local solar time), \
|
||||
area 5 m^2, Cd 2.2"
|
||||
.to_string(),
|
||||
],
|
||||
vec![
|
||||
"Spacecraft mass 500 kg, starting at 400 km.".to_string(),
|
||||
"The space weather is written out of the same embedded GFZ record \
|
||||
astrokit reads, into the CelesTrak-shaped CSV nyx wants, so the two \
|
||||
tools see numerically identical drivers. nyx takes the 81-day \
|
||||
centered average from a column rather than computing it, so that \
|
||||
column carries astrokit's own value -- which removes the 0.028% that \
|
||||
CSSI's rounding cost the equivalent Orekit comparison."
|
||||
.to_string(),
|
||||
"Geomagnetic mode is the 57-hour ap history rather than nyx's default \
|
||||
single daily Ap, matching what astrokit passes. Local solar time is \
|
||||
mean on both sides, which is already nyx's default and was not \
|
||||
Orekit's."
|
||||
.to_string(),
|
||||
"Altitude is geodetic on both sides and the relative wind is taken in \
|
||||
EARTH_ITRF93, but the two ellipsoids differ and could not be made to \
|
||||
agree: nyx measures against the planetary constants kernel's Earth, \
|
||||
6378.1366 km equatorial, and astrokit against WGS84's 6378.137. That \
|
||||
places the same spacecraft 0.400 m apart in altitude, worth about \
|
||||
7e-6 of the density at 400 km, which is two orders below the model \
|
||||
difference below. Setting the frame's ellipsoid does not help -- \
|
||||
ANISE's transform replaces it from the loaded kernel."
|
||||
.to_string(),
|
||||
"No field, no third body, no radiation pressure.".to_string(),
|
||||
],
|
||||
)
|
||||
}
|
||||
other => panic!("unknown scenario {other}"),
|
||||
};
|
||||
|
||||
@@ -259,6 +375,10 @@ fn main() {
|
||||
area_m2: 5.0,
|
||||
coeff_reflectivity: 1.3,
|
||||
})
|
||||
.drag(DragData {
|
||||
area_m2: 5.0,
|
||||
coeff_drag: 2.2,
|
||||
})
|
||||
.build();
|
||||
|
||||
let options = IntegratorOptions {
|
||||
@@ -371,7 +491,7 @@ echo "building against nyx-space $NYX"
|
||||
( cd "$work/gen" && cargo build --release --quiet )
|
||||
gen="$work/gen/target/release/nyxgen"
|
||||
|
||||
scenarios="two-body gravity-8 gravity-20 third-body radiation"
|
||||
scenarios="two-body gravity-8 gravity-20 third-body radiation drag-msis"
|
||||
|
||||
if [ "$mode" = "time" ]; then
|
||||
for scenario in $scenarios; do
|
||||
|
||||
+1
@@ -0,0 +1 @@
|
||||
{"name":"drag-msis-low-earth","setup":{"tool":"nyx-space","version":"2.6.0","generated":"2026-10-02","forces":["Earth point mass, mu = 398600.4415 km^3/s^2 from the frame","Drag(AtmDensity::NRLMSISE00, ExtendedHistory57h, mean local solar time), area 5 m^2, Cd 2.2"],"integrator":"RK89, IntegratorOptions tolerance 1e-14, error control RSSCartesianStep, max step 2700 s","notes":["Spacecraft mass 500 kg, starting at 400 km.","The space weather is written out of the same embedded GFZ record astrokit reads, into the CelesTrak-shaped CSV nyx wants, so the two tools see numerically identical drivers. nyx takes the 81-day centered average from a column rather than computing it, so that column carries astrokit's own value -- which removes the 0.028% that CSSI's rounding cost the equivalent Orekit comparison.","Geomagnetic mode is the 57-hour ap history rather than nyx's default single daily Ap, matching what astrokit passes. Local solar time is mean on both sides, which is already nyx's default and was not Orekit's.","Altitude is geodetic on both sides and the relative wind is taken in EARTH_ITRF93, but the two ellipsoids differ and could not be made to agree: nyx measures against the planetary constants kernel's Earth, 6378.1366 km equatorial, and astrokit against WGS84's 6378.137. That places the same spacecraft 0.400 m apart in altitude, worth about 7e-6 of the density at 400 km, which is two orders below the model difference below. Setting the frame's ellipsoid does not help -- ANISE's transform replaces it from the loaded kernel.","No field, no third body, no radiation pressure.","States are nyx's own accepted integrator steps rather than a fixed grid, so nothing on this side is interpolated. nyx's trajectory interpolation is a five-point Hermite fitted after the fact over stored states, and it is the accuracy floor rather than the integration: on the point-mass arc at this tolerance the stored states agree with Kepler to 2.7e-6 m while reading the same trajectory on a 60 s grid is 8.1e-4 m. Sampling the steps themselves removes that term from the comparison."]},"epoch":"2026-01-01T00:00:00 UTC","frame":"EARTH_J2000","states":[{"elapsed":0.0,"position":[6778000.0,0.0,0.0],"velocity":[0.0,7668.0,0.0]},{"elapsed":38.947560749,"position":[6771420.469044267,298553.25204069447,9.388406105286336e-10],"velocity":[-337.81149086106865,7660.556399405986,2.1297036006966623e-10]},{"elapsed":77.31991415,"position":[6752081.460175954,592133.1783678555,2.8372554616698456e-8],"velocity":[-669.9960059433871,7638.677898282628,1.3703633673012948e-9]},{"elapsed":115.769381177,"position":[6719940.755038541,885183.4538164686,1.1788202509895291e-7],"velocity":[-1001.5818513644057,7602.316691684378,3.432552223843646e-9]},{"elapsed":154.298599822,"position":[6674978.963277502,1177161.1022591616,3.0384802000049697e-7],"velocity":[-1331.9549618880153,7551.450774966874,6.360513184628212e-9]},{"elapsed":192.651291909,"position":[6617627.741624049,1465580.6195304208,6.166910034842156e-7],"velocity":[-1658.30328113725,7486.5683935662055,1.0083864057738247e-8]},{"elapsed":231.066337677,"position":[6547694.495856405,1751702.892464092,1.0879177641160105e-6],"velocity":[-1982.0537934155443,7407.451543226486,1.457178110755673e-8]},{"elapsed":269.458332408,"position":[6465447.346956375,2034348.0254953366,1.7448337544218233e-6],"velocity":[-2301.8717137354042,7314.403410010821,1.9762845295084924e-8]},{"elapsed":307.975606211,"position":[6370675.1886256365,2314058.294870179,2.6169313246292616e-6],"velocity":[-2618.3709101507784,7207.184994543471,2.5625195607037797e-8]},{"elapsed":346.409215949,"position":[6264051.236265277,2588781.6537735225,3.723868349643218e-6],"velocity":[-2929.229693237649,7086.557732832857,3.207206790645895e-8]},{"elapsed":384.714313392,"position":[6145999.982730295,2857717.0991808283,5.084118442785547e-6],"velocity":[-3233.541911345344,6953.001661692889,3.9034441394667434e-8]},{"elapsed":423.385385278,"position":[6015113.409241519,3123776.4082465433,6.737446957219513e-6],"velocity":[-3534.6026072660866,6804.923594569369,4.654934685858022e-8]},{"elapsed":461.833813847,"position":[5873565.709582861,3382377.239049698,8.67775979145647e-6],"velocity":[-3827.2268192736133,6644.783043702982,5.44473852834631e-8]},{"elapsed":500.337651938,"position":[5720676.581457687,3634936.7799456194,0.000010932469507062726],"velocity":[-4113.018340332627,6471.810086824766,6.272545425254571e-8]},{"elapsed":538.793374807,"position":[5557143.128239333,3880296.1519273133,0.00001350867725915659],"velocity":[-4390.665952485439,6286.7931209919025,7.130486541987151e-8]},{"elapsed":577.119812228,"position":[5383693.02773496,4117523.1681592413,0.00001640955292428148],"velocity":[-4659.1148917907885,6090.555140123361,8.011030070997805e-8]},{"elapsed":615.540594817,"position":[5199656.359711009,4347563.366040146,0.000019660338161310214],"velocity":[-4919.435229616536,5882.337913448691,8.913912275666822e-8]},{"elapsed":653.946969434,"position":[5005869.4262301Line truncated
|
||||
@@ -214,3 +214,78 @@ fn test_radiation_pressure_agrees_with_nyx() {
|
||||
"worst velocity disagreement {velocity:.3e} m/s over the day"
|
||||
);
|
||||
}
|
||||
|
||||
/// Drag through NRLMSISE-00, where the drivers match exactly and the models do not.
|
||||
///
|
||||
/// This is the one comparison here with no driver mismatch at all. The fixture's space
|
||||
/// weather is written out of the same embedded GFZ record this crate reads, into the
|
||||
/// CelesTrak-shaped CSV nyx wants, and nyx takes the 81-day centered average from a
|
||||
/// column rather than computing it, so that column carries this crate's own value.
|
||||
/// Printed side by side, both tools feed the model identical numbers: f107 170.5,
|
||||
/// f107a 154.3333333333331 and all seven ap values the same. The Orekit equivalent
|
||||
/// could not manage that and lost 0.028% to CSSI publishing its average rounded.
|
||||
///
|
||||
/// # What is left is the two implementations, and the difference is nyx's
|
||||
///
|
||||
/// With the drivers identical the densities still differ by 2.6e-4, consistently and in
|
||||
/// the same direction. Run against case 0 of the seventeen standard NRLMSISE-00 cases,
|
||||
/// which pin the NRL reference to sixteen digits: **this crate reproduces the published
|
||||
/// total mass density exactly, and nyx is 1.08e-3 away from it**. So the disagreement
|
||||
/// below is nyx departing from the reference, not this crate.
|
||||
///
|
||||
/// One element is unmatched and is ours: nyx takes geodetic altitude against the
|
||||
/// planetary constants kernel's Earth and this crate against WGS84, which puts the same
|
||||
/// spacecraft 0.400 m apart. That is worth about 7e-6 of the density, two orders below
|
||||
/// the model difference, and it could not be removed -- ANISE's transform replaces a
|
||||
/// frame's ellipsoid from the loaded kernel.
|
||||
///
|
||||
/// Bounded as a fraction of what drag does over the arc rather than in meters, since
|
||||
/// the meters depend on the arc and the fraction does not.
|
||||
#[cfg(feature = "msis")]
|
||||
#[test]
|
||||
fn test_msis_drag_agrees_with_nyx() {
|
||||
let Some(almanac) = almanac() else {
|
||||
eprintln!("skipped: kernels not fetched, run scripts/fetch-kernels.sh");
|
||||
return;
|
||||
};
|
||||
let arc = fixture("drag-msis");
|
||||
assert!(arc.setup.forces[1].contains("NRLMSISE00"));
|
||||
|
||||
let (position, velocity) = compare(
|
||||
&arc,
|
||||
&[
|
||||
TwoBody { body: Body::EARTH },
|
||||
Drag {
|
||||
density: DensityModel::Nrlmsise00,
|
||||
},
|
||||
],
|
||||
Some(Spacecraft::cannonball(500.0, 5.0, 2.2, 1.3)),
|
||||
Some(almanac.clone()),
|
||||
);
|
||||
|
||||
// What drag moved the orbit over the day, so the disagreement reads against it.
|
||||
let start = &arc.states[0];
|
||||
let epoch = Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0);
|
||||
let initial = Orbit
|
||||
.with_position(Vector3::from_column_slice(&start.position))
|
||||
.with_velocity(Vector3::from_column_slice(&start.velocity));
|
||||
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();
|
||||
|
||||
// Measured 2026-10-03: 8.55 m, 0.0286% of the drag effect, 9.66e-3 m/s. That is the
|
||||
// 2.6e-4 density difference integrated along-track, and it is nyx's departure from
|
||||
// the reference rather than this crate's.
|
||||
assert!(
|
||||
position < 1.0e-3 * effect,
|
||||
"worst disagreement {position:.3e} m is {:.4}% of the {effect:.0} m drag moves it, \
|
||||
and {velocity:.3e} m/s",
|
||||
position / effect * 100.0
|
||||
);
|
||||
}
|
||||
Reference in new issue
Block a user