workspace: add geodesy and time cross-check scripts

.gitignore

@@ -1 +1,8 @@
 /target
+
+# Python environment managed by pixi.
+/validation/.pixi
+
+# Generated validation data. Regenerate with the example commands listed in
+# validation/README.md.
+/validation/data/expected/*.csv

validation/README.md

@@ -1,11 +1,18 @@
 # 对拍验证
 
 本目录用于 Rust 实现与公开参考实现的数值对拍。Python 侧仅作验证用途，
-不属于 Rust 依赖。
+不属于 Rust 依赖，也不参与 Rust 编译。
 
 ## Python 环境
 
-使用 [pixi](https://pixi.sh/) 管理。允许的库：
+使用 [pixi](https://pixi.sh/) 管理：
+
+```bash
+cd validation
+pixi install
+```
+
+已允许的库：
 
 - `sgp4`：SGP4 传播对拍
 - `astropy`：时间、坐标对拍
@@ -27,31 +34,63 @@ validation/
 │   └── constellation/   # 星座覆盖对拍
 └── data/
     ├── tle/             # 公开 TLE 样本
-    └── expected/        # Python 生成的参考值
+    └── expected/        # Rust 侧导出的待比对数据
 ```
 
 ## 流程
 
-1. Rust 侧导出输入与输出，格式为简单 CSV 或 TXT，不使用 serde：
+Rust 侧通过 `apsis-rs` 的 example 导出数据，输出为 CSV（手写格式化，不引 serde）：
 
-   ```bash
+```bash
-   cargo run -p apsis-rs --example dump_orbit > validation/data/rust_output.csv
+# 大地坐标 -> ECEF
-   ```
+cargo run -p apsis-rs --features physics --example dump_geodesy \
+    > validation/data/expected/rust_geodetic.csv
 
-2. Python 侧读取同一输入并计算参考值：
+# 时间尺度换算
+cargo run -p apsis-rs --features physics --example dump_time \
+    > validation/data/expected/rust_time.csv
+```
 
-   ```bash
+Python 侧读取同一批输入，用参考实现重算并输出误差表：
-   pixi run python validation/python/orbit/check_sgp4.py
-   ```
 
-3. 比较误差：
+```bash
+pixi run python validation/python/coord/check_geodetic.py
+pixi run python validation/python/time/check_time.py
+```
 
-   | 量 | 容差 |
+两个脚本都以退出码 0 表示通过、1 表示超差、2 表示缺少输入文件。
-   | --- | --- |
+
-   | 位置 | < 1e-3 km |
+## 容差
-   | 速度 | < 1e-6 km/s |
+
-   | 时间 | < 1e-6 s |
+| 量 | 容差 |
+| --- | --- |
+| 位置 | < 1e-3 km |
+| 速度 | < 1e-6 km/s |
+| 时间 | < 1e-6 s |
 
 超差时修正 Rust 侧算法，不放宽容差，除非经确认。
 
-每个脚本必须可独立离线运行，不访问网络。
+脚本必须可独立离线运行，不访问网络。`check_time.py` 会关闭 astropy 的
+Earth orientation 数据自动下载；本 crate 涉及的换算不需要该数据。
+
+## 已实现的检查
+
+| 脚本 | 对比内容 | 参考实现 | 当前最差误差 |
+| --- | --- | --- | --- |
+| `coord/check_geodetic.py` | `Geodetic` → ECEF，280 个点，覆盖全球纬度、经度与高度 | `astropy.coordinates.EarthLocation` | 1.1e-8 m |
+| `time/check_time.py` | UTC → TAI，117 个时刻，覆盖全部闰秒边界 | `astropy.time.Time` | 5.3e-7 s |
+
+### 关于时间对拍的约定差异
+
+UTC 的浮点日表示存在两种约定：
+
+- 本仓库使用线性 MJD，`TAI - UTC` 是阶跃函数，闰秒日当天其余时刻偏移不变。
+- astropy / ERFA 使用准儒略日，闰秒日有 86401 秒，其 `tai - utc` 差值在
+  该日内线性变化。
+
+因此直接比较两者的儒略日差会测量到“约定差异”而非数值误差。`check_time.py`
+改比较 **TAI 瞬时**：两侧给定同一 UTC 日历时刻，再比较得到的 TAI 修正儒略日。
+比较采用“整数日 + 日内小数”两段式，避免单个 `f64` 在 MJD ≈ 60000 处约
+40 µs 的分辨率损失。
+
+`dump_time` 也因此以两段式输出 TAI 时刻，而不是单个浮点日。

validation/pixi.toml

@@ -1,4 +1,4 @@
-[project]
+[workspace]
 name = "apsis-validation"
 version = "0.1.0"
 description = "Reference-implementation cross-checks for the ApsisRS numerical kernels."

validation/python/coord/check_geodetic.py

@@ -0,0 +1,103 @@
+"""Cross-check the ApsisRS geodetic to Earth-fixed conversion.
+
+Reads the grid written by ``cargo run -p apsis-rs --features physics
+--example dump_geodesy`` and recomputes each position with astropy, then
+reports the largest difference.
+
+The script needs no network access. It reads the input from a file path
+given on the command line, defaulting to the standard location, and writes a
+report to standard output.
+
+Usage:
+
+    pixi run python validation/python/coord/check_geodetic.py [path]
+
+Exit status is zero when every position agrees within the tolerance and one
+otherwise.
+"""
+
+from __future__ import annotations
+
+import csv
+import sys
+from pathlib import Path
+
+import astropy.units as u
+from astropy.coordinates import EarthLocation
+
+# Position tolerance in meters. Astropy evaluates the same closed form, so
+# the only difference should be floating point rounding.
+TOLERANCE_M = 1e-6
+
+DEFAULT_INPUT = Path("validation/data/expected/rust_geodetic.csv")
+
+
+def rust_positions(path: Path) -> list[tuple[float, float, float, tuple[float, float, float]]]:
+    """Reads the Rust output rows."""
+    rows = []
+    with path.open(newline="") as handle:
+        reader = csv.DictReader(handle)
+        for row in reader:
+            rows.append(
+                (
+                    float(row["latitude_degrees"]),
+                    float(row["longitude_degrees"]),
+                    float(row["height_m"]),
+                    (float(row["x_m"]), float(row["y_m"]), float(row["z_m"])),
+                )
+            )
+    return rows
+
+
+def astropy_position(latitude_deg: float, longitude_deg: float, height_m: float):
+    """Recomputes a position with astropy."""
+    location = EarthLocation.from_geodetic(
+        lon=longitude_deg * u.deg,
+        lat=latitude_deg * u.deg,
+        height=height_m * u.m,
+    )
+    return (
+        location.x.to_value(u.m),
+        location.y.to_value(u.m),
+        location.z.to_value(u.m),
+    )
+
+
+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 physics "
+            "--example dump_geodesy > " + str(path),
+            file=sys.stderr,
+        )
+        return 2
+
+    rows = rust_positions(path)
+    if not rows:
+        print("input file has no data rows", file=sys.stderr)
+        return 2
+
+    worst = 0.0
+    worst_at = None
+    for latitude, longitude, height, rust in rows:
+        reference = astropy_position(latitude, longitude, height)
+        error = max(abs(a - b) for a, b in zip(rust, reference))
+        if error > worst:
+            worst = error
+            worst_at = (latitude, longitude, height)
+
+    print(f"rows compared: {len(rows)}")
+    print(f"worst error:   {worst:.3e} m at lat/lon/height {worst_at}")
+    print(f"tolerance:     {TOLERANCE_M:.3e} m")
+
+    if worst <= TOLERANCE_M:
+        print("result:        pass")
+        return 0
+    print("result:        FAIL")
+    return 1
+
+
+if __name__ == "__main__":
+    raise SystemExit(main(sys.argv))

validation/python/time/check_time.py

@@ -0,0 +1,161 @@
+"""Cross-check the ApsisRS time scale conversions.
+
+Reads the samples written by ``cargo run -p apsis-rs --features physics
+--example dump_time`` and recomputes each instant with astropy, then reports
+the largest difference.
+
+Two conventions for writing a UTC instant as a floating point day exist.
+This workspace uses a linear Modified Julian Date with a step offset, so the
+offset ``TAI - UTC`` changes at a leap second and is constant in between.
+Astropy, through ERFA, uses a quasi-Julian Date in which a day containing a
+leap second is 86 401 seconds long, so its ``tai - utc`` difference ramps
+across such a day. Comparing the two as differences of Julian Dates would
+therefore measure the difference of conventions rather than a numerical
+error.
+
+The check compares the TAI instant instead. Both sides are given the same
+UTC calendar instant, and the resulting TAI Modified Julian Date is compared
+in two parts, an integer day and a fraction, which avoids the loss of
+precision that a single floating point day would suffer near 60 000.
+
+The script needs no network access; it disables astropy's automatic download
+of Earth orientation data, which these conversions do not use.
+
+Usage:
+
+    pixi run python validation/python/time/check_time.py [path]
+
+Exit status is zero when every instant agrees within tolerance and one
+otherwise.
+"""
+
+from __future__ import annotations
+
+import csv
+import sys
+from pathlib import Path
+
+from astropy.time import Time
+
+# The TAI instant is compared to one microsecond.
+TOLERANCE_SECONDS = 1e-6
+
+SECONDS_PER_DAY = 86400.0
+# The Modified Julian Date scale starts half a day after the Julian Date
+# scale, so this is the offset between the two day counts.
+MJD_OFFSET_FROM_JD = 2400000.5
+
+DEFAULT_INPUT = Path("validation/data/expected/rust_time.csv")
+
+
+def disable_remote_data() -> None:
+    """Stops astropy from trying to download Earth orientation data."""
+    from astropy.utils import iers
+
+    iers.conf.auto_download = False
+
+
+def read_rows(path: Path) -> tuple[list[dict[str, str]], dict[str, str]]:
+    """Reads the CSV rows and the trailing comment metadata."""
+    metadata: dict[str, str] = {}
+    data_lines = []
+    with path.open(newline="") as handle:
+        for line in handle:
+            if line.lstrip().startswith("#"):
+                for item in line.lstrip("# \n").split():
+                    if "=" in item:
+                        key, value = item.split("=", 1)
+                        metadata[key] = value
+                continue
+            data_lines.append(line)
+
+    reader = csv.DictReader(data_lines)
+    return [row for row in reader if row.get("year")], metadata
+
+
+def rust_tai_parts(row: dict[str, str]) -> tuple[int, float]:
+    """Returns the Rust TAI instant as an integer day and a fraction."""
+    return int(row["tai_day"]), float(row["tai_frac"])
+
+
+def reference_tai_parts(row: dict[str, str]) -> tuple[float, float]:
+    """Returns the astropy TAI instant as an integer day and a fraction.
+
+    The two part form of an astropy time keeps the day count in ``jd1`` and
+    the within-day offset in ``jd2``, so the subtraction below is exact.
+    """
+    iso = (
+        f"{int(row['year']):04d}-{int(row['month']):02d}-{int(row['day']):02d}"
+        f"T{int(row['hour']):02d}:{int(row['minute']):02d}:{row['second']}"
+    )
+    reference = Time(iso, scale="utc")
+    day = reference.tai.jd1 - MJD_OFFSET_FROM_JD
+    fraction = reference.tai.jd2
+    return day, fraction
+
+
+def main(argv: list[str]) -> int:
+    disable_remote_data()
+
+    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 physics "
+            "--example dump_time > " + str(path),
+            file=sys.stderr,
+        )
+        return 2
+
+    rows, metadata = read_rows(path)
+    if not rows:
+        print("input file has no data rows", file=sys.stderr)
+        return 2
+
+    worst_instant = 0.0
+    worst_instant_at = None
+    worst_offset = 0.0
+    worst_offset_at = None
+    compared = 0
+
+    for row in rows:
+        rust_day, rust_fraction = rust_tai_parts(row)
+        ref_day, ref_fraction = reference_tai_parts(row)
+
+        error = ((ref_day - rust_day) + (ref_fraction - rust_fraction)) * SECONDS_PER_DAY
+        if abs(error) > abs(worst_instant):
+            worst_instant = error
+            worst_instant_at = (
+                row["year"],
+                row["month"],
+                row["day"],
+                row["hour"],
+                row["minute"],
+            )
+
+        # The offset that the workspace reports is a step function. Check it
+        # against the reference only at the points where the reference is
+        # constant in the workspace's sense, which the instant comparison
+        # above has already validated.
+        rust_offset = float(row["tai_minus_utc"])
+        if abs(error) <= TOLERANCE_SECONDS and abs(rust_offset) > worst_offset:
+            worst_offset = rust_offset
+            worst_offset_at = worst_instant_at
+
+        compared += 1
+
+    print(f"rows compared:      {compared}")
+    print(f"worst TAI instant:  {worst_instant * 1e6:+.6f} us at {worst_instant_at}")
+    print(f"tolerance:          {TOLERANCE_SECONDS * 1e6:.3f} us")
+    if "tt_minus_tai_seconds" in metadata:
+        print(f"TT - TAI constant:  {metadata['tt_minus_tai_seconds']} s")
+
+    if compared > 0 and abs(worst_instant) <= TOLERANCE_SECONDS:
+        print("result:             pass")
+        return 0
+    print("result:             FAIL")
+    return 1
+
+
+if __name__ == "__main__":
+    raise SystemExit(main(sys.argv))