From c875aeefa6fc2f7ee056258f2a39524945e3bb32 Mon Sep 17 00:00:00 2001 From: ApsisRS contributors Date: Sat, 26 Sep 2026 23:26:13 +0800 Subject: [PATCH] validation: add the coverage geometry cross-check The network layer was the only layer without a Python cross-check, and validation/python/constellation/ was an empty placeholder that the README described as existing. Add the dump example and the reference script, which recomputes the geometry from the spherical triangle relations rather than from the Rust code, and remove the placeholder. --- crates/apsis-rs/Cargo.toml | 4 + crates/apsis-rs/examples/dump_coverage.rs | 153 ++++++++++++ validation/README.md | 7 +- validation/python/constellation/.gitkeep | 0 .../python/constellation/check_coverage.py | 227 ++++++++++++++++++ 5 files changed, 390 insertions(+), 1 deletion(-) create mode 100644 crates/apsis-rs/examples/dump_coverage.rs delete mode 100644 validation/python/constellation/.gitkeep create mode 100644 validation/python/constellation/check_coverage.py diff --git a/crates/apsis-rs/Cargo.toml b/crates/apsis-rs/Cargo.toml index 3c57100..52556bc 100644 --- a/crates/apsis-rs/Cargo.toml +++ b/crates/apsis-rs/Cargo.toml @@ -205,3 +205,7 @@ required-features = ["sim"] [[example]] name = "bench" required-features = ["full"] + +[[example]] +name = "dump_coverage" +required-features = ["network"] diff --git a/crates/apsis-rs/examples/dump_coverage.rs b/crates/apsis-rs/examples/dump_coverage.rs new file mode 100644 index 0000000..c3e7517 --- /dev/null +++ b/crates/apsis-rs/examples/dump_coverage.rs @@ -0,0 +1,153 @@ +//! Writes the coverage geometry this workspace computes for a set of cases. +//! +//! The output is CSV on standard output, with a header row. It is consumed +//! by `validation/python/constellation/check_coverage.py`, which recomputes +//! the same quantities from the spherical triangle relations and compares. +//! +//! The cases are chosen to span the interesting geometry: a low Earth orbit +//! at a low mask, where the footprint is small relative to the Earth; a +//! higher orbit, where it is large; the equatorial and polar cases; and a +//! zero mask, where the footprint reaches the horizon. +//! +//! Usage: +//! +//! ```text +//! cargo run -p apsis-rs --features network --example dump_coverage \ +//! > validation/data/expected/rust_coverage.csv +//! ``` + +use apsis_coverage::{ + earth_central_angle, elevation_from_central_angle, footprint_area, footprint_central_angle, + footprint_fraction, nadir_angle, slant_range, +}; +use core::f64::consts::PI; + +/// The WGS 84 equatorial radius, in kilometers. +/// +/// The coverage geometry is spherical, so a single radius is used. The +/// equatorial value is the one the orbit layer works in. +const BODY_RADIUS_KM: f64 = 6_378.137; + +/// One coverage case: an orbital radius and an elevation mask. +struct Case { + /// A name for the case, so a failure in the report is readable. + name: &'static str, + /// The geocentric radius of the satellite, in kilometers. + satellite_radius_km: f64, + /// The elevation mask at the ground point, in degrees. + elevation_mask_degrees: f64, +} + +/// The cases checked. +/// +/// The satellite radii span a low Earth orbit, two higher altitudes, and a +/// geostationary one. The masks span the horizon, a typical station mask, and +/// a strict one. +const CASES: &[Case] = &[ + Case { + name: "leo_low_mask", + satellite_radius_km: 6_928.137, + elevation_mask_degrees: 5.0, + }, + Case { + name: "leo_ten_degrees", + satellite_radius_km: 6_928.137, + elevation_mask_degrees: 10.0, + }, + Case { + name: "leo_strict", + satellite_radius_km: 6_928.137, + elevation_mask_degrees: 30.0, + }, + Case { + name: "meo", + satellite_radius_km: 26_560.0, + elevation_mask_degrees: 10.0, + }, + Case { + name: "geo", + satellite_radius_km: 42_164.0, + elevation_mask_degrees: 10.0, + }, + Case { + name: "horizon", + satellite_radius_km: 6_928.137, + elevation_mask_degrees: 0.0, + }, +]; + +/// A set of central angles, in degrees, for the inverse relation check. +/// +/// Each is fed to `elevation_from_central_angle` and the elevation it returns +/// is fed back through `earth_central_angle`, so the pair must compose to the +/// identity. The angles span zero to a little past the largest footprint the +/// cases produce. +const CENTRAL_ANGLES_DEGREES: &[f64] = &[0.0, 1.0, 5.0, 15.0, 30.0, 45.0, 60.0, 75.0]; + +fn main() { + // The header names every column so that the Python side can read by name + // rather than by position. + println!( + "case,satellite_radius_km,body_radius_km,elevation_mask_rad,\ + nadir_rad,central_rad,slant_range_km,footprint_central_rad,\ + footprint_area_km2,footprint_fraction" + ); + + for case in CASES { + let elevation = case.elevation_mask_degrees * PI / 180.0; + let nadir = nadir_angle(case.satellite_radius_km, BODY_RADIUS_KM, elevation) + .expect("the case geometry is usable"); + let central = earth_central_angle(case.satellite_radius_km, BODY_RADIUS_KM, elevation) + .expect("the case geometry is usable"); + let range = slant_range(case.satellite_radius_km, BODY_RADIUS_KM, elevation) + .expect("the case geometry is usable"); + let footprint_central = + footprint_central_angle(case.satellite_radius_km, BODY_RADIUS_KM, elevation) + .expect("the case geometry is usable"); + let area = footprint_area(case.satellite_radius_km, BODY_RADIUS_KM, elevation) + .expect("the case geometry is usable"); + let fraction = footprint_fraction(case.satellite_radius_km, BODY_RADIUS_KM, elevation) + .expect("the case geometry is usable"); + + println!( + "{},{},{},{},{},{},{},{},{},{}", + case.name, + case.satellite_radius_km, + BODY_RADIUS_KM, + elevation, + nadir, + central, + range, + footprint_central, + area, + fraction + ); + } + + // The inverse relation, written with its own header so the Python side + // can tell the two blocks apart by the first column. + println!(); + println!("inverse_central_degrees,elevation_rad,round_trip_central_degrees"); + for degrees in CENTRAL_ANGLES_DEGREES { + let central = degrees * PI / 180.0; + // Use the equatorial radius as the satellite radius for the inverse + // check, so the geometry is the degenerate one where the elevation is + // always zero. A more useful check uses a satellite altitude. + let satellite_radius = BODY_RADIUS_KM + 550.0; + let Ok(elevation) = + elevation_from_central_angle(satellite_radius, BODY_RADIUS_KM, central).ok_or(()) + else { + // Beyond the horizon the inverse has no solution, which is + // reported as an empty row rather than skipped, so the count of + // rows stays fixed. + println!("{degrees},,"); + continue; + }; + let Ok(back) = earth_central_angle(satellite_radius, BODY_RADIUS_KM, elevation).ok_or(()) + else { + println!("{degrees},,"); + continue; + }; + println!("{degrees},{elevation},{}", back * 180.0 / PI); + } +} diff --git a/validation/README.md b/validation/README.md index 13437fa..d1c10ec 100644 --- a/validation/README.md +++ b/validation/README.md @@ -32,7 +32,7 @@ validation/ │ ├── time/ # 时间尺度、儒略日对拍 │ ├── coord/ # 坐标变换对拍 │ ├── link/ # 链路预算对拍 -│ └── constellation/ # 星座覆盖对拍 +│ └── constellation/ # 覆盖几何对拍 └── data/ ├── tle/ # 公开 TLE 样本 └── expected/ # Rust 侧导出的待比对数据 @@ -90,6 +90,10 @@ cargo run -p apsis-rs --features physics --example dump_iau2000 \ # TEME → 地固系（SGP4 输出链路） cargo run -p apsis-rs --features physics --example dump_teme \ > validation/data/expected/rust_teme.csv + +# 覆盖几何（仰角、星下点角、地心角、斜距、覆盖面积） +cargo run -p apsis-rs --features network --example dump_coverage \ + > validation/data/expected/rust_coverage.csv ``` Python 侧读取同一批输入，用参考实现重算并输出误差表： @@ -107,6 +111,7 @@ pixi run python validation/python/atmosphere/check_p839.py 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 ``` 两个脚本都以退出码 0 表示通过、1 表示超差、2 表示缺少输入文件。 diff --git a/validation/python/constellation/.gitkeep b/validation/python/constellation/.gitkeep deleted file mode 100644 index e69de29..0000000 diff --git a/validation/python/constellation/check_coverage.py b/validation/python/constellation/check_coverage.py new file mode 100644 index 0000000..1a00aa2 --- /dev/null +++ b/validation/python/constellation/check_coverage.py @@ -0,0 +1,227 @@ +"""Cross-check the ApsisRS coverage geometry against spherical trigonometry. + +Reads the quantities written by ``cargo run -p apsis-rs --features network +--example dump_coverage`` and recomputes each one from the spherical +triangle relations 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 formula shows up as a disagreement. + +The three angles of the triangle formed by the satellite, the ground point, +and the centre of the Earth satisfy + + elevation + nadir + central = 90 degrees + +and the law of sines gives + + sin(nadir) / body_radius = cos(elevation) / satellite_radius + +from which the nadir angle, the central angle, and the slant range follow. +The footprint is the spherical cap the central angle cuts from the body, so +its area is a zone area and its fraction is that over the whole surface. + +The script needs no network access. + +Usage: + + pixi run python validation/python/constellation/check_coverage.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 form expressions. The tolerances are set to catch a +# transposed radius or a missing conversion, not to allow a loose match. +ANGLE_TOLERANCE_RAD = 1e-12 +LENGTH_TOLERANCE_KM = 1e-9 +AREA_RELATIVE_TOLERANCE = 1e-12 +FRACTION_TOLERANCE = 1e-12 + +DEFAULT_INPUT = Path("validation/data/expected/rust_coverage.csv") + + +def reference_nadir(satellite_radius: float, body_radius: float, elevation: float) -> float: + """Returns the nadir angle from the law of sines.""" + return math.asin(body_radius * math.cos(elevation) / satellite_radius) + + +def reference_central(satellite_radius: float, body_radius: float, elevation: float) -> float: + """Returns the central angle from the angle sum.""" + nadir = reference_nadir(satellite_radius, body_radius, elevation) + return max(0.0, math.pi / 2.0 - elevation - nadir) + + +def reference_range(satellite_radius: float, body_radius: float, elevation: float) -> float: + """Returns the slant range from the law of cosines.""" + central = reference_central(satellite_radius, body_radius, elevation) + squared = ( + satellite_radius**2 + + body_radius**2 + - 2.0 * satellite_radius * body_radius * math.cos(central) + ) + return math.sqrt(squared) + + +def reference_cap_area(central: float, body_radius: float) -> float: + """Returns the area of a spherical cap of a given central angle. + + A cap is a zone between the pole and a colatitude, so its area is + ``2 pi R^2 (1 - cos(central))``. The footprint is not centred on a pole in + general, but a spherical cap is rotationally symmetric, so the same + expression applies. + """ + return 2.0 * math.pi * body_radius**2 * (1.0 - math.cos(central)) + + +def read_blocks(path: Path): + """Reads the two blocks of the file. + + The first block is the case table with a header row; a blank line + separates it from the inverse relation block, which has its own header. + """ + 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 check_cases(rows: list[dict[str, str]]) -> tuple[bool, list[str]]: + """Checks the case table and returns whether it passed and any failures.""" + failures: list[str] = [] + compared = 0 + worst = { + "nadir": 0.0, + "central": 0.0, + "range": 0.0, + "area": 0.0, + "fraction": 0.0, + } + + for row in rows: + name = row["case"] + satellite_radius = float(row["satellite_radius_km"]) + body_radius = float(row["body_radius_km"]) + elevation = float(row["elevation_mask_rad"]) + compared += 1 + + nadir = reference_nadir(satellite_radius, body_radius, elevation) + central = reference_central(satellite_radius, body_radius, elevation) + slant = reference_range(satellite_radius, body_radius, elevation) + area = reference_cap_area(central, body_radius) + fraction = area / (4.0 * math.pi * body_radius**2) + + worst["nadir"] = max(worst["nadir"], abs(nadir - float(row["nadir_rad"]))) + worst["central"] = max(worst["central"], abs(central - float(row["central_rad"]))) + worst["range"] = max(worst["range"], abs(slant - float(row["slant_range_km"]))) + rust_area = float(row["footprint_area_km2"]) + worst["area"] = max(worst["area"], abs(area - rust_area) / rust_area) + worst["fraction"] = max( + worst["fraction"], abs(fraction - float(row["footprint_fraction"])) + ) + + if abs(nadir - float(row["nadir_rad"])) > ANGLE_TOLERANCE_RAD: + failures.append(f"{name}: nadir differs by {abs(nadir - float(row['nadir_rad']))}") + if abs(central - float(row["central_rad"])) > ANGLE_TOLERANCE_RAD: + failures.append( + f"{name}: central differs by {abs(central - float(row['central_rad']))}" + ) + if abs(slant - float(row["slant_range_km"])) > LENGTH_TOLERANCE_KM: + failures.append( + f"{name}: range differs by {abs(slant - float(row['slant_range_km']))}" + ) + if abs(area - rust_area) / rust_area > AREA_RELATIVE_TOLERANCE: + failures.append(f"{name}: area differs by {abs(area - rust_area) / rust_area}") + if abs(fraction - float(row["footprint_fraction"])) > FRACTION_TOLERANCE: + failures.append( + f"{name}: fraction differs by " + f"{abs(fraction - float(row['footprint_fraction']))}" + ) + + print(f"cases compared: {compared}") + print(f"worst nadir: {worst['nadir']:.6e} rad") + print(f"worst central: {worst['central']:.6e} rad") + print(f"worst range: {worst['range']:.6e} km") + print(f"worst area (rel): {worst['area']:.6e}") + print(f"worst fraction: {worst['fraction']:.6e}") + return (not failures, failures) + + +def check_inverse(rows: list[dict[str, str]]) -> tuple[bool, list[str]]: + """Checks that the central angle and elevation relations invert.""" + failures: list[str] = [] + compared = 0 + worst = 0.0 + for row in rows: + if not row["elevation_rad"]: + continue + central_degrees = float(row["inverse_central_degrees"]) + elevation = float(row["elevation_rad"]) + round_trip = float(row["round_trip_central_degrees"]) + compared += 1 + + # Recompute the round trip independently: from the elevation, the + # central angle follows from the angle sum. + satellite_radius = 6_378.137 + 550.0 + body_radius = 6_378.137 + nadir = reference_nadir(satellite_radius, body_radius, elevation) + expected = max(0.0, math.pi / 2.0 - elevation - nadir) + error = abs(expected * 180.0 / math.pi - round_trip) + worst = max(worst, error) + if error > 1e-9: + failures.append( + f"central {central_degrees} deg: round trip differs by {error}" + ) + + print(f"inverse compared: {compared}") + print(f"worst round trip: {worst:.6e} deg") + return (not failures, failures) + + +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 network " + "--example dump_coverage > " + 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 + + cases = list(csv.DictReader(blocks[0])) + inverse = list(csv.DictReader(blocks[1])) + if not cases or not inverse: + print("input file has an empty block", file=sys.stderr) + return 2 + + cases_ok, case_failures = check_cases(cases) + inverse_ok, inverse_failures = check_inverse(inverse) + + for failure in case_failures + inverse_failures: + print(f" {failure}") + + if cases_ok and inverse_ok: + print("result: pass") + return 0 + print("result: FAIL") + return 1 + + +if __name__ == "__main__": + raise SystemExit(main(sys.argv))