workspace: add SGP4 cross-check against the reference library
This commit is contained in:
ApsisRS contributors
parent
6eaaa86e30
commit
74594e6a04
|
|
@ -22,6 +22,10 @@ required-features = ["physics"]
|
|||
name = "dump_time"
|
||||
required-features = ["physics"]
|
||||
|
||||
[[example]]
|
||||
name = "dump_orbit"
|
||||
required-features = ["orbit"]
|
||||
|
||||
[lints]
|
||||
workspace = true
|
||||
|
||||
|
|
|
|||
|
|
@ -0,0 +1,82 @@
|
|||
//! Writes propagated states for a set of element sets and times.
|
||||
//!
|
||||
//! The output is CSV on standard output, with a header row that names the
|
||||
//! element set index, the minutes from epoch, and the position and velocity
|
||||
//! components. It is consumed by `validation/python/orbit/check_sgp4.py`,
|
||||
//! which propagates the same element sets with the reference implementation
|
||||
//! and compares.
|
||||
//!
|
||||
//! Usage:
|
||||
//!
|
||||
//! ```text
|
||||
//! cargo run -p apsis-rs --features orbit --example dump_orbit \
|
||||
//! > validation/data/expected/rust_sgp4.csv
|
||||
//! ```
|
||||
|
||||
use apsis_sgp4::{GravityModel, Satellite};
|
||||
use apsis_tle::TwoLineElement;
|
||||
|
||||
/// The element sets used by the cross-check. They are the historical test
|
||||
/// cases distributed with the SGP4 reference implementation, which cover
|
||||
/// low, medium, and high eccentricity and the deep space branch.
|
||||
const ELEMENT_SETS: &[&str] = &[
|
||||
// International Space Station, near circular low Earth orbit.
|
||||
concat!(
|
||||
"1 25544U 98067A 08264.51782528 -.00002182 00000-0 -11606-4 0 2927\r\n",
|
||||
"2 25544 51.6416 247.4627 0006703 130.5360 325.0288 15.72125391563537\r\n",
|
||||
),
|
||||
// A highly eccentric orbit, which exercises the eccentric anomaly
|
||||
// iteration and the drag terms away from circular motion.
|
||||
concat!(
|
||||
"1 00005U 58002B 00179.78495062 .00000023 00000-0 28098-4 0 4753\r\n",
|
||||
"2 00005 34.2682 348.7242 1859667 331.7664 19.3264 10.82419157413667\r\n",
|
||||
),
|
||||
// A geostationary transfer orbit with a low perigee.
|
||||
concat!(
|
||||
"1 27651U 03004A 06185.99374200 .00000109 00000-0 19909-4 0 445\r\n",
|
||||
"2 27651 7.1743 182.9443 7867251 269.0386 18.6601 1.03461646 9884\r\n",
|
||||
),
|
||||
// A deep space orbit, whose period exceeds 225 minutes.
|
||||
concat!(
|
||||
"1 09880U 77021A 06174.85710247 .00000090 00000-0 10000-3 0 7322\r\n",
|
||||
"2 09880 64.1526 313.1296 6546925 103.4683 305.3980 2.00754276 13047\r\n",
|
||||
),
|
||||
];
|
||||
|
||||
/// The times to sample, in minutes from the element epoch.
|
||||
const TIMES_MINUTES: &[f64] = &[
|
||||
0.0, 1.0, 60.0, 360.0, 720.0, 1_440.0, 2_880.0, 10_080.0, -1_440.0,
|
||||
];
|
||||
|
||||
fn main() {
|
||||
println!("set,minutes,error,x_km,y_km,z_km,vx_kms,vy_kms,vz_kms");
|
||||
|
||||
for (index, text) in ELEMENT_SETS.iter().enumerate() {
|
||||
let Ok(element) = TwoLineElement::parse(text) else {
|
||||
eprintln!("element set {index} did not parse");
|
||||
continue;
|
||||
};
|
||||
let Some(satellite) = Satellite::from_element(&element, GravityModel::WGS72) else {
|
||||
eprintln!("element set {index} was rejected");
|
||||
continue;
|
||||
};
|
||||
for minutes in TIMES_MINUTES {
|
||||
match satellite.propagate_minutes(*minutes) {
|
||||
Ok(state) => {
|
||||
println!(
|
||||
"{index},{minutes:.1},0,{:.9},{:.9},{:.9},{:.12},{:.12},{:.12}",
|
||||
state.position.x,
|
||||
state.position.y,
|
||||
state.position.z,
|
||||
state.velocity.x,
|
||||
state.velocity.y,
|
||||
state.velocity.z
|
||||
);
|
||||
}
|
||||
Err(error) => {
|
||||
println!("{index},{minutes:.1},{error},0,0,0,0,0,0");
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,167 @@
|
|||
"""Cross-check the ApsisRS SGP4 propagator against the reference library.
|
||||
|
||||
Reads the states written by ``cargo run -p apsis-rs --features orbit
|
||||
--example dump_orbit``, propagates the same element sets with the ``sgp4``
|
||||
reference implementation, and reports the largest position and velocity
|
||||
difference for each element set.
|
||||
|
||||
The script needs no network access.
|
||||
|
||||
Usage:
|
||||
|
||||
pixi run python validation/python/orbit/check_sgp4.py [path]
|
||||
|
||||
Exit status is zero when every state agrees within tolerance and one
|
||||
otherwise.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import csv
|
||||
import sys
|
||||
from pathlib import Path
|
||||
|
||||
from sgp4.earth_gravity import wgs72
|
||||
from sgp4.model import Satrec, WGS72
|
||||
|
||||
# Tolerances for the cross-check.
|
||||
#
|
||||
# These are the limits the workspace commits to for propagated position and
|
||||
# velocity. They are far above the observed agreement, which is at the
|
||||
# floating point level, and far below the accuracy of the SGP4 theory
|
||||
# itself, which is on the order of a kilometer per day of propagation.
|
||||
POSITION_TOLERANCE_KM = 1e-3
|
||||
VELOCITY_TOLERANCE_KMS = 1e-6
|
||||
|
||||
DEFAULT_INPUT = Path("validation/data/expected/rust_sgp4.csv")
|
||||
|
||||
# The element sets, in the same order as the Rust example. They are the
|
||||
# historical test cases distributed with the reference implementation.
|
||||
ELEMENT_SETS = [
|
||||
(
|
||||
"1 25544U 98067A 08264.51782528 -.00002182 00000-0 -11606-4 0 2927",
|
||||
"2 25544 51.6416 247.4627 0006703 130.5360 325.0288 15.72125391563537",
|
||||
),
|
||||
(
|
||||
"1 00005U 58002B 00179.78495062 .00000023 00000-0 28098-4 0 4753",
|
||||
"2 00005 34.2682 348.7242 1859667 331.7664 19.3264 10.82419157413667",
|
||||
),
|
||||
(
|
||||
"1 27651U 03004A 06185.99374200 .00000109 00000-0 19909-4 0 445",
|
||||
"2 27651 7.1743 182.9443 7867251 269.0386 18.6601 1.03461646 9884",
|
||||
),
|
||||
(
|
||||
"1 09880U 77021A 06174.85710247 .00000090 00000-0 10000-3 0 7322",
|
||||
"2 09880 64.1526 313.1296 6546925 103.4683 305.3980 2.00754276 13047",
|
||||
),
|
||||
]
|
||||
|
||||
|
||||
def read_rows(path: Path) -> list[dict[str, str]]:
|
||||
"""Reads the Rust output rows."""
|
||||
with path.open(newline="") as handle:
|
||||
reader = csv.DictReader(handle)
|
||||
return [row for row in reader if row.get("set")]
|
||||
|
||||
|
||||
def reference_state(satellite: Satrec, minutes: float):
|
||||
"""Returns the reference position and velocity, or None on error."""
|
||||
error, position, velocity = satellite.sgp4_tsince(minutes)
|
||||
if error != 0:
|
||||
return None
|
||||
return position, velocity
|
||||
|
||||
|
||||
def main(argv: list[str]) -> int:
|
||||
path = Path(argv[1]) if len(argv) > 1 else DEFAULT_INPUT
|
||||
if not path.exists():
|
||||
print(f"missing input file: {path}", file=sys.stderr)
|
||||
print(
|
||||
"generate it with: cargo run -p apsis-rs --features orbit "
|
||||
"--example dump_orbit > " + str(path),
|
||||
file=sys.stderr,
|
||||
)
|
||||
return 2
|
||||
|
||||
rows = read_rows(path)
|
||||
if not rows:
|
||||
print("input file has no data rows", file=sys.stderr)
|
||||
return 2
|
||||
|
||||
satellites = [
|
||||
Satrec.twoline2rv(first, second, WGS72) for first, second in ELEMENT_SETS
|
||||
]
|
||||
|
||||
worst_position = 0.0
|
||||
worst_position_at = None
|
||||
worst_velocity = 0.0
|
||||
worst_velocity_at = None
|
||||
compared = 0
|
||||
rust_errors = 0
|
||||
reference_errors = 0
|
||||
mismatched_errors = 0
|
||||
|
||||
for row in rows:
|
||||
index = int(row["set"])
|
||||
minutes = float(row["minutes"])
|
||||
rust_error = row["error"] != "0"
|
||||
reference = reference_state(satellites[index], minutes)
|
||||
|
||||
if rust_error and reference is None:
|
||||
# Both sides reject the instant, which is agreement.
|
||||
continue
|
||||
if rust_error or reference is None:
|
||||
mismatched_errors += 1
|
||||
if rust_error:
|
||||
rust_errors += 1
|
||||
else:
|
||||
reference_errors += 1
|
||||
continue
|
||||
|
||||
position, velocity = reference
|
||||
rust_position = (float(row["x_km"]), float(row["y_km"]), float(row["z_km"]))
|
||||
rust_velocity = (
|
||||
float(row["vx_kms"]),
|
||||
float(row["vy_kms"]),
|
||||
float(row["vz_kms"]),
|
||||
)
|
||||
position_error = max(abs(a - b) for a, b in zip(rust_position, position))
|
||||
velocity_error = max(abs(a - b) for a, b in zip(rust_velocity, velocity))
|
||||
|
||||
if position_error > worst_position:
|
||||
worst_position = position_error
|
||||
worst_position_at = (index, minutes)
|
||||
if velocity_error > worst_velocity:
|
||||
worst_velocity = velocity_error
|
||||
worst_velocity_at = (index, minutes)
|
||||
compared += 1
|
||||
|
||||
print(f"states compared: {compared}")
|
||||
print(f"rust only errors: {rust_errors}")
|
||||
print(f"reference errors: {reference_errors}")
|
||||
print(f"mismatched errors: {mismatched_errors}")
|
||||
print(
|
||||
f"worst position: {worst_position:.6e} km at (set, minutes) {worst_position_at}"
|
||||
)
|
||||
print(
|
||||
f"worst velocity: {worst_velocity:.6e} km/s at (set, minutes) {worst_velocity_at}"
|
||||
)
|
||||
print(f"position tolerance: {POSITION_TOLERANCE_KM:.6e} km")
|
||||
print(f"velocity tolerance: {VELOCITY_TOLERANCE_KMS:.6e} km/s")
|
||||
|
||||
if mismatched_errors > 0:
|
||||
print("result: FAIL (error handling disagrees)")
|
||||
return 1
|
||||
if (
|
||||
compared > 0
|
||||
and worst_position <= POSITION_TOLERANCE_KM
|
||||
and worst_velocity <= VELOCITY_TOLERANCE_KMS
|
||||
):
|
||||
print("result: pass")
|
||||
return 0
|
||||
print("result: FAIL")
|
||||
return 1
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
raise SystemExit(main(sys.argv))
|
||||
Loading…
Reference in New Issue