"""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. 用球面三角交叉核对 ApsisRS 的覆盖几何。 读取由 ``cargo run -p apsis-rs --features network --example dump_coverage`` 写出的量，并只用标准库与 numpy 从球面三角形关系逐个重算。要点在于 Python 侧是依据几何而非依据 Rust 代码写成的，故 Rust 公式中的错误会表现为不一致。 由卫星、地面点与地心构成的三角形的三个角满足 elevation + nadir + central = 90 degrees 正弦定理给出 sin(nadir) / body_radius = cos(elevation) / satellite_radius 由此得出天底角、中心角与斜距。覆盖区是中心角从天体上截出的球冠，故其面积 是带面积，其比例即该面积除以整个表面。 本脚本不需要网络访问。 用法： pixi run python validation/python/constellation/check_coverage.py [path] 每个量都在容差内一致时退出状态为零，否则为一。 """ 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. 返回给定中心角的球冠面积。 球冠是极点与某余纬之间的带，故其面积为 ``2 pi R^2 (1 - cos(central))``。 覆盖区通常并不以极点为中心，但球冠具有旋转对称性，故同一表达式适用。 """ 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))