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))