diff --git a/crates/apsis-rs/Cargo.toml b/crates/apsis-rs/Cargo.toml index 52556bc..3229e4a 100644 --- a/crates/apsis-rs/Cargo.toml +++ b/crates/apsis-rs/Cargo.toml @@ -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"] diff --git a/crates/apsis-rs/examples/dump_kepler.rs b/crates/apsis-rs/examples/dump_kepler.rs new file mode 100644 index 0000000..0699ae7 --- /dev/null +++ b/crates/apsis-rs/examples/dump_kepler.rs @@ -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() + ); + } +} diff --git a/crates/apsis-rs/examples/dump_tle.rs b/crates/apsis-rs/examples/dump_tle.rs new file mode 100644 index 0000000..289e6db --- /dev/null +++ b/crates/apsis-rs/examples/dump_tle.rs @@ -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(",")); + } +} diff --git a/validation/README.md b/validation/README.md index d29e0be..3f397a7 100644 --- a/validation/README.md +++ b/validation/README.md @@ -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 格式独立解析器 | 逐字段精确一致 | ### 关于射频层对拍的说明 diff --git a/validation/python/orbit/check_kepler.py b/validation/python/orbit/check_kepler.py new file mode 100644 index 0000000..fb929c3 --- /dev/null +++ b/validation/python/orbit/check_kepler.py @@ -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)) diff --git a/validation/python/orbit/check_tle.py b/validation/python/orbit/check_tle.py new file mode 100644 index 0000000..8de85e8 --- /dev/null +++ b/validation/python/orbit/check_tle.py @@ -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))