rf: add the ITU-R P.839-4 rain height model

crates/apsis-rs/Cargo.toml

@@ -173,3 +173,7 @@ required-features = ["rf"]
 [[example]]
 name = "dump_s465"
 required-features = ["rf"]
+
+[[example]]
+name = "dump_p839"
+required-features = ["rf"]

crates/apsis-rs/examples/dump_p839.rs

@@ -0,0 +1,43 @@
+//! Writes the rain height computed by this workspace, for cross-checking
+//! against ITU-R P.839-4.
+//!
+//! Usage:
+//!
+//! ```text
+//! cargo run -p apsis-rs --features rf --example dump_p839 \
+//!     > validation/data/expected/rust_p839.csv
+//! ```
+
+use apsis_atmosphere::{rain_height, zero_isotherm_height};
+
+/// The locations sampled, as latitude and longitude in degrees.
+const POINTS: &[(f64, f64)] = &[
+    (0.0, 0.0),
+    (45.0, 0.0),
+    (-45.0, 90.0),
+    (51.5, -0.1),
+    (35.7, 139.7),
+    (-33.9, 151.2),
+    (28.6, -80.6),
+    (0.0, 359.9),
+    (0.0, -0.1),
+    (60.0, 30.0),
+    (-60.0, 300.0),
+    (23.5, 45.0),
+    (89.9, 0.0),
+    (-89.9, 180.0),
+    (10.0, 100.0),
+];
+
+fn main() {
+    println!("latitude_deg,longitude_deg,isotherm_km,rain_height_km");
+    for &(latitude, longitude) in POINTS {
+        let Some(isotherm) = zero_isotherm_height(latitude, longitude) else {
+            continue;
+        };
+        let Some(height) = rain_height(latitude, longitude) else {
+            continue;
+        };
+        println!("{latitude},{longitude},{isotherm:.12},{height:.12}");
+    }
+}

crates/rf/apsis-atmosphere/src/lib.rs

@@ -50,6 +50,8 @@
 mod gaseous;
 mod gaseous_path;
 mod geometry;
