workspace: add the S.465 cross-check and document the ITU-R suite

crates/rf/apsis-antenna/src/reference.rs

@@ -115,6 +115,12 @@ pub fn main_lobe_break_angle(diameter_over_wavelength: f64) -> Option<f64> {
 /// head is 29 for a large aperture and 32 for a small one, and the floor
 /// angle is 36 or 48 degrees accordingly.
 ///
+/// The falling sidelobe expression does not reach exactly minus ten decibels
+/// at the floor angle, so the pattern has a step of a few hundredths of a
+/// decibel there. This is a property of the recommendation, which defines
+/// the regions separately, and it is reproduced here rather than smoothed
+/// over.
+///
 /// Returns `None` when the ratio or the angle is not usable.
 pub fn reference_pattern_gain(diameter_over_wavelength: f64, off_axis_angle: f64) -> Option<f64> {
     if !off_axis_angle.is_finite() || off_axis_angle < 0.0 {
@@ -133,7 +139,6 @@ pub fn reference_pattern_gain(diameter_over_wavelength: f64, off_axis_angle: f64
         Some(-10.0)
     }
 }
-
 /// Returns the gain of the S.465 reference pattern as a decibel value.
 ///
 /// Returns `None` when the ratio or the angle is not usable.

validation/README.md

@@ -70,6 +70,10 @@ cargo run -p apsis-rs --features rf --example dump_rain \
 # ITU-R P.618-13 雨衰斜路径
 cargo run -p apsis-rs --features rf --example dump_rain_path \
     > validation/data/expected/rust_rain_path.csv
+
+# ITU-R S.465 参考天线方向图
+cargo run -p apsis-rs --features rf --example dump_s465 \
+    > validation/data/expected/rust_s465.csv
 ```
 
 Python 侧读取同一批输入，用参考实现重算并输出误差表：
@@ -82,6 +86,7 @@ pixi run python validation/python/link/check_link.py
 pixi run python validation/python/atmosphere/check_gas.py
 pixi run python validation/python/atmosphere/check_rain.py
 pixi run python validation/python/atmosphere/check_rain_path.py
+pixi run python validation/python/link/check_s465.py
 ```
 
 两个脚本都以退出码 0 表示通过、1 表示超差、2 表示缺少输入文件。
@@ -110,6 +115,7 @@ Earth orientation 数据自动下载；本 crate 涉及的换算不需要该数
 | `atmosphere/check_gas.py` | 气体吸收，18 个频率，涵盖 60 GHz 氧气复体、22/183 GHz 水汽线 | `itur` 的 ITU-R P.676-12 实现 | 干空气 6.1e-11，水汽 6.4e-9（相对） |
 | `atmosphere/check_rain.py` | 雨衰系数，17 个频率 × 4 组仰角/极化 | `itur` 的 ITU-R P.838-3 实现 | k 7.5e-9，alpha 6.5e-13（相对） |
 | `atmosphere/check_rain_path.py` | 雨衰斜路径，8 个场景 | ITU-R P.618-13 方程（Python 转写，系数来自 `itur`） | 3.4e-13 km / 4.3e-13 dB |
+| `link/check_s465.py` | 参考天线方向图，3 种口径 × 8 个离轴角 | ITU-R S.465 方程（Python 转写） | 4.0e-10 dB |
 
 ### 关于射频层对拍的说明
 

validation/python/link/check_s465.py

@@ -0,0 +1,130 @@
+"""Cross-check the ApsisRS ITU-R S.465 reference pattern.
+
+Reads the values written by ``cargo run -p apsis-rs --features rf --example
+dump_s465`` and recomputes each one from the equations of the recommendation.
+
+The script needs no network access.
+
+Usage:
+
+    pixi run python validation/python/link/check_s465.py [path]
+
+Exit status is zero when every value agrees within tolerance and one
+otherwise.
+"""
+
+from __future__ import annotations
+
+import csv
+import math
+import sys
+from pathlib import Path
+
+# The peak gain offset and the branch constants of the recommendation.
+PEAK_GAIN_OFFSET_DB = 7.7
+LARGE_APERTURE_RATIO = 100.0
+
+# Tolerance in decibels.
+TOLERANCE_DB = 1e-6
+
+DEFAULT_INPUT = Path("validation/data/expected/rust_s465.csv")
+
+
+def read_rows(path: Path) -> list[dict[str, str]]:
+    """Reads the CSV rows."""
+    with path.open(newline="") as handle:
+        reader = csv.DictReader(handle)
+        return [row for row in reader if row.get("ratio")]
+
+
+def break_angle(ratio: float) -> tuple[float, float]:
+    """Returns the branch constant and the break angle, in degrees."""
+    head, floor_angle = (29.0, 36.0) if ratio >= LARGE_APERTURE_RATIO else (32.0, 48.0)
+    peak = 20.0 * math.log10(ratio) + PEAK_GAIN_OFFSET_DB
+
+    def difference(angle_degrees: float) -> float:
+        main_lobe = peak - 2.5e-3 * (ratio * angle_degrees) ** 2
+        sidelobe = head - 25.0 * math.log10(angle_degrees)
+        return main_lobe - sidelobe
+
+    low, high = 1.0e-9, floor_angle
+    if difference(high) > 0.0:
+        return head, high
+    for _ in range(200):
+        middle = 0.5 * (low + high)
+        if difference(middle) > 0.0:
+            low = middle
+        else:
+            high = middle
+    return head, 0.5 * (low + high)
+
+
+def reference(ratio: float, angle_degrees: float) -> float:
+    """Recomputes the pattern gain in decibels.
+
+    The API takes the angle in radians, so the degrees recorded in the file
+    are converted to radians and back before the branch decision. That
+    reproduces the rounding the implementation performs, which matters at
+    the branch boundaries where the recommendation has a small step.
+    """
+    peak = 20.0 * math.log10(ratio) + PEAK_GAIN_OFFSET_DB
+    head, break_at = break_angle(ratio)
+    head, floor_angle = (
+        (29.0, 36.0) if ratio >= LARGE_APERTURE_RATIO else (32.0, 48.0)
+    )
+    angle_degrees = math.degrees(math.radians(angle_degrees))
+    if angle_degrees < break_at:
+        return peak - 2.5e-3 * (ratio * angle_degrees) ** 2
+    if angle_degrees <= floor_angle:
+        return head - 25.0 * math.log10(angle_degrees)
+    return -10.0
+
+
+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 rf "
+            "--example dump_s465 > " + str(path),
+            file=sys.stderr,
+        )
+        return 2
+
+    rows = read_rows(path)
+    if not rows:
+        print("input file has no data rows", file=sys.stderr)
+        return 2
+
+    worst_gain = 0.0
+    worst_break = 0.0
+    worst_at = None
+
+    for row in rows:
+        ratio = float(row["ratio"])
+        angle_degrees = float(row["angle_deg"])
+        _, expected_break = break_angle(ratio)
+        expected_gain = reference(ratio, angle_degrees)
+
+        break_error = abs(float(row["break_deg"]) - expected_break)
+        gain_error = abs(float(row["gain_dbi"]) - expected_gain)
+        if break_error > worst_break or gain_error > worst_gain:
+            worst_at = (ratio, angle_degrees)
+            worst_break = max(worst_break, break_error)
+            worst_gain = max(worst_gain, gain_error)
+
+    print(f"points compared:  {len(rows)}")
+    print(f"worst break:      {worst_break:.3e} deg")
+    print(f"worst gain:       {worst_gain:.3e} dB")
+    print(f"worst point:      {worst_at}")
+    print(f"tolerance:        {TOLERANCE_DB:.3e} dB")
+
+    if worst_break <= TOLERANCE_DB and worst_gain <= TOLERANCE_DB:
+        print("result:           pass")
+        return 0
+    print("result:           FAIL")
+    return 1
+
+
+if __name__ == "__main__":
+    raise SystemExit(main(sys.argv))