rf: add link budget cross-check example and script
This commit is contained in:
ApsisRS contributors
parent
c6064cc8e5
commit
9b1f8df76b
|
|
@ -26,6 +26,10 @@ required-features = ["physics"]
|
|||
name = "dump_orbit"
|
||||
required-features = ["orbit"]
|
||||
|
||||
[[example]]
|
||||
name = "dump_link"
|
||||
required-features = ["rf"]
|
||||
|
||||
[lints]
|
||||
workspace = true
|
||||
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
);
|
||||
}
|
||||
}
|
||||
|
|
@ -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 对拍的约定
|
||||
|
||||
|
|
|
|||
|
|
@ -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))
|
||||
Loading…
Reference in New Issue