diff --git a/crates/apsis-rs/Cargo.toml b/crates/apsis-rs/Cargo.toml index 7b26d0c..e28715c 100644 --- a/crates/apsis-rs/Cargo.toml +++ b/crates/apsis-rs/Cargo.toml @@ -165,3 +165,7 @@ required-features = ["rf"] [[example]] name = "dump_rain" required-features = ["rf"] + +[[example]] +name = "dump_rain_path" +required-features = ["rf"] diff --git a/crates/apsis-rs/examples/dump_rain_path.rs b/crates/apsis-rs/examples/dump_rain_path.rs new file mode 100644 index 0000000..bda3fff --- /dev/null +++ b/crates/apsis-rs/examples/dump_rain_path.rs @@ -0,0 +1,50 @@ +//! Writes rain path attenuation computed by this workspace, for +//! cross-checking the ITU-R P.618-13 slant path method. +//! +//! Usage: +//! +//! ```text +//! cargo run -p apsis-rs --features rf --example dump_rain_path \ +//! > validation/data/expected/rust_rain_path.csv +//! ``` + +use apsis_atmosphere::{DEFAULT_RAIN_HEIGHT_KM, rain_attenuation}; + +/// The scenarios sampled. Each is a rain rate in millimeters per hour, a +/// frequency in gigahertz, an elevation in degrees, and a latitude in +/// degrees. The station height is zero and the rain height is the default. +const SCENARIOS: &[(f64, f64, f64, f64)] = &[ + (20.0, 12.0, 30.0, 45.0), + (20.0, 12.0, 10.0, 45.0), + (20.0, 20.0, 30.0, 45.0), + (5.0, 12.0, 30.0, 45.0), + (50.0, 12.0, 30.0, 45.0), + (20.0, 12.0, 30.0, 10.0), + (20.0, 12.0, 60.0, 45.0), + (20.0, 6.0, 30.0, 45.0), +]; + +fn main() { + println!( + "rain_rate_mm_per_h,frequency_ghz,elevation_deg,latitude_deg,\ + specific_db_per_km,path_length_km,attenuation_db" + ); + for &(rain_rate, frequency_ghz, elevation_degrees, latitude) in SCENARIOS { + let elevation = elevation_degrees.to_radians(); + let Some(result) = rain_attenuation( + rain_rate, + frequency_ghz * 1.0e9, + 0.0, + elevation, + 0.0, + DEFAULT_RAIN_HEIGHT_KM, + latitude, + ) else { + continue; + }; + println!( + "{rain_rate},{frequency_ghz},{elevation_degrees},{latitude},{:.12},{:.12},{:.12}", + result.specific_db_per_km, result.path_length_km, result.attenuation_db + ); + } +} diff --git a/crates/rf/apsis-antenna/src/lib.rs b/crates/rf/apsis-antenna/src/lib.rs index b6f70af..3f649e1 100644 --- a/crates/rf/apsis-antenna/src/lib.rs +++ b/crates/rf/apsis-antenna/src/lib.rs @@ -27,9 +27,10 @@ //! - Half power beam width of a circular aperture: `theta = k lambda / D` //! with `k` about 1.02 for uniform illumination and close to 1.2 for //! typical illumination tapers. -//! - The idealised element pattern: ITU-R Recommendation M.1853, which gives -//! `G = -min(12 (theta / theta_0)^2, G_0)` for the co-polar pattern of a -//! simple Earth station antenna. +//! - The co-polar reference pattern of a simple Earth station antenna, which +//! falls as the square of the off-axis angle down to a sidelobe floor. This +//! is the form used by ITU-R Recommendation S.465 to bound interference +//! from one satellite system into another. //! //! ``` //! use apsis_antenna::ParabolicAperture; diff --git a/crates/rf/apsis-antenna/src/pattern.rs b/crates/rf/apsis-antenna/src/pattern.rs index 5a07155..82c0957 100644 --- a/crates/rf/apsis-antenna/src/pattern.rs +++ b/crates/rf/apsis-antenna/src/pattern.rs @@ -7,8 +7,9 @@ use core::f64::consts::PI; /// Returns the gain of an antenna at an angle from boresight, given a peak /// gain and a half power beam width. /// -/// The model is the one used for interference analysis in ITU-R -/// Recommendation M.1853 and adopted in the Radio Regulations. The gain +/// This is the squared angle model used for interference analysis, in the +/// form adopted by the Radio Regulations and by ITU-R Recommendation S.465 +/// for the reference pattern of a simple Earth station antenna. The gain /// falls as the square of the off-axis angle until it reaches the peak gain /// minus a fixed floor, and stays there beyond. In decibels, /// diff --git a/crates/rf/apsis-atmosphere/src/rain.rs b/crates/rf/apsis-atmosphere/src/rain.rs index 6193d8f..486752b 100644 --- a/crates/rf/apsis-atmosphere/src/rain.rs +++ b/crates/rf/apsis-atmosphere/src/rain.rs @@ -184,7 +184,12 @@ pub struct RainAttenuation { pub attenuation_db: f64, /// The horizontal reduction factor, dimensionless and at most one. pub horizontal_reduction: f64, - /// The vertical adjustment factor, dimensionless and at most one. + /// The vertical adjustment factor, dimensionless. + /// + /// The factor is the ratio of the effective path length to the longer + /// of the two geometric path estimates. It is at most one when the path + /// meets the rain layer at a steep angle and can exceed one at a shallow + /// angle, which is the behaviour Recommendation P.618 gives. pub vertical_adjustment: f64, } @@ -392,17 +397,22 @@ pub fn rain_attenuation( }; // Equation 9: the vertical adjustment factor. + // + // The exponential takes the elevation in degrees, not radians, because + // the coefficients of 31 and 0.45 were fitted with degrees. The sine + // takes the elevation in radians, as usual. let latitude_magnitude = latitude.abs(); let xi = if latitude_magnitude < 36.0 { 36.0 - latitude_magnitude } else { 0.0 }; + let elevation_degrees = elevation.to_degrees(); let vertical_adjustment = 1.0 / (1.0 + sine.sqrt() * (31.0 - * (1.0 - (-(elevation / (1.0 + xi))).exp()) + * (1.0 - (-(elevation_degrees / (1.0 + xi))).exp()) * (path_length * specific).sqrt() / frequency_ghz.powi(2) - 0.45)); diff --git a/validation/python/atmosphere/check_rain_path.py b/validation/python/atmosphere/check_rain_path.py new file mode 100644 index 0000000..30954e5 --- /dev/null +++ b/validation/python/atmosphere/check_rain_path.py @@ -0,0 +1,167 @@ +"""Cross-check the ApsisRS rain path against ITU-R P.618-13. + +Reads the values written by ``cargo run -p apsis-rs --features rf --example +dump_rain_path`` and recomputes each one from the equations of the +recommendation. The specific attenuation coefficients come from the ``itur`` +implementation of ITU-R P.838-3, and the path geometry and the two reduction +factors are transcribed from Recommendation P.618-13. + +The script needs no network access. + +Usage: + + pixi run python validation/python/atmosphere/check_rain_path.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 + +from itur.models.itu838 import _ITU838_3_ + +# Tolerance on the path length and the attenuation, in their own units. +LENGTH_TOLERANCE_KM = 1e-9 +ATTENUATION_TOLERANCE_DB = 1e-9 + +# The effective radius of the Earth used by the path geometry, in +# kilometers, from Recommendation P.618-13. +EFFECTIVE_EARTH_RADIUS_KM = 8500.0 + +# The rain layer height used by the Rust example, in kilometers. +RAIN_HEIGHT_KM = 3.0 + +DEFAULT_INPUT = Path("validation/data/expected/rust_rain_path.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("rain_rate_mm_per_h")] + + +def reference(row: dict[str, str], coefficients) -> tuple[float, float]: + """Recomputes the path length and the attenuation.""" + rain_rate = float(row["rain_rate_mm_per_h"]) + frequency = float(row["frequency_ghz"]) + elevation = float(row["elevation_deg"]) + latitude = float(row["latitude_deg"]) + station_height = 0.0 + + k, alpha = coefficients.rain_specific_attenuation_coefficients( + frequency, elevation, 0.0 + ) + specific = float(k) * rain_rate ** float(alpha) + + # Equation 2: the slant path length. + height_difference = RAIN_HEIGHT_KM - station_height + sine = math.sin(math.radians(elevation)) + slant_length = 2.0 * height_difference / ( + math.sqrt(sine * sine + 2.0 * height_difference / EFFECTIVE_EARTH_RADIUS_KM) + + sine + ) + + # Equation 3: the horizontal projection. + horizontal_length = abs(slant_length * math.cos(math.radians(elevation))) + + # Equation 6: the horizontal reduction factor. + horizontal_reduction = 1.0 / ( + 1.0 + + 0.78 * math.sqrt(horizontal_length * specific / frequency) + - 0.38 * (1.0 - math.exp(-2.0 * horizontal_length)) + ) + + # Equation 7: the angle at which the path meets the rain layer. + eta = math.degrees( + math.atan2(height_difference, horizontal_length * horizontal_reduction) + ) + + # Equation 8: the longer of the two path estimates. + if eta > elevation: + path_length = ( + horizontal_length * horizontal_reduction / math.cos(math.radians(elevation)) + ) + else: + path_length = height_difference / sine + + # Equation 9: the vertical adjustment factor. The exponential takes the + # elevation in degrees. + xi = 36.0 - abs(latitude) if abs(latitude) < 36.0 else 0.0 + vertical_adjustment = 1.0 / ( + 1.0 + + math.sqrt(sine) + * ( + 31.0 + * (1.0 - math.exp(-(elevation / (1.0 + xi)))) + * math.sqrt(path_length * specific) + / frequency**2 + - 0.45 + ) + ) + + effective_length = path_length * vertical_adjustment + return effective_length, specific * effective_length + + +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_rain_path > " + 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 + + coefficients = _ITU838_3_() + worst_length = 0.0 + worst_attenuation = 0.0 + worst_at = None + + for row in rows: + length, attenuation = reference(row, coefficients) + length_error = abs(float(row["path_length_km"]) - length) + attenuation_error = abs(float(row["attenuation_db"]) - attenuation) + if length_error > worst_length or attenuation_error > worst_attenuation: + worst_at = ( + row["rain_rate_mm_per_h"], + row["frequency_ghz"], + row["elevation_deg"], + row["latitude_deg"], + ) + worst_length = max(worst_length, length_error) + worst_attenuation = max(worst_attenuation, attenuation_error) + + print(f"scenarios compared: {len(rows)}") + print(f"worst path length: {worst_length:.3e} km") + print(f"worst attenuation: {worst_attenuation:.3e} dB") + print(f"worst scenario: {worst_at}") + print( + f"tolerance: {LENGTH_TOLERANCE_KM:.3e} km or " + f"{ATTENUATION_TOLERANCE_DB:.3e} dB" + ) + + if ( + worst_length <= LENGTH_TOLERANCE_KM + and worst_attenuation <= ATTENUATION_TOLERANCE_DB + ): + print("result: pass") + return 0 + print("result: FAIL") + return 1 + + +if __name__ == "__main__": + raise SystemExit(main(sys.argv))