From 9b1f8df76b9e6de4f5f5af4c6cc59b904c78ab35 Mon Sep 17 00:00:00 2001 From: ApsisRS contributors Date: Sat, 26 Sep 2026 14:46:36 +0800 Subject: [PATCH] rf: add link budget cross-check example and script --- crates/apsis-rs/Cargo.toml | 4 + crates/apsis-rs/examples/dump_link.rs | 152 ++++++++++++++++++++++++++ validation/README.md | 18 +++ validation/python/link/check_link.py | 120 ++++++++++++++++++++ 4 files changed, 294 insertions(+) create mode 100644 crates/apsis-rs/examples/dump_link.rs create mode 100644 validation/python/link/check_link.py diff --git a/crates/apsis-rs/Cargo.toml b/crates/apsis-rs/Cargo.toml index 884146f..80512e6 100644 --- a/crates/apsis-rs/Cargo.toml +++ b/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 diff --git a/crates/apsis-rs/examples/dump_link.rs b/crates/apsis-rs/examples/dump_link.rs new file mode 100644 index 0000000..e36e2ca --- /dev/null +++ b/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| 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) + ); + } +} diff --git a/validation/README.md b/validation/README.md index 909db7e..26364ac 100644 --- a/validation/README.md +++ b/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 对拍的约定 diff --git a/validation/python/link/check_link.py b/validation/python/link/check_link.py new file mode 100644 index 0000000..c7ca661 --- /dev/null +++ b/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))