rf: add link budget cross-check example and script

crates/apsis-rs/Cargo.toml

@@ -26,6 +26,10 @@ required-features = ["physics"]
 name = "dump_orbit"
 required-features = ["orbit"]
 
+[[example]]
+name = "dump_link"
+required-features = ["rf"]
+
 [lints]
 workspace = true
 

crates/apsis-rs/examples/dump_link.rs

@@ -0,0 +1,152 @@
+//! Writes link budget inputs and the figures this workspace computes for
+//! them.
+//!
+//! The output is CSV on standard output, with a header row. It is consumed
+//! by `validation/python/link/check_link.py`, which recomputes the same
+//! quantities from the link equation and compares.
+//!
+//! Usage:
+//!
+//! ```text
+//! cargo run -p apsis-rs --features rf --example dump_link \
+//!     > validation/data/expected/rust_link.csv
+//! ```
+
+use apsis_link::{LinkBudget, LinkGeometry};
+use apsis_rf::{Decibel, Gain, Power};
+
+/// One link scenario, named so that the fields stay readable at each use.
+///
+/// The unit of each field is in its name: powers and gains in decibels,
+/// temperature in kelvin, and everything else in SI base units.
+struct Scenario {
+    transmit_dbw: f64,
+    transmit_dbi: f64,
+    transmit_loss_db: f64,
+    receive_dbi: f64,
+    receive_loss_db: f64,
+    temperature_k: f64,
+    bandwidth_hz: f64,
+    distance_m: f64,
+    frequency_hz: f64,
+}
+
+/// The scenarios checked. They cover a Ku band downlink, an S band link, a
+/// wideband Ka link, and an L band link with a warmer receiver.
+const SCENARIOS: &[Scenario] = &[
+    Scenario {
+        transmit_dbw: 10.0,
+        transmit_dbi: 30.0,
+        transmit_loss_db: 1.0,
+        receive_dbi: 35.0,
+        receive_loss_db: 0.5,
+        temperature_k: 150.0,
+        bandwidth_hz: 36.0e6,
+        distance_m: 38.0e6,
+        frequency_hz: 12.0e9,
+    },
+    Scenario {
+        transmit_dbw: 0.0,
+        transmit_dbi: 20.0,
+        transmit_loss_db: 2.0,
+        receive_dbi: 25.0,
+        receive_loss_db: 1.0,
+        temperature_k: 290.0,
+        bandwidth_hz: 1.0e6,
+        distance_m: 1_000_000.0,
+        frequency_hz: 2.2e9,
+    },
+    Scenario {
+        transmit_dbw: 20.0,
+        transmit_dbi: 45.0,
+        transmit_loss_db: 0.5,
+        receive_dbi: 50.0,
+        receive_loss_db: 0.2,
+        temperature_k: 80.0,
+        bandwidth_hz: 500.0e6,
+        distance_m: 400_000.0,
+        frequency_hz: 20.0e9,
+    },
+    Scenario {
+        transmit_dbw: 5.0,
+        transmit_dbi: 15.0,
+        transmit_loss_db: 3.0,
+        receive_dbi: 10.0,
+        receive_loss_db: 2.0,
+        temperature_k: 500.0,
+        bandwidth_hz: 10.0e6,
+        distance_m: 36_000_000.0,
+        frequency_hz: 1.5e9,
+    },
+];
+
+fn main() {
+    println!(
+        "transmit_dbw,transmit_dbi,transmit_loss_db,receive_dbi,receive_loss_db,\
+         temperature_k,bandwidth_hz,distance_m,frequency_hz,eirp_dbw,carrier_dbw,\
+         cn0_db,cn_db,ebn0_db"
+    );
+
+    for scenario in SCENARIOS {
+        let Scenario {
+            transmit_dbw,
+            transmit_dbi,
+            transmit_loss_db: transmit_loss,
+            receive_dbi,
+            receive_loss_db: receive_loss,
+            temperature_k: temperature,
+            bandwidth_hz: bandwidth,
+            distance_m: distance,
+            frequency_hz: frequency,
+        } = *scenario;
+
+        let Some(transmit_power) = Power::from_dbw(transmit_dbw) else {
+            continue;
+        };
+        let Some(transmit_gain) = Gain::from_dbi(transmit_dbi) else {
+            continue;
+        };
+        let Some(receive_gain) = Gain::from_dbi(receive_dbi) else {
+            continue;
+        };
+
+        let budget = LinkBudget {
+            transmit_power,
+            transmit_gain,
+            transmit_loss: Decibel::from_db(transmit_loss),
+            receive_gain,
+            receive_loss: Decibel::from_db(receive_loss),
+            system_temperature: apsis_link::NoiseTemperature::from_kelvin(temperature),
+            bandwidth,
+            geometry: LinkGeometry {
+                distance,
+                frequency,
+            },
+        };
+
+        let eirp = budget.eirp().and_then(|power| power.dbw());
+        let carrier = budget.received_power().and_then(|power| power.dbw());
+        let cn0 = budget.carrier_to_noise_density().map(|value| value.db());
+        let cn = budget.carrier_to_noise().map(|value| value.db());
+        // The energy per bit ratio is reported at a fixed data rate equal to
+        // the bandwidth, which makes it the carrier to noise ratio.
+        let ebn0 = budget
+            .energy_per_bit_over_noise(bandwidth)
+            .map(|value| value.db());
+
+        let field = |value: Option<f64>| match value {
+            Some(value) => format!("{value:.12}"),
+            None => "nan".to_string(),
+        };
+        println!(
+            "{transmit_dbw},{transmit_dbi},{transmit_loss},{receive_dbi},{receive_loss},\
+             {temperature},{bandwidth:.6},{distance:.6},{frequency:.6},\
+             {},{},{},{},{}",
+            field(eirp),
+            field(carrier),
+            field(cn0),
+            field(cn),
+            field(ebn0)
+        );
+    }
+}

validation/README.md

@@ -53,6 +53,10 @@ cargo run -p apsis-rs --features physics --example dump_time \
 # SGP4/SDP4 传播
 cargo run -p apsis-rs --features orbit --example dump_orbit \
     > validation/data/expected/rust_sgp4.csv
+
+# 链路预算
+cargo run -p apsis-rs --features rf --example dump_link \
+    > validation/data/expected/rust_link.csv
 ```
 
 Python 侧读取同一批输入，用参考实现重算并输出误差表：
@@ -61,6 +65,7 @@ Python 侧读取同一批输入，用参考实现重算并输出误差表：
 pixi run python validation/python/coord/check_geodetic.py
 pixi run python validation/python/time/check_time.py
 pixi run python validation/python/orbit/check_sgp4.py
+pixi run python validation/python/link/check_link.py
 ```
 
 两个脚本都以退出码 0 表示通过、1 表示超差、2 表示缺少输入文件。
@@ -85,6 +90,19 @@ 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 |
 | `orbit/check_sgp4.py` | SGP4/SDP4 传播，4 组 TLE（近地、高偏心、深空、GTO）× 9 个时刻 | `sgp4` 参考库 | 位置 1.2e-8 km，速度 5.3e-12 km/s |
+| `link/check_link.py` | 链路预算，4 个场景（Ku 下行、S 波段、Ka 宽带、L 波段） | 链路方程闭式解（Python 手写） | 5.0e-13 dB |
+
+### 关于射频层对拍的说明
+
+射频层没有可直接对拍的权威 Python 库，因此采用两种策略：
+
+- 闭式解手写对拍：链路预算的每一步都是初等运算，`check_link.py` 用 Python 独立
+  重写同一方程，常数取自 SI 定义，与 Rust 侧一致。
+- 解析性质断言：各 crate 的单元测试验证模型必须满足的性质，例如 3 dB 点、平方
+  律下降、Friis 公式的级联行为、以及 ITU-R P.838 系数表在表列频率处的取值。
+
+ITU-R 建议书（P.676、P.838、P.618）与 ITU-R M.1853 的公式均以数值表或闭式
+形式写在源码注释中，未引入外部数据文件。
 
 ### 关于 SGP4 对拍的约定
 

validation/python/link/check_link.py

@@ -0,0 +1,120 @@
+"""Cross-check the ApsisRS link budget against the link equation.
+
+Reads the scenarios written by ``cargo run -p apsis-rs --features rf
+--example dump_link`` and recomputes each figure from the link equation in
+plain Python, then reports the largest difference in each column.
+
+The reference here is the link equation itself rather than a library, because
+the equation is closed form and short enough to transcribe without
+ambiguity. The constants are taken from the 2019 SI definitions, so they
+match the values the Rust side uses.
+
+The script needs no network access.
+
+Usage:
+
+    pixi run python validation/python/link/check_link.py [path]
+
+Exit status is zero when every figure agrees within tolerance and one
+otherwise.
+"""
+
+from __future__ import annotations
+
+import csv
+import math
+import sys
+from pathlib import Path
+
+# Defining constants, exact in the SI.
+BOLTZMANN = 1.380649e-23
+SPEED_OF_LIGHT = 299792458.0
+
+# Tolerance in decibels. The two sides evaluate the same expressions with
+# double precision, so the only difference should be rounding.
+TOLERANCE_DB = 1e-9
+
+DEFAULT_INPUT = Path("validation/data/expected/rust_link.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("transmit_dbw")]
+
+
+def reference(row: dict[str, str]) -> dict[str, float]:
+    """Recomputes a budget from the link equation."""
+    transmit_dbw = float(row["transmit_dbw"])
+    transmit_dbi = float(row["transmit_dbi"])
+    transmit_loss = float(row["transmit_loss_db"])
+    receive_dbi = float(row["receive_dbi"])
+    receive_loss = float(row["receive_loss_db"])
+    temperature = float(row["temperature_k"])
+    bandwidth = float(row["bandwidth_hz"])
+    distance = float(row["distance_m"])
+    frequency = float(row["frequency_hz"])
+
+    eirp_dbw = transmit_dbw + transmit_dbi - transmit_loss
+    wavelength = SPEED_OF_LIGHT / frequency
+    path_loss_db = 20.0 * math.log10(4.0 * math.pi * distance / wavelength)
+    carrier_dbw = eirp_dbw - path_loss_db + receive_dbi - receive_loss
+    noise_density_dbw = 10.0 * math.log10(BOLTZMANN * temperature)
+    cn0_db = carrier_dbw - noise_density_dbw
+    cn_db = cn0_db - 10.0 * math.log10(bandwidth)
+    # The example reports the energy per bit ratio at a data rate equal to
+    # the bandwidth.
+    ebn0_db = cn0_db - 10.0 * math.log10(bandwidth)
+
+    return {
+        "eirp_dbw": eirp_dbw,
+        "carrier_dbw": carrier_dbw,
+        "cn0_db": cn0_db,
+        "cn_db": cn_db,
+        "ebn0_db": ebn0_db,
+    }
+
+
+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_link > " + 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
+
+    columns = ["eirp_dbw", "carrier_dbw", "cn0_db", "cn_db", "ebn0_db"]
+    worst = {column: 0.0 for column in columns}
+    worst_at = {column: None for column in columns}
+
+    for row in rows:
+        expected = reference(row)
+        for column in columns:
+            difference = abs(float(row[column]) - expected[column])
+            if difference > worst[column]:
+                worst[column] = difference
+                worst_at[column] = row["distance_m"]
+
+    print(f"scenarios compared: {len(rows)}")
+    for column in columns:
+        print(f"{column:>12}: worst {worst[column]:.3e} dB at distance {worst_at[column]} m")
+    print(f"{'tolerance':>12}:       {TOLERANCE_DB:.3e} dB")
+
+    if all(value <= TOLERANCE_DB for value in worst.values()):
+        print("result:       pass")
+        return 0
+    print("result:       FAIL")
+    return 1
+
+
+if __name__ == "__main__":
+    raise SystemExit(main(sys.argv))