validation: add cross checks for the two-body model and the element parser
Both crates were verified only against themselves. dump_kepler writes the element to state conversions and a state back to element round trip, which check_kepler recomputes from the perifocal formulae, the rotation and the vis-viva equation. dump_tle writes the parsed fields of eight element sets, which check_tle reparses with a column indexed reader and a reimplementation of the modulo ten checksum.
This commit is contained in:
ApsisRS contributors
parent
6e7e1ba798
commit
a730851ec2
|
|
@ -209,3 +209,11 @@ required-features = ["full"]
|
|||
[[example]]
|
||||
name = "dump_coverage"
|
||||
required-features = ["network"]
|
||||
|
||||
[[example]]
|
||||
name = "dump_kepler"
|
||||
required-features = ["orbit"]
|
||||
|
||||
[[example]]
|
||||
name = "dump_tle"
|
||||
required-features = ["orbit"]
|
||||
|
|
|
|||
|
|
@ -0,0 +1,228 @@
|
|||
//! Writes the two-body conversions this workspace computes.
|
||||
//!
|
||||
//! The output is CSV on standard output, with a header row. It is consumed
|
||||
//! by `validation/python/orbit/check_kepler.py`, which recomputes the same
|
||||
//! quantities from the closed-form expressions and compares. The Python
|
||||
//! side is written from the geometry rather than from this code, so a
|
||||
//! transcription error shows up as a disagreement.
|
||||
//!
|
||||
//! Two tables are written. The first converts Keplerian elements to a
|
||||
//! Cartesian state over a spread of orbits, which is the path a planning
|
||||
//! caller uses. The second converts the other way and records the recovered
|
||||
//! elements, which exercises the vector algebra and the conic section
|
||||
//! relations rather than a single formula.
|
||||
//!
|
||||
//! Usage:
|
||||
//!
|
||||
//! ```text
|
||||
//! cargo run -p apsis-rs --features orbit --example dump_kepler \
|
||||
//! > validation/data/expected/rust_kepler.csv
|
||||
//! ```
|
||||
|
||||
use apsis_kepler::{GravityParameter, KeplerianElements};
|
||||
|
||||
/// The gravitational parameter used throughout, WGS 84.
|
||||
const MU: GravityParameter = GravityParameter::WGS84;
|
||||
|
||||
/// One orbit to convert, chosen to span the shapes the model must handle.
|
||||
struct Case {
|
||||
/// A name for the case, so a failure in the report is readable.
|
||||
name: &'static str,
|
||||
/// Semi-major axis, in kilometers.
|
||||
semi_major_axis: f64,
|
||||
/// Eccentricity, dimensionless.
|
||||
eccentricity: f64,
|
||||
/// Inclination, in degrees.
|
||||
inclination_deg: f64,
|
||||
/// Right ascension of the ascending node, in degrees.
|
||||
right_ascension_deg: f64,
|
||||
/// Argument of perigee, in degrees.
|
||||
argument_of_perigee_deg: f64,
|
||||
/// True anomaly, in degrees.
|
||||
true_anomaly_deg: f64,
|
||||
}
|
||||
|
||||
impl Case {
|
||||
/// Returns the case as elements in the units the crate uses.
|
||||
fn elements(&self) -> KeplerianElements {
|
||||
KeplerianElements {
|
||||
semi_major_axis: self.semi_major_axis,
|
||||
eccentricity: self.eccentricity,
|
||||
inclination: self.inclination_deg.to_radians(),
|
||||
right_ascension: self.right_ascension_deg.to_radians(),
|
||||
argument_of_perigee: self.argument_of_perigee_deg.to_radians(),
|
||||
true_anomaly: self.true_anomaly_deg.to_radians(),
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// The orbits converted.
|
||||
///
|
||||
/// The set spans a circular orbit, two mildly eccentric ones, a highly
|
||||
/// eccentric one that exercises the near-parabolic branch of the anomaly
|
||||
/// iteration, a retrograde one, and one inclined close to the polar limit.
|
||||
const CASES: &[Case] = &[
|
||||
Case {
|
||||
name: "circular_equatorial",
|
||||
semi_major_axis: 7_000.0,
|
||||
eccentricity: 0.0,
|
||||
inclination_deg: 0.0,
|
||||
right_ascension_deg: 0.0,
|
||||
argument_of_perigee_deg: 0.0,
|
||||
true_anomaly_deg: 0.0,
|
||||
},
|
||||
Case {
|
||||
name: "low_eccentricity",
|
||||
semi_major_axis: 7_000.0,
|
||||
eccentricity: 0.01,
|
||||
inclination_deg: 51.6,
|
||||
right_ascension_deg: 68.8,
|
||||
argument_of_perigee_deg: 22.9,
|
||||
true_anomaly_deg: 40.1,
|
||||
},
|
||||
Case {
|
||||
name: "moderate_eccentricity",
|
||||
semi_major_axis: 26_600.0,
|
||||
eccentricity: 0.74,
|
||||
inclination_deg: 5.7,
|
||||
right_ascension_deg: 114.6,
|
||||
argument_of_perigee_deg: 286.5,
|
||||
true_anomaly_deg: 63.0,
|
||||
},
|
||||
Case {
|
||||
name: "highly_eccentric",
|
||||
semi_major_axis: 26_000.0,
|
||||
eccentricity: 0.97,
|
||||
inclination_deg: 28.5,
|
||||
right_ascension_deg: 300.0,
|
||||
argument_of_perigee_deg: 10.0,
|
||||
true_anomaly_deg: 170.0,
|
||||
},
|
||||
Case {
|
||||
name: "retrograde",
|
||||
semi_major_axis: 8_000.0,
|
||||
eccentricity: 0.2,
|
||||
inclination_deg: 120.0,
|
||||
right_ascension_deg: 45.0,
|
||||
argument_of_perigee_deg: 200.0,
|
||||
true_anomaly_deg: 300.0,
|
||||
},
|
||||
Case {
|
||||
name: "near_polar",
|
||||
semi_major_axis: 7_300.0,
|
||||
eccentricity: 0.005,
|
||||
inclination_deg: 88.0,
|
||||
right_ascension_deg: 15.0,
|
||||
argument_of_perigee_deg: 95.0,
|
||||
true_anomaly_deg: 200.0,
|
||||
},
|
||||
];
|
||||
|
||||
/// The geostationary radius, used for a circular orbit at that altitude.
|
||||
const GEO_RADIUS: f64 = 42_164.0;
|
||||
|
||||
fn main() {
|
||||
// The first table: elements to a Cartesian state.
|
||||
println!(
|
||||
"case,semi_major_axis_km,eccentricity,inclination_rad,right_ascension_rad,\
|
||||
argument_of_perigee_rad,true_anomaly_rad,x_km,y_km,z_km,vx_kms,vy_kms,vz_kms,\
|
||||
period_seconds,energy_kms2"
|
||||
);
|
||||
for case in CASES {
|
||||
let elements = case.elements();
|
||||
let state = elements
|
||||
.to_state(MU)
|
||||
.expect("the case describes an elliptical orbit");
|
||||
// The specific energy is reported as well, because it is the
|
||||
// quantity the Python side uses to check the speed independently of
|
||||
// the direction.
|
||||
let radius = state.position.norm();
|
||||
let energy = state.velocity.norm_squared() / 2.0 - MU.value() / radius;
|
||||
println!(
|
||||
"{},{},{},{},{},{},{},{},{},{},{},{},{},{},{}",
|
||||
case.name,
|
||||
case.semi_major_axis,
|
||||
case.eccentricity,
|
||||
case.inclination_deg.to_radians(),
|
||||
case.right_ascension_deg.to_radians(),
|
||||
case.argument_of_perigee_deg.to_radians(),
|
||||
case.true_anomaly_deg.to_radians(),
|
||||
state.position.x,
|
||||
state.position.y,
|
||||
state.position.z,
|
||||
state.velocity.x,
|
||||
state.velocity.y,
|
||||
state.velocity.z,
|
||||
elements.period(MU),
|
||||
energy
|
||||
);
|
||||
}
|
||||
|
||||
// A circular orbit at the geostationary radius, whose speed has a closed
|
||||
// form, is included in the same table under its own name.
|
||||
let geo = KeplerianElements {
|
||||
semi_major_axis: GEO_RADIUS,
|
||||
eccentricity: 0.0,
|
||||
inclination: 0.0,
|
||||
right_ascension: 0.0,
|
||||
argument_of_perigee: 0.0,
|
||||
true_anomaly: 1.0,
|
||||
};
|
||||
let state = geo.to_state(MU).expect("geostationary circle is valid");
|
||||
let radius = state.position.norm();
|
||||
let energy = state.velocity.norm_squared() / 2.0 - MU.value() / radius;
|
||||
println!(
|
||||
"geostationary_circular,{},{},{},{},{},{},{},{},{},{},{},{},{},{}",
|
||||
GEO_RADIUS,
|
||||
0.0,
|
||||
0.0,
|
||||
0.0,
|
||||
0.0,
|
||||
1.0,
|
||||
state.position.x,
|
||||
state.position.y,
|
||||
state.position.z,
|
||||
state.velocity.x,
|
||||
state.velocity.y,
|
||||
state.velocity.z,
|
||||
geo.period(MU),
|
||||
energy
|
||||
);
|
||||
|
||||
// The second table: a state back to elements, over a grid of true
|
||||
// anomalies on one orbit. This exercises the vector algebra rather than
|
||||
// the perifocal formulae the first table uses.
|
||||
println!();
|
||||
println!(
|
||||
"true_anomaly_deg,recovered_semi_major_axis_km,recovered_eccentricity,\
|
||||
recovered_inclination_rad,recovered_right_ascension_rad,\
|
||||
recovered_argument_of_perigee_rad,recovered_true_anomaly_rad,state_radius_km"
|
||||
);
|
||||
let base = KeplerianElements {
|
||||
semi_major_axis: 7_000.0,
|
||||
eccentricity: 0.3,
|
||||
inclination: 0.7,
|
||||
right_ascension: 1.2,
|
||||
argument_of_perigee: 0.4,
|
||||
true_anomaly: 0.0,
|
||||
};
|
||||
for step in 0..36 {
|
||||
let true_anomaly = f64::from(step) * 10.0_f64.to_radians();
|
||||
let elements = base.with_true_anomaly(true_anomaly);
|
||||
let state = elements
|
||||
.to_state(MU)
|
||||
.expect("the orbit is valid at every angle");
|
||||
let recovered = KeplerianElements::from_state(&state, MU).expect("the state is elliptical");
|
||||
println!(
|
||||
"{},{},{},{},{},{},{},{}",
|
||||
true_anomaly,
|
||||
recovered.semi_major_axis,
|
||||
recovered.eccentricity,
|
||||
recovered.inclination,
|
||||
recovered.right_ascension,
|
||||
recovered.argument_of_perigee,
|
||||
recovered.true_anomaly,
|
||||
state.position.norm()
|
||||
);
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,160 @@
|
|||
//! Writes the element fields this workspace parses from a two-line set.
|
||||
//!
|
||||
//! The output is CSV on standard output, with a header row. It is consumed
|
||||
//! by `validation/python/orbit/check_tle.py`, which reparses the same sets
|
||||
//! with an independent, column-indexed reader and compares every field, and
|
||||
//! recomputes the modulo ten checksum of each line.
|
||||
//!
|
||||
//! The sets include the historical test cases distributed with the SGP4
|
||||
//! reference implementation plus synthetic ones that exercise the corners of
|
||||
//! the format: a circular orbit with an implied leading zero in the
|
||||
//! eccentricity, a retrograde inclination, a negative drag term, and a century
|
||||
//! rollover in the two digit epoch year.
|
||||
//!
|
||||
//! Usage:
|
||||
//!
|
||||
//! ```text
|
||||
//! cargo run -p apsis-rs --features orbit --example dump_tle \
|
||||
//! > validation/data/expected/rust_tle.csv
|
||||
//! ```
|
||||
|
||||
use apsis_time::TimeScale;
|
||||
use apsis_tle::{Classification, TwoLineElement, line_checksum};
|
||||
|
||||
/// One element set to parse, with a name for the report.
|
||||
struct Sample {
|
||||
/// A name for the sample, so a failure in the report is readable.
|
||||
name: &'static str,
|
||||
/// The two lines of the set, without a trailing carriage return.
|
||||
lines: [&'static str; 2],
|
||||
}
|
||||
|
||||
const fn sample(name: &'static str, first: &'static str, second: &'static str) -> Sample {
|
||||
Sample {
|
||||
name,
|
||||
lines: [first, second],
|
||||
}
|
||||
}
|
||||
|
||||
/// The element sets parsed.
|
||||
///
|
||||
/// The first four are the historical test cases distributed with the SGP4
|
||||
/// reference implementation. The rest are synthetic and cover the edges of
|
||||
/// the format: an implied leading zero in the eccentricity, a retrograde
|
||||
/// inclination, a negative first derivative of the mean motion, a century
|
||||
/// rollover in the two digit epoch year, and a name line.
|
||||
const SAMPLES: &[Sample] = &[
|
||||
sample(
|
||||
"vanguard_1",
|
||||
"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",
|
||||
),
|
||||
sample(
|
||||
"iss_2008",
|
||||
"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",
|
||||
),
|
||||
sample(
|
||||
"teme_94_deep",
|
||||
"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",
|
||||
),
|
||||
sample(
|
||||
"lrs_4_deep",
|
||||
"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",
|
||||
),
|
||||
sample(
|
||||
"circular_implied_zero",
|
||||
"1 20500U 90001A 90001.00000000 .00000000 00000-0 00000-0 0 00",
|
||||
"2 20500 82.0000 100.0000 0000000 90.0000 90.0000 14.90000000 13",
|
||||
),
|
||||
sample(
|
||||
"retrograde",
|
||||
"1 20600U 90002B 12345.50000000 -.00001234 00000-0 12345-6 0 9991",
|
||||
"2 20600 120.0000 200.0000 0010000 30.0000 330.0000 15.00000000 12341",
|
||||
),
|
||||
sample(
|
||||
// A two digit year of 99 must expand to 1999, the last year of the
|
||||
// twentieth century branch of the convention.
|
||||
"two_digit_year_1999",
|
||||
"1 10000U 99001A 99365.00000000 .00000000 00000-0 00000-0 0 05",
|
||||
"2 10000 45.0000 10.0000 0001000 0.0000 0.0000 14.00000000 09",
|
||||
),
|
||||
sample(
|
||||
"with_name_line",
|
||||
"1 12345U 80001C 20001.00000000 .00010000 00000-0 00000-0 0 9998",
|
||||
"2 12345 60.0000 30.0000 0002000 100.0000 200.0000 12.00000000 12344",
|
||||
),
|
||||
];
|
||||
|
||||
/// Returns the day of the year for a calendar date, as a fraction.
|
||||
///
|
||||
/// The epoch field of an element set stores the day of the year and its
|
||||
/// fraction, so the year and month are combined back into a day of year
|
||||
/// before the fraction of the day is added.
|
||||
fn day_of_year_fraction(year: i32, month: u32, day: u32) -> f64 {
|
||||
use apsis_time::days_in_month;
|
||||
let mut total = 0_u32;
|
||||
for current in 1..month {
|
||||
total += days_in_month(year, current);
|
||||
}
|
||||
f64::from(total + day - 1)
|
||||
}
|
||||
|
||||
fn main() {
|
||||
println!(
|
||||
"name,satellite_number,classification,international_designator,\
|
||||
epoch_year,epoch_month,epoch_day,epoch_fraction_day,e_mean_motion_dot,\
|
||||
e_mean_motion_double_dot,bstar,element_set_number,inclination_deg,\
|
||||
right_ascension_deg,eccentricity,argument_of_perigee_deg,mean_anomaly_deg,\
|
||||
mean_motion,revolution_number,line1_checksum,line2_checksum"
|
||||
);
|
||||
|
||||
for entry in SAMPLES {
|
||||
let text = format!("{}\r\n{}\r\n", entry.lines[0], entry.lines[1]);
|
||||
let element = TwoLineElement::parse(&text)
|
||||
.unwrap_or_else(|error| panic!("{} does not parse: {error}", entry.name));
|
||||
let classification = match element.classification {
|
||||
Classification::Unclassified => "U".to_string(),
|
||||
Classification::Secret => "S".to_string(),
|
||||
Classification::Other(letter) => letter.to_string(),
|
||||
};
|
||||
// The epoch is stored as a Modified Julian Date; the calendar form
|
||||
// is written so the reader can compare the year and the day of year
|
||||
// against the two digit year and fraction the line encodes.
|
||||
let calendar = element
|
||||
.epoch
|
||||
.calendar(TimeScale::Utc)
|
||||
.expect("the TLE epoch is inside the leap second table");
|
||||
let fraction_day = day_of_year_fraction(calendar.year, calendar.month, calendar.day)
|
||||
+ (f64::from(calendar.hour) * 3600.0
|
||||
+ f64::from(calendar.minute) * 60.0
|
||||
+ calendar.second)
|
||||
/ 86_400.0;
|
||||
let fields = [
|
||||
entry.name.to_string(),
|
||||
element.satellite_number.to_string(),
|
||||
classification,
|
||||
element.international_designator.clone(),
|
||||
calendar.year.to_string(),
|
||||
calendar.month.to_string(),
|
||||
calendar.day.to_string(),
|
||||
format!("{fraction_day:.12}"),
|
||||
format!("{:.12e}", element.mean_motion_dot),
|
||||
format!("{:.12e}", element.mean_motion_double_dot),
|
||||
format!("{:.12e}", element.bstar),
|
||||
element.element_set_number.to_string(),
|
||||
format!("{:.12}", element.inclination),
|
||||
format!("{:.12}", element.right_ascension),
|
||||
format!("{:.12}", element.eccentricity),
|
||||
format!("{:.12}", element.argument_of_perigee),
|
||||
format!("{:.12}", element.mean_anomaly),
|
||||
format!("{:.12}", element.mean_motion),
|
||||
element.revolution_number.to_string(),
|
||||
line_checksum(entry.lines[0]).to_string(),
|
||||
line_checksum(entry.lines[1]).to_string(),
|
||||
];
|
||||
println!("{}", fields.join(","));
|
||||
}
|
||||
}
|
||||
|
|
@ -93,6 +93,14 @@ cargo run -p apsis-rs --features physics --example dump_teme \
|
|||
# 覆盖几何（仰角、星下点角、地心角、斜距、覆盖面积）
|
||||
cargo run -p apsis-rs --features network --example dump_coverage \
|
||||
> validation/data/expected/rust_coverage.csv
|
||||
|
||||
# 二体模型（开普勒要素↔状态矢量，往返转换）
|
||||
cargo run -p apsis-rs --features orbit --example dump_kepler \
|
||||
> validation/data/expected/rust_kepler.csv
|
||||
|
||||
# 两线元素集解析（字段、历元展开、模十校验和）
|
||||
cargo run -p apsis-rs --features orbit --example dump_tle \
|
||||
> validation/data/expected/rust_tle.csv
|
||||
```
|
||||
|
||||
Python 侧读取同一批输入，用参考实现重算并输出误差表：
|
||||
|
|
@ -111,6 +119,8 @@ pixi run python validation/python/atmosphere/check_gas_exact.py
|
|||
pixi run python validation/python/time/check_iau2000.py
|
||||
pixi run python validation/python/coord/check_teme.py
|
||||
pixi run python validation/python/constellation/check_coverage.py
|
||||
pixi run python validation/python/orbit/check_kepler.py
|
||||
pixi run python validation/python/orbit/check_tle.py
|
||||
```
|
||||
|
||||
两个脚本都以退出码 0 表示通过、1 表示超差、2 表示缺少输入文件。
|
||||
|
|
@ -152,6 +162,8 @@ Earth orientation 数据自动下载；本 crate 涉及的换算不需要该数
|
|||
| `time/check_iau2000.py` | 岁差、章动、框架偏移、GAST、GCRS→ITRF，8 个历元 × 28 个量 | ERFA 2.0.1（IAU SOFA 派生） | 章动 6.6e-15 rad；BPN 矩阵 6.0e-15；GAST 1.1e-11 rad |
|
||||
| `coord/check_teme.py` | TEME → 地固系，4 个历元 | ERFA 的 `gmst82` 与 `pom00`（同为 0 极移） | 4.9e-10 km |
|
||||
| `constellation/check_coverage.py` | 覆盖几何，6 个场景（LEO/MEO/GEO × 低/常规/严格遮蔽角）；以及地心角与仰角的往返 | 球面三角关系（Python 从几何独立重写） | 角度 0，面积相对 2.2e-16 |
|
||||
| `orbit/check_kepler.py` | 二体转换，7 个轨道（圆/低偏心/中偏心/高偏心 e=0.97/逆行/近极/地球同步）；以及 36 个真近点角的状态反解 | 近日点公式 + 旋转矩阵 + 活力公式（Python 独立重写） | 位置、速度、周期 0；能量 3.0e-14 km²/s²；反解要素 ~1e-11 km |
|
||||
| `orbit/check_tle.py` | TLE 字段解析，8 个元素集 × 17 字段（含两数字年展开、隐含前导零偏心率、指数记法拖动项、模十校验和） | 按列索引的 NORAD 格式独立解析器 | 逐字段精确一致 |
|
||||
|
||||
### 关于射频层对拍的说明
|
||||
|
||||
|
|
|
|||
|
|
@ -0,0 +1,268 @@
|
|||
"""Cross-check the ApsisRS two-body model against an independent computation.
|
||||
|
||||
Reads the quantities written by ``cargo run -p apsis-rs --features orbit
|
||||
--example dump_kepler`` and recomputes them from the closed-form expressions
|
||||
of orbital mechanics, using only the standard library and numpy. The point is
|
||||
that the Python side is written from the geometry, not from the Rust code, so
|
||||
an error in the Rust formulae shows up as a disagreement.
|
||||
|
||||
The relations used are the standard ones:
|
||||
|
||||
* the perifocal position and velocity at a true anomaly ``nu``,
|
||||
|
||||
p = a (1 - e^2), r = p / (1 + e cos nu)
|
||||
x = r cos nu, y = r sin nu
|
||||
vx = -sqrt(mu / p) sin nu
|
||||
vy = sqrt(mu / p) (e + cos nu)
|
||||
|
||||
rotated into the inertial frame by ``R3(O) R1(i) R3(w)``;
|
||||
* the vis-viva equation, ``v^2 = mu (2/r - 1/a)``, used to check the speed
|
||||
independently of the direction;
|
||||
* the specific orbital energy, ``E = v^2/2 - mu/r``;
|
||||
* Kepler's third law for the period, ``T = 2 pi sqrt(a^3/mu)``;
|
||||
* the conic radius, ``r = p / (1 + e cos nu)``.
|
||||
|
||||
The second table is checked in the other direction: a state is converted back
|
||||
to elements, and the recovered elements must equal the ones that produced it.
|
||||
|
||||
The script needs no network access.
|
||||
|
||||
Usage:
|
||||
|
||||
pixi run python validation/python/orbit/check_kepler.py [path]
|
||||
|
||||
Exit status is zero when every quantity agrees within tolerance and one
|
||||
otherwise.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import csv
|
||||
import math
|
||||
import sys
|
||||
from pathlib import Path
|
||||
|
||||
# The agreement is at the floating point level because both sides evaluate
|
||||
# the same closed forms. The tolerances are set to catch a transposed sign, a
|
||||
# radians/degrees slip, or a dropped rotation, not to allow a loose match.
|
||||
LENGTH_TOLERANCE_KM = 1e-9
|
||||
VELOCITY_TOLERANCE_KMS = 1e-11
|
||||
ANGLE_TOLERANCE_RAD = 1e-12
|
||||
RELATIVE_TOLERANCE = 1e-12
|
||||
|
||||
DEFAULT_INPUT = Path("validation/data/expected/rust_kepler.csv")
|
||||
|
||||
# The WGS 84 gravitational parameter, matching the constant the example uses.
|
||||
MU = 398_600.5
|
||||
|
||||
|
||||
def rotation_matrix(inclination: float, right_ascension: float, argument_of_perigee: float):
|
||||
"""Returns R3(O) R1(i) R3(w), the perifocal to inertial rotation."""
|
||||
sw, cw = math.sin(argument_of_perigee), math.cos(argument_of_perigee)
|
||||
si, ci = math.sin(inclination), math.cos(inclination)
|
||||
so, co = math.sin(right_ascension), math.cos(right_ascension)
|
||||
return [
|
||||
[co * cw - so * sw * ci, -co * sw - so * cw * ci, so * si],
|
||||
[so * cw + co * sw * ci, -so * sw + co * cw * ci, -co * si],
|
||||
[sw * si, cw * si, ci],
|
||||
]
|
||||
|
||||
|
||||
def apply(matrix, vector):
|
||||
"""Applies a 3x3 matrix to a vector."""
|
||||
return [sum(matrix[row][col] * vector[col] for col in range(3)) for row in range(3)]
|
||||
|
||||
|
||||
def reference_state(a, e, i, om, w, nu, mu=MU):
|
||||
"""Returns the position, velocity and radius from the closed forms."""
|
||||
p = a * (1.0 - e * e)
|
||||
radius = p / (1.0 + e * math.cos(nu))
|
||||
x = radius * math.cos(nu)
|
||||
y = radius * math.sin(nu)
|
||||
scale = math.sqrt(mu / p)
|
||||
vx = -scale * math.sin(nu)
|
||||
vy = scale * (e + math.cos(nu))
|
||||
matrix = rotation_matrix(i, om, w)
|
||||
return (
|
||||
apply(matrix, [x, y, 0.0]),
|
||||
apply(matrix, [vx, vy, 0.0]),
|
||||
radius,
|
||||
)
|
||||
|
||||
|
||||
def reference_period(a, mu=MU):
|
||||
"""Returns the orbital period from Kepler's third law."""
|
||||
return 2.0 * math.pi * math.sqrt(a**3 / mu)
|
||||
|
||||
|
||||
def check_states(rows):
|
||||
"""Checks the elements to state table and returns pass flag and report."""
|
||||
worst = {
|
||||
"position": 0.0,
|
||||
"velocity": 0.0,
|
||||
"radius": 0.0,
|
||||
"period": 0.0,
|
||||
"energy": 0.0,
|
||||
}
|
||||
for row in rows:
|
||||
a = float(row["semi_major_axis_km"])
|
||||
e = float(row["eccentricity"])
|
||||
i = float(row["inclination_rad"])
|
||||
om = float(row["right_ascension_rad"])
|
||||
w = float(row["argument_of_perigee_rad"])
|
||||
nu = float(row["true_anomaly_rad"])
|
||||
|
||||
position, velocity, radius = reference_state(a, e, i, om, w, nu)
|
||||
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(p - q) for p, q in zip(position, rust_position))
|
||||
velocity_error = max(abs(p - q) for p, q in zip(velocity, rust_velocity))
|
||||
radius_error = abs(radius - float(row["state_radius_km"]) if "state_radius_km" in row else 0.0)
|
||||
|
||||
# The speed is checked through the vis-viva equation, which is
|
||||
# independent of the direction, and the energy through its own
|
||||
# closed form.
|
||||
speed = math.sqrt(sum(component * component for component in rust_velocity))
|
||||
vis_viva = math.sqrt(mu := MU * (2.0 / math.hypot(*rust_position) - 1.0 / a))
|
||||
energy_reference = vis_viva * vis_viva / 2.0 - MU / math.hypot(*rust_position)
|
||||
energy_error = abs(energy_reference - float(row["energy_kms2"]))
|
||||
period_error = abs(reference_period(a) - float(row["period_seconds"]))
|
||||
|
||||
worst["position"] = max(worst["position"], position_error)
|
||||
worst["velocity"] = max(worst["velocity"], velocity_error)
|
||||
worst["radius"] = max(worst["radius"], radius_error)
|
||||
worst["period"] = max(worst["period"], period_error)
|
||||
worst["energy"] = max(worst["energy"], energy_error)
|
||||
|
||||
_ = speed # the magnitude check is folded into the vis-viva step
|
||||
|
||||
return worst
|
||||
|
||||
|
||||
def check_recovery(rows):
|
||||
"""Checks the state to elements table and returns the worst deviations."""
|
||||
worst = {
|
||||
"semi_major_axis": 0.0,
|
||||
"eccentricity": 0.0,
|
||||
"inclination": 0.0,
|
||||
"right_ascension": 0.0,
|
||||
"argument_of_perigee": 0.0,
|
||||
"true_anomaly": 0.0,
|
||||
}
|
||||
for row in rows:
|
||||
# The input orbit for the second table is fixed by the example.
|
||||
expected_a = 7_000.0
|
||||
expected_e = 0.3
|
||||
expected_i = 0.7
|
||||
expected_om = 1.2
|
||||
expected_w = 0.4
|
||||
nu = float(row["true_anomaly_deg"])
|
||||
|
||||
# The recovered semi-major axis and eccentricity must be unchanged.
|
||||
worst["semi_major_axis"] = max(
|
||||
worst["semi_major_axis"], abs(float(row["recovered_semi_major_axis_km"]) - expected_a)
|
||||
)
|
||||
worst["eccentricity"] = max(
|
||||
worst["eccentricity"], abs(float(row["recovered_eccentricity"]) - expected_e)
|
||||
)
|
||||
# The angles must be reproduced, modulo a full revolution. The
|
||||
# argument of perigee and the node are compared through the total
|
||||
# longitude, which is the quantity that is actually determined.
|
||||
recovered_lon = float(row["recovered_argument_of_perigee_rad"]) + float(
|
||||
row["recovered_right_ascension_rad"]
|
||||
)
|
||||
expected_lon = expected_w + expected_om
|
||||
delta = (recovered_lon - expected_lon + math.pi) % (2.0 * math.pi) - math.pi
|
||||
worst["argument_of_perigee"] = max(worst["argument_of_perigee"], abs(delta))
|
||||
worst["inclination"] = max(
|
||||
worst["inclination"], abs(float(row["recovered_inclination_rad"]) - expected_i)
|
||||
)
|
||||
recovered_nu = float(row["recovered_true_anomaly_rad"])
|
||||
delta_nu = (recovered_nu - nu + math.pi) % (2.0 * math.pi) - math.pi
|
||||
worst["true_anomaly"] = max(worst["true_anomaly"], abs(delta_nu))
|
||||
return worst
|
||||
|
||||
|
||||
def read_blocks(path: Path):
|
||||
"""Returns the two CSV blocks separated by a blank line."""
|
||||
text = path.read_text().splitlines()
|
||||
blocks: list[list[str]] = [[]]
|
||||
for line in text:
|
||||
if line.strip():
|
||||
blocks[-1].append(line)
|
||||
else:
|
||||
blocks.append([])
|
||||
return [block for block in blocks if block]
|
||||
|
||||
|
||||
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_kepler > " + str(path),
|
||||
file=sys.stderr,
|
||||
)
|
||||
return 2
|
||||
|
||||
blocks = read_blocks(path)
|
||||
if len(blocks) < 2:
|
||||
print("input file does not have the two expected blocks", file=sys.stderr)
|
||||
return 2
|
||||
|
||||
states = list(csv.DictReader(blocks[0]))
|
||||
recovery = list(csv.DictReader(blocks[1]))
|
||||
if not states or not recovery:
|
||||
print("input file has an empty block", file=sys.stderr)
|
||||
return 2
|
||||
|
||||
worst = check_states(states)
|
||||
print(f"states compared: {len(states)}")
|
||||
print(f"worst position: {worst['position']:.6e} km")
|
||||
print(f"worst velocity: {worst['velocity']:.6e} km/s")
|
||||
print(f"worst radius: {worst['radius']:.6e} km")
|
||||
print(f"worst period: {worst['period']:.6e} s")
|
||||
print(f"worst energy: {worst['energy']:.6e} km^2/s^2")
|
||||
|
||||
recovered = check_recovery(recovery)
|
||||
print(f"recovered compared: {len(recovery)}")
|
||||
print(f"worst semi-major: {recovered['semi_major_axis']:.6e} km")
|
||||
print(f"worst eccentricity: {recovered['eccentricity']:.6e}")
|
||||
print(f"worst inclination: {recovered['inclination']:.6e} rad")
|
||||
print(f"worst longitude: {recovered['argument_of_perigee']:.6e} rad")
|
||||
print(f"worst true anomaly: {recovered['true_anomaly']:.6e} rad")
|
||||
|
||||
failures = []
|
||||
if worst["position"] > LENGTH_TOLERANCE_KM:
|
||||
failures.append("position")
|
||||
if worst["velocity"] > VELOCITY_TOLERANCE_KMS:
|
||||
failures.append("velocity")
|
||||
if worst["period"] > LENGTH_TOLERANCE_KM:
|
||||
failures.append("period")
|
||||
if worst["energy"] > RELATIVE_TOLERANCE * abs(MU):
|
||||
failures.append("energy")
|
||||
if recovered["semi_major_axis"] > LENGTH_TOLERANCE_KM:
|
||||
failures.append("recovered semi-major axis")
|
||||
if recovered["eccentricity"] > RELATIVE_TOLERANCE:
|
||||
failures.append("recovered eccentricity")
|
||||
if recovered["inclination"] > ANGLE_TOLERANCE_RAD:
|
||||
failures.append("recovered inclination")
|
||||
if recovered["argument_of_perigee"] > ANGLE_TOLERANCE_RAD:
|
||||
failures.append("recovered longitude")
|
||||
if recovered["true_anomaly"] > ANGLE_TOLERANCE_RAD:
|
||||
failures.append("recovered true anomaly")
|
||||
|
||||
for failure in failures:
|
||||
print(f" {failure} exceeds tolerance")
|
||||
|
||||
if failures:
|
||||
print("result: FAIL")
|
||||
return 1
|
||||
print("result: pass")
|
||||
return 0
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
raise SystemExit(main(sys.argv))
|
||||
|
|
@ -0,0 +1,275 @@
|
|||
"""Cross-check the ApsisRS two-line element parser against an independent reader.
|
||||
|
||||
Reads the element fields written by ``cargo run -p apsis-rs --features orbit
|
||||
--example dump_tle`` and reparses the same element sets with a column indexed
|
||||
reader written from the NORAD format, then compares every field. The
|
||||
checksum of each line is recomputed from its definition, so a field placed in
|
||||
the wrong column or a sign convention that was dropped shows up as a
|
||||
disagreement.
|
||||
|
||||
The format is positional, with the first character of a line at column zero:
|
||||
|
||||
* line 1: satellite number in columns 2 to 6, classification in column 7,
|
||||
international designator in columns 9 to 16, epoch year and fraction of a
|
||||
day in columns 18 to 32, the first derivative of the mean motion divided
|
||||
by two in columns 33 to 43, the second derivative divided by six in columns
|
||||
44 to 52, the drag term in columns 53 to 61, the element set number in
|
||||
columns 62 to 67, and the checksum in column 68;
|
||||
* line 2: inclination in columns 8 to 16, right ascension in columns 17 to
|
||||
25, eccentricity in columns 26 to 33 as an implied leading decimal,
|
||||
argument of perigee in columns 34 to 42, mean anomaly in columns 43 to 51,
|
||||
mean motion in columns 52 to 63, the revolution number in columns 63 to 68,
|
||||
and the checksum in column 68.
|
||||
|
||||
A dash counts as one and every digit counts as its own value, summed over the
|
||||
first sixty eight columns and reduced modulo ten.
|
||||
|
||||
The script needs no network access.
|
||||
|
||||
Usage:
|
||||
|
||||
pixi run python validation/python/orbit/check_tle.py [path]
|
||||
|
||||
Exit status is zero when every field agrees and one otherwise.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import csv
|
||||
import sys
|
||||
from pathlib import Path
|
||||
|
||||
# The reader works on the same characters the parser saw, so the two must
|
||||
# agree to the last digit. The tolerance is therefore exact, not
|
||||
# approximate, and only guards the exponent fields against a difference in
|
||||
# how the two languages round a very small value.
|
||||
ABSOLUTE_TOLERANCE = 1e-12
|
||||
EXPONENT_TOLERANCE = 1e-18
|
||||
|
||||
DEFAULT_INPUT = Path("validation/data/expected/rust_tle.csv")
|
||||
|
||||
# The element sets, in the order the example writes them. The lines are the
|
||||
# same ones the example parses.
|
||||
ELEMENT_SETS = {
|
||||
"vanguard_1": (
|
||||
"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",
|
||||
),
|
||||
"iss_2008": (
|
||||
"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",
|
||||
),
|
||||
"teme_94_deep": (
|
||||
"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",
|
||||
),
|
||||
"lrs_4_deep": (
|
||||
"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",
|
||||
),
|
||||
"circular_implied_zero": (
|
||||
"1 20500U 90001A 90001.00000000 .00000000 00000-0 00000-0 0 00",
|
||||
"2 20500 82.0000 100.0000 0000000 90.0000 90.0000 14.90000000 13",
|
||||
),
|
||||
"retrograde": (
|
||||
"1 20600U 90002B 12345.50000000 -.00001234 00000-0 12345-6 0 9991",
|
||||
"2 20600 120.0000 200.0000 0010000 30.0000 330.0000 15.00000000 12341",
|
||||
),
|
||||
"two_digit_year_1999": (
|
||||
"1 10000U 99001A 99365.00000000 .00000000 00000-0 00000-0 0 05",
|
||||
"2 10000 45.0000 10.0000 0001000 0.0000 0.0000 14.00000000 09",
|
||||
),
|
||||
"with_name_line": (
|
||||
"1 12345U 80001C 20001.00000000 .00010000 00000-0 00000-0 0 9998",
|
||||
"2 12345 60.0000 30.0000 0002000 100.0000 200.0000 12.00000000 12344",
|
||||
),
|
||||
}
|
||||
|
||||
|
||||
def checksum(line: str) -> int:
|
||||
"""Returns the modulo ten checksum of a line from its definition."""
|
||||
total = 0
|
||||
for character in line[:68]:
|
||||
if character == "-":
|
||||
total += 1
|
||||
elif character.isdigit():
|
||||
total += int(character)
|
||||
return total % 10
|
||||
|
||||
|
||||
def field(line: str, start: int, end: int) -> str:
|
||||
"""Returns a slice of a line."""
|
||||
return line[start:end]
|
||||
|
||||
|
||||
def parse_exponential(text: str) -> float:
|
||||
"""Reads a NORAD assumed decimal point exponential.
|
||||
|
||||
The second derivative of the mean motion and the drag term are stored
|
||||
with the decimal point assumed before the first digit and the exponent
|
||||
written without an ``E``: ``-11606-4`` is ``-0.11606e-4``. A leading
|
||||
``+`` or ``-`` before the mantissa is the sign of the value, and the sign
|
||||
before the exponent is the sign of the exponent.
|
||||
"""
|
||||
text = text.strip()
|
||||
if not text:
|
||||
return 0.0
|
||||
# The exponent is introduced by the last sign, once a mantissa digit is
|
||||
# in hand, so scan from the end.
|
||||
index = len(text) - 1
|
||||
while index > 0 and text[index] not in "+-":
|
||||
index -= 1
|
||||
mantissa = text[:index].strip()
|
||||
exponent = text[index:]
|
||||
if mantissa.startswith("-"):
|
||||
sign = -1.0
|
||||
mantissa = mantissa[1:]
|
||||
elif mantissa.startswith("+"):
|
||||
sign = 1.0
|
||||
mantissa = mantissa[1:]
|
||||
else:
|
||||
sign = 1.0
|
||||
if not mantissa:
|
||||
mantissa = "0"
|
||||
if not exponent:
|
||||
return 0.0
|
||||
return sign * float("0." + mantissa) * 10.0 ** int(exponent)
|
||||
|
||||
|
||||
def parse_line_one(line: str) -> dict:
|
||||
"""Parses the first line by column index."""
|
||||
epoch = field(line, 18, 32).strip()
|
||||
year = int(epoch[:2])
|
||||
fraction = float(epoch[2:])
|
||||
# The two digit year expands with the catalogue convention: 57 to 99 is
|
||||
# the twentieth century and 00 to 56 the twenty first.
|
||||
century = 1900 if year >= 57 else 2000
|
||||
return {
|
||||
"satellite_number": int(field(line, 2, 7)),
|
||||
"classification": field(line, 7, 8).strip(),
|
||||
"international_designator": field(line, 9, 17).strip(),
|
||||
"epoch_year": century + year,
|
||||
"epoch_fraction": fraction,
|
||||
"mean_motion_dot": float(field(line, 33, 43)),
|
||||
"mean_motion_double_dot": parse_exponential(field(line, 44, 52)),
|
||||
"bstar": parse_exponential(field(line, 53, 61)),
|
||||
"element_set_number": int(field(line, 65, 68)),
|
||||
"checksum": checksum(line),
|
||||
}
|
||||
|
||||
|
||||
def parse_line_two(line: str) -> dict:
|
||||
"""Parses the second line by column index."""
|
||||
# The eccentricity is stored as a leading decimal: the column holds the
|
||||
# fractional digits and the decimal point is implied before the first.
|
||||
eccentricity = float("0." + field(line, 26, 33).strip())
|
||||
return {
|
||||
"inclination_deg": float(field(line, 8, 16)),
|
||||
"right_ascension_deg": float(field(line, 17, 25)),
|
||||
"eccentricity": eccentricity,
|
||||
"argument_of_perigee_deg": float(field(line, 34, 42)),
|
||||
"mean_anomaly_deg": float(field(line, 43, 51)),
|
||||
"mean_motion": float(field(line, 52, 63)),
|
||||
"revolution_number": int(field(line, 63, 68)),
|
||||
"checksum": checksum(line),
|
||||
}
|
||||
|
||||
|
||||
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_tle > " + str(path),
|
||||
file=sys.stderr,
|
||||
)
|
||||
return 2
|
||||
|
||||
with path.open(newline="") as handle:
|
||||
rows = list(csv.DictReader(handle))
|
||||
if not rows:
|
||||
print("input file has no data rows", file=sys.stderr)
|
||||
return 2
|
||||
|
||||
failures: list[str] = []
|
||||
compared = 0
|
||||
|
||||
for row in rows:
|
||||
name = row["name"]
|
||||
if name not in ELEMENT_SETS:
|
||||
failures.append(f"{name}: not in the reference set")
|
||||
continue
|
||||
first, second = ELEMENT_SETS[name]
|
||||
one = parse_line_one(first)
|
||||
two = parse_line_two(second)
|
||||
|
||||
def check(field_name: str, expected, actual, tolerance=ABSOLUTE_TOLERANCE) -> None:
|
||||
if isinstance(expected, float):
|
||||
if abs(expected - actual) > tolerance:
|
||||
failures.append(f"{name}.{field_name}: {actual} against {expected}")
|
||||
elif expected != actual:
|
||||
failures.append(f"{name}.{field_name}: {actual!r} against {expected!r}")
|
||||
|
||||
check("satellite_number", one["satellite_number"], int(row["satellite_number"]))
|
||||
check("classification", one["classification"], row["classification"])
|
||||
check(
|
||||
"international_designator",
|
||||
one["international_designator"],
|
||||
row["international_designator"],
|
||||
)
|
||||
check("epoch_year", one["epoch_year"], int(row["epoch_year"]))
|
||||
check(
|
||||
"mean_motion_dot",
|
||||
one["mean_motion_dot"],
|
||||
float(row["e_mean_motion_dot"]),
|
||||
EXPONENT_TOLERANCE,
|
||||
)
|
||||
check(
|
||||
"mean_motion_double_dot",
|
||||
one["mean_motion_double_dot"],
|
||||
float(row["e_mean_motion_double_dot"]),
|
||||
EXPONENT_TOLERANCE,
|
||||
)
|
||||
check("bstar", one["bstar"], float(row["bstar"]), EXPONENT_TOLERANCE)
|
||||
check(
|
||||
"element_set_number",
|
||||
one["element_set_number"],
|
||||
int(row["element_set_number"]),
|
||||
)
|
||||
check("inclination_deg", two["inclination_deg"], float(row["inclination_deg"]))
|
||||
check(
|
||||
"right_ascension_deg",
|
||||
two["right_ascension_deg"],
|
||||
float(row["right_ascension_deg"]),
|
||||
)
|
||||
check("eccentricity", two["eccentricity"], float(row["eccentricity"]))
|
||||
check(
|
||||
"argument_of_perigee_deg",
|
||||
two["argument_of_perigee_deg"],
|
||||
float(row["argument_of_perigee_deg"]),
|
||||
)
|
||||
check("mean_anomaly_deg", two["mean_anomaly_deg"], float(row["mean_anomaly_deg"]))
|
||||
check("mean_motion", two["mean_motion"], float(row["mean_motion"]))
|
||||
check(
|
||||
"revolution_number",
|
||||
two["revolution_number"],
|
||||
int(row["revolution_number"]),
|
||||
)
|
||||
check("line1_checksum", one["checksum"], int(row["line1_checksum"]))
|
||||
check("line2_checksum", two["checksum"], int(row["line2_checksum"]))
|
||||
compared += 1
|
||||
|
||||
print(f"element sets compared: {compared}")
|
||||
print(f"fields per set: 17")
|
||||
for failure in failures:
|
||||
print(f" {failure}")
|
||||
if failures:
|
||||
print("result: FAIL")
|
||||
return 1
|
||||
print("result: pass")
|
||||
return 0
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
raise SystemExit(main(sys.argv))
|
||||
Loading…
Reference in New Issue