+mod p839;
+mod p839_table;
 mod rain;
 
 pub use gaseous::{
@@ -60,6 +62,7 @@ pub use gaseous_path::{
     EquivalentHeights, equivalent_heights, slant_path_attenuation, zenith_attenuation,
 };
 pub use geometry::{air_mass, elevation_from_zenith, zenith_to_elevation};
+pub use p839::{RAIN_HEIGHT_OFFSET_KM, grid_step_degrees, rain_height, zero_isotherm_height};
 pub use rain::{
     DEFAULT_RAIN_HEIGHT_KM, RainAttenuation, polarisation_tilt, rain_attenuation,
     rain_coefficients, specific_rain_attenuation,

crates/rf/apsis-atmosphere/src/p839.rs

@@ -0,0 +1,94 @@
+//! The rain height model of ITU-R Recommendation P.839.
+//!
+//! The height of the rain layer is taken as the height of the zero degree
+//! isotherm above sea level. The recommendation tabulates that height on a
+//! global grid, and the value at a given location is found by bilinear
+//! interpolation between the four surrounding grid points.
+//!
+//! The rain height is the height at which the falling snow melts, below
+//! which the precipitation is liquid. Above it, the attenuation is much
+//! smaller, so it is the effective thickness of the attenuating layer and an
+//! input to the path calculation of Recommendation P.618.
+//!
+//! Source: ITU-R Recommendation P.839-4, *Rain height model for prediction
+//! methods*.
+
+use crate::p839_table::ZERO_ISOTHERM_HEIGHT_METRES;
+
+/// The grid resolution of the tabulated isotherm height, in degrees.
+const STEP_DEGREES: f64 = 1.5;
+
+/// The number of rows in the table, running from north to south.
+const ROWS: usize = 121;
+
+/// The number of columns in the table, running from east to east.
+const COLUMNS: usize = 241;
+
+/// The offset added to the isotherm height to obtain the rain height, in
+/// kilometers.
+///
+/// The recommendation states that the rain height is the isotherm height
+/// plus 0.36 km, which represents the depth of the melting layer below the
+/// isotherm.
+pub const RAIN_HEIGHT_OFFSET_KM: f64 = 0.36;
+
+/// Returns the height of the zero degree isotherm at a location, in
+/// kilometers above mean sea level.
+///
+/// The latitude is in degrees, positive north, and the longitude in degrees,
+/// positive east. The value is interpolated bilinearly from the tabulated
+/// grid, and longitudes are reduced into `[0, 360)` first so that any
+/// convention is accepted.
+///
+/// Returns `None` when either coordinate is not finite, or when the latitude
+/// is outside the range of the grid.
+///
+/// Source: ITU-R Recommendation P.839-4, equation 1 and the tabulated
+/// isotherm height.
+pub fn zero_isotherm_height(latitude: f64, longitude: f64) -> Option<f64> {
+    if !latitude.is_finite() || !longitude.is_finite() {
+        return None;
+    }
+    if !(-90.0..=90.0).contains(&latitude) {
+        return None;
+    }
+
+    // The rows run from north to south, so the row index grows as the
+    // latitude falls.
+    let row_position = (90.0 - latitude) / STEP_DEGREES;
+    let row_lower = (row_position.floor() as usize).min(ROWS - 1);
+    let row_upper = (row_lower + 1).min(ROWS - 1);
+    let row_fraction = row_position - row_lower as f64;
+
+    // The columns run from 0 to 360 degrees east, wrapping at the seam.
+    let longitude = longitude.rem_euclid(360.0);
+    let column_position = longitude / STEP_DEGREES;
+    let column_lower = (column_position.floor() as usize).min(COLUMNS - 1);
+    let column_upper = (column_lower + 1) % COLUMNS;
+    let column_fraction = column_position - column_lower as f64;
+
+    let sample =
+        |row: usize, column: usize| f64::from(ZERO_ISOTHERM_HEIGHT_METRES[row][column]) / 1000.0;
+
+    // Bilinear interpolation.
+    let top = sample(row_lower, column_lower)
+        + (sample(row_lower, column_upper) - sample(row_lower, column_lower)) * column_fraction;
+    let bottom = sample(row_upper, column_lower)
+        + (sample(row_upper, column_upper) - sample(row_upper, column_lower)) * column_fraction;
+    Some(top + (bottom - top) * row_fraction)
+}
+
+/// Returns the rain height at a location, in kilometers above mean sea level.
+///
+/// This is the zero degree isotherm height plus
+/// [`RAIN_HEIGHT_OFFSET_KM`], as the recommendation defines it.
+///
+/// Returns `None` when the coordinates are not usable.
+pub fn rain_height(latitude: f64, longitude: f64) -> Option<f64> {
+    Some(zero_isotherm_height(latitude, longitude)? + RAIN_HEIGHT_OFFSET_KM)
+}
+
+/// Returns the grid resolution of the tabulated isotherm height, in degrees.
+pub const fn grid_step_degrees() -> f64 {
+    STEP_DEGREES
+}

crates/rf/apsis-atmosphere/tests/p839.rs

@@ -0,0 +1,147 @@
+//! Tests for the ITU-R P.839-4 rain height model.
+//!
+//! The expected values are those of the recommendation's tabulated grid,
+//! read at the grid points and interpolated between them. They were produced
+//! with the reference implementation of the recommendation.
+
+use apsis_atmosphere::{
+    RAIN_HEIGHT_OFFSET_KM, grid_step_degrees, rain_height, zero_isotherm_height,
+};
+
+/// The tolerance on an interpolated height, in kilometers.
+///
+/// The table is stored in metres, so the interpolation is exact to the
+/// resolution the recommendation publishes.
+const TOLERANCE: f64 = 1e-9;
+
+#[test]
+fn the_grid_resolution_matches_the_recommendation() {
+    assert_eq!(grid_step_degrees(), 1.5);
+}
+
+#[test]
+fn the_rain_height_offset_matches_the_recommendation() {
+    assert_eq!(RAIN_HEIGHT_OFFSET_KM, 0.36);
+}
+
+#[test]
+fn values_at_grid_points_match_the_recommendation() {
+    // These are the tabulated values at exact grid points, so no
+    // interpolation is involved.
+    let cases = [
+        (0.0, 0.0, 4.566),
+        (45.0, 0.0, 2.356),
+        (-45.0, 90.0, 1.158),
+        (90.0, 0.0, 2.096),
+        (-90.0, 0.0, 2.880),
+        (60.0, 30.0, 2.086),
+        (-60.0, 300.0, 0.544),
+    ];
+    for (latitude, longitude, expected) in cases {
+        let actual = zero_isotherm_height(latitude, longitude).expect("valid coordinates");
+        assert!(
+            (actual - expected).abs() < TOLERANCE,
+            "({latitude}, {longitude}): got {actual}, expected {expected}"
+        );
+    }
+}
+
+#[test]
+fn interpolated_values_match_the_recommendation() {
+    let cases = [
+        (35.7, 139.7, 3.866_48),
+        (-33.9, 151.2, 3.057_68),
+        (23.5, 45.0, 3.553_666_666_667),
+        (89.9, 0.0, 2.118_666_666_667),
+        (-89.9, 180.0, 2.895_933_333_333),
+        (10.0, 100.0, 4.676),
+    ];
+    for (latitude, longitude, expected) in cases {
+        let actual = zero_isotherm_height(latitude, longitude).expect("valid coordinates");
+        assert!(
+            (actual - expected).abs() < 1e-6,
+            "({latitude}, {longitude}): got {actual}, expected {expected}"
+        );
+    }
+}
+
+#[test]
+fn the_rain_height_is_the_isotherm_height_plus_the_offset() {
+    for (latitude, longitude) in [(0.0, 0.0), (45.0, 30.0), (-20.0, 200.0)] {
+        let isotherm = zero_isotherm_height(latitude, longitude).expect("valid coordinates");
+        let rain = rain_height(latitude, longitude).expect("valid coordinates");
+        assert!((rain - isotherm - RAIN_HEIGHT_OFFSET_KM).abs() < TOLERANCE);
+    }
+}
+
+#[test]
+fn longitude_wraps_at_the_seam() {
+    // The grid runs to 360 degrees, which is the same meridian as zero, so
+    // the height must be continuous across the seam.
+    let at_zero = zero_isotherm_height(10.0, 0.0).expect("valid");
+    let at_360 = zero_isotherm_height(10.0, 360.0).expect("valid");
+    assert!((at_zero - at_360).abs() < TOLERANCE);
+
+    // A negative longitude is the same as its positive equivalent.
+    let negative = zero_isotherm_height(10.0, -10.0).expect("valid");
+    let positive = zero_isotherm_height(10.0, 350.0).expect("valid");
+    assert!((negative - positive).abs() < TOLERANCE);
+}
+
+#[test]
+fn the_height_is_bounded_by_the_table() {
+    // Every tabulated value lies between the minimum and maximum of the
+    // recommendation's global grid, which are about 0.006 km and 6.281 km.
+    for latitude in [-90.0, -45.0, 0.0, 45.0, 90.0] {
+        for longitude in [0.0, 90.0, 180.0, 270.0] {
+            let height = zero_isotherm_height(latitude, longitude).expect("valid coordinates");
+            assert!(
+                (0.0..=7.0).contains(&height),
+                "({latitude}, {longitude}) gave {height}"
+            );
+        }
+    }
+}
+
+#[test]
+fn the_height_is_continuous_across_the_grid() {
+    // The interpolation must not jump between neighbouring grid cells.
+    let mut previous = zero_isotherm_height(0.0, 0.0).expect("valid");
+    for step in 1..=240 {
+        let longitude = step as f64 * 1.5;
+        let current = zero_isotherm_height(0.0, longitude).expect("valid coordinates");
+        assert!(
+            (current - previous).abs() < 1.0,
+            "jump at longitude {longitude}: {previous} to {current}"
+        );
+        previous = current;
+    }
+}
+
+#[test]
+fn the_tropics_are_higher_than_the_poles() {
+    // The isotherm is highest near the equator, where the atmosphere is
+    // warmest, and lowest at high latitudes.
+    let equator = zero_isotherm_height(0.0, 100.0).expect("valid");
+    let mid_latitude = zero_isotherm_height(50.0, 100.0).expect("valid");
+    let pole = zero_isotherm_height(80.0, 100.0).expect("valid");
+    assert!(equator > mid_latitude, "{equator} against {mid_latitude}");
+    assert!(mid_latitude > pole, "{mid_latitude} against {pole}");
+}
+
+#[test]
+fn invalid_coordinates_are_rejected() {
+    assert_eq!(zero_isotherm_height(f64::NAN, 0.0), None);
+    assert_eq!(zero_isotherm_height(0.0, f64::NAN), None);
+    assert_eq!(zero_isotherm_height(91.0, 0.0), None);
+    assert_eq!(zero_isotherm_height(-91.0, 0.0), None);
+    assert_eq!(rain_height(91.0, 0.0), None);
+    assert_eq!(rain_height(f64::INFINITY, 0.0), None);
+}
+
+#[test]
+fn the_poles_are_accepted() {
+    assert!(zero_isotherm_height(90.0, 0.0).is_some());
+    assert!(zero_isotherm_height(-90.0, 0.0).is_some());
+    assert!(zero_isotherm_height(90.0, 359.9).is_some());
+}

validation/README.md

@@ -74,6 +74,10 @@ cargo run -p apsis-rs --features rf --example dump_rain_path \
 # ITU-R S.465 参考天线方向图
 cargo run -p apsis-rs --features rf --example dump_s465 \
     > validation/data/expected/rust_s465.csv
+
+# ITU-R P.839-4 零度等温线高度
+cargo run -p apsis-rs --features rf --example dump_p839 \
+    > validation/data/expected/rust_p839.csv
 ```
 
 Python 侧读取同一批输入，用参考实现重算并输出误差表：
@@ -87,6 +91,7 @@ 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
+pixi run python validation/python/atmosphere/check_p839.py
 ```
 
 两个脚本都以退出码 0 表示通过、1 表示超差、2 表示缺少输入文件。
@@ -116,6 +121,7 @@ Earth orientation 数据自动下载；本 crate 涉及的换算不需要该数
 | `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 |
+| `atmosphere/check_p839.py` | 零度等温线高度，15 个全球点 | `itur` 的 ITU-R P.839-4 实现 | 3.3e-13 km |
 
 ### 关于射频层对拍的说明
 
@@ -126,7 +132,7 @@ Earth orientation 数据自动下载；本 crate 涉及的换算不需要该数
 - 其余以解析性质断言为主：各 crate 的单元测试验证模型必须满足的性质，例如 3 dB 点、
   平方律下降、Friis 级联行为、P.838-3 系数在表列频率处的取值。
 
-所有建议书均采用**现行版本**：P.676-12、P.838-3、P.618-13、S.465。
+所有建议书均采用**现行版本**：P.676-12、P.838-3、P.618-13、P.839-4、S.465。
 P.838-3 的系数是闭式曲线拟合，不是已被取代的 P.838-1 数值表。
 
 ### 关于 SGP4 对拍的约定

validation/python/atmosphere/check_p839.py

@@ -0,0 +1,105 @@
+"""Cross-check the ApsisRS rain height against ITU-R P.839-4.
+
+Reads the values written by ``cargo run -p apsis-rs --features rf --example
+dump_p839`` and recomputes each one with the ``itur`` package, which
+implements the recommendation directly.
+
+The reference implementation returns a not-a-number value outside the region
+its bilinear interpolator covers, which happens near the poles and at the
+seam. Those points are reported separately rather than counted as failures,
+because the tabulated grid does cover them and the Rust model interpolates
+from it.
+
+The script needs no network access.
+
+Usage:
+
+    pixi run python validation/python/atmosphere/check_p839.py [path]
+
+Exit status is zero when every compared value agrees within tolerance and one
+otherwise.
+"""
+
+from __future__ import annotations
+
+import csv
+import math
+import sys
+from pathlib import Path
+
+import numpy as np
+
+from itur.models.itu839 import _ITU839_4_
+
+# The table is stored in metres, so the interpolation reproduces the
+# tabulated values to that resolution.
+TOLERANCE_KM = 1e-9
+
+DEFAULT_INPUT = Path("validation/data/expected/rust_p839.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("latitude_deg")]
+
+
+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_p839 > " + 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
+
+    model = _ITU839_4_()
+    worst = 0.0
+    worst_at = None
+    compared = 0
+    outside = 0
+
+    for row in rows:
+        latitude = float(row["latitude_deg"])
+        longitude = float(row["longitude_deg"])
+        reference = float(
+            np.asarray(
+                model.isoterm_0(
+                    np.array([[latitude]]), np.array([[longitude]])
+                )
+            ).ravel()[0]
+        )
+        if math.isnan(reference):
+            # The reference interpolator does not cover this point.
+            outside += 1
+            continue
+
+        actual = float(row["isotherm_km"])
+        error = abs(actual - reference)
+        if error > worst:
+            worst = error
+            worst_at = (latitude, longitude)
+        compared += 1
+
+    print(f"points compared:      {compared}")
+    print(f"outside the reference: {outside}")
+    print(f"worst difference:     {worst:.3e} km at {worst_at}")
+    print(f"tolerance:            {TOLERANCE_KM:.3e} km")
+
+    if compared > 0 and worst <= TOLERANCE_KM:
+        print("result:               pass")
+        return 0
+    print("result:               FAIL")
+    return 1
+
+
+if __name__ == "__main__":
+    raise SystemExit(main(sys.argv))

validation/tools/generate_p839_table.py

@@ -0,0 +1,65 @@
+#!/usr/bin/env python3
+"""Generate the Rust source for the ITU-R P.839-4 rain height table.
+
+The recommendation gives the height of the zero degree isotherm as a
+tabulated global grid at 1.5 degree resolution. This script reads the grid
+from the ``itur`` reference package and emits a Rust module that embeds it as
+integer tenths of a kilometre, which is the resolution the recommendation
+publishes.
+
+Usage:
+
+    python validation/tools/generate_p839_table.py > /tmp/p839_table.rs
+"""
+
+from __future__ import annotations
+
+import sys
+
+import numpy as np
+
+BASE = (
+    "validation/.pixi/envs/default/lib/python3.12/site-packages/itur/data/839/"
+)
+
+# The grid resolution, in degrees.
+STEP_DEGREES = 1.5
+
+
+def main() -> int:
+    heights = np.load(BASE + "v4_esa0height.npz")["arr_0"]
+    latitudes = np.load(BASE + "v4_esalat.npz")["arr_0"]
+    longitudes = np.load(BASE + "v4_esalon.npz")["arr_0"]
+
+    rows = heights.shape[0]
+    columns = heights.shape[1]
+
+    # The grid is regular: latitude runs from 90 down to -90 and longitude
+    # from 0 to 360, both in steps of 1.5 degrees.
+    assert np.allclose(latitudes[:, 0], np.linspace(90.0, -90.0, rows)), "latitude grid"
+    assert np.allclose(longitudes[0, :], np.linspace(0.0, 360.0, columns)), "longitude grid"
+
+    # The recommendation publishes the height to the nearest metre, so the
+    # table is stored in metres to avoid losing precision.
+    metres = np.round(heights * 1000.0).astype(np.uint16)
+
+    print("/// The height of the zero degree isotherm, in metres.")
+    print("///")
+    print("/// The grid runs from 90 degrees north to 90 degrees south in the row")
+    print("/// direction and from 0 to 360 degrees east in the column direction,")
+    print("/// both in steps of 1.5 degrees, so there are 121 rows and 241 columns.")
+    print("/// The value at a fractional latitude or longitude is found by bilinear")
+    print("/// interpolation between the four surrounding points.")
+    print("///")
+    print("/// Source: ITU-R Recommendation P.839-4, *Rain height model for prediction")
+    print("/// methods*, which tabulates the isotherm height on this grid.")
+    print(f"pub static ZERO_ISOTHERM_HEIGHT_METRES: [[u16; {columns}]; {rows}] = [")
+    for row in range(rows):
+        values = ", ".join(str(int(value)) for value in metres[row])
+        print(f"    [{values}],")
+    print("];")
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())