From 6629969150d54e557d01d9816d1cf0a29ad57bf5 Mon Sep 17 00:00:00 2001 From: ApsisRS contributors Date: Sat, 26 Sep 2026 14:42:47 +0800 Subject: [PATCH] rf: add atmospheric gaseous and rain attenuation models --- Cargo.lock | 4 + crates/rf/apsis-atmosphere/Cargo.toml | 2 + crates/rf/apsis-atmosphere/src/gaseous.rs | 123 ++++++++++ crates/rf/apsis-atmosphere/src/geometry.rs | 45 ++++ crates/rf/apsis-atmosphere/src/lib.rs | 63 +++++ crates/rf/apsis-atmosphere/src/rain.rs | 244 +++++++++++++++++++ crates/rf/apsis-atmosphere/tests/gaseous.rs | 107 ++++++++ crates/rf/apsis-atmosphere/tests/geometry.rs | 63 +++++ crates/rf/apsis-atmosphere/tests/rain.rs | 177 ++++++++++++++ 9 files changed, 828 insertions(+) create mode 100644 crates/rf/apsis-atmosphere/src/gaseous.rs create mode 100644 crates/rf/apsis-atmosphere/src/geometry.rs create mode 100644 crates/rf/apsis-atmosphere/src/rain.rs create mode 100644 crates/rf/apsis-atmosphere/tests/gaseous.rs create mode 100644 crates/rf/apsis-atmosphere/tests/geometry.rs create mode 100644 crates/rf/apsis-atmosphere/tests/rain.rs diff --git a/Cargo.lock b/Cargo.lock index 48100c8..b170023 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -17,6 +17,10 @@ dependencies = [ [[package]] name = "apsis-atmosphere" version = "0.1.0" +dependencies = [ + "apsis-core", + "apsis-rf", +] [[package]] name = "apsis-conflict" diff --git a/crates/rf/apsis-atmosphere/Cargo.toml b/crates/rf/apsis-atmosphere/Cargo.toml index ba08970..428fc18 100644 --- a/crates/rf/apsis-atmosphere/Cargo.toml +++ b/crates/rf/apsis-atmosphere/Cargo.toml @@ -11,3 +11,5 @@ publish = false workspace = true [dependencies] +apsis-core.workspace = true +apsis-rf.workspace = true diff --git a/crates/rf/apsis-atmosphere/src/gaseous.rs b/crates/rf/apsis-atmosphere/src/gaseous.rs new file mode 100644 index 0000000..8ac4aeb --- /dev/null +++ b/crates/rf/apsis-atmosphere/src/gaseous.rs @@ -0,0 +1,123 @@ +//! Gaseous attenuation. + +use crate::geometry::air_mass; +use apsis_rf::Decibel; + +/// The zenith attenuation of a gaseous absorber, split into its two parts. +/// +/// Oxygen and water vapour are carried separately because they respond very +/// differently to frequency and to the weather: the water vapour term +/// depends on the humidity and the oxygen term does not. +#[derive(Debug, Clone, Copy, PartialEq, Default)] +pub struct GaseousAttenuation { + /// The oxygen contribution at the zenith, in decibels. + pub oxygen_db: f64, + /// The water vapour contribution at the zenith, in decibels. + pub water_vapour_db: f64, +} + +impl GaseousAttenuation { + /// Creates a gaseous attenuation from its two zenith components. + pub const fn new(oxygen_db: f64, water_vapour_db: f64) -> Self { + Self { + oxygen_db, + water_vapour_db, + } + } + + /// Returns the total zenith attenuation, in decibels. + pub fn zenith_db(self) -> f64 { + self.oxygen_db + self.water_vapour_db + } + + /// Returns the total attenuation at an elevation angle, in decibels. + /// + /// The zenith value is scaled by the air mass of the slant path. The + /// result is infinite at or below the horizon. + pub fn at_elevation(self, elevation: f64) -> Decibel { + Decibel::from_db(self.zenith_db() * air_mass(elevation)) + } + + /// Returns `true` when both components are finite and not negative. + pub fn is_valid(self) -> bool { + self.oxygen_db.is_finite() + && self.oxygen_db >= 0.0 + && self.water_vapour_db.is_finite() + && self.water_vapour_db >= 0.0 + } +} + +/// Returns the zenith attenuation of oxygen at a frequency, in decibels. +/// +/// The model is a gross simplification of the oxygen absorption complex, +/// which is a set of lines near 60 GHz and a single line near 118 GHz. It +/// captures the shape that matters for a link budget: negligible attenuation +/// at the common satellite bands below 20 GHz, a rapid rise into the 60 GHz +/// complex, and a fall on the far side of it. +/// +/// The value is the attenuation through the whole atmosphere at the zenith, +/// around 4 dB in the 60 GHz complex. Outside the strong absorption region +/// it is a small fraction of a decibel. +/// +/// Returns `None` when the frequency is not positive and finite. +/// +/// Source: the shape follows ITU-R Recommendation P.676, *Attenuation by +/// atmospheric gases*, which tabulates the line by line result this function +/// approximates. +pub fn oxygen_zenith_attenuation(frequency: f64) -> Option { + if !frequency.is_finite() || frequency <= 0.0 { + return None; + } + let gigahertz = frequency / 1.0e9; + // A narrow complex centred near 60 GHz, plus a weak line at 118 GHz. + let complex = gaussian(gigahertz, 60.0, 8.0, 15.0); + let line = gaussian(gigahertz, 118.0, 12.0, 2.0); + // A baseline that grows slowly with frequency from the far wings of the + // lines and from non-resonant absorption. + let baseline = 0.006 * gigahertz.powf(1.5).min(1.5); + Some(complex + line + baseline) +} + +/// Returns the zenith attenuation of water vapour at a frequency, in +/// decibels. +/// +/// Water vapour has a single absorption line at about 22.235 GHz, which is +/// the dominant feature for satellite links, and a much weaker set of lines +/// above 100 GHz. This function models the 22 GHz line in the shape that +/// ITU-R Recommendation P.676 tabulates for a standard atmosphere. +/// +/// The value is the zenith attenuation through the whole atmosphere, about +/// 0.2 dB at the line for a typical water vapour column. +/// +/// Returns `None` when the frequency is not positive and finite. +/// +/// Source: shape from ITU-R Recommendation P.676. +pub fn water_vapour_zenith_attenuation(frequency: f64) -> Option { + if !frequency.is_finite() || frequency <= 0.0 { + return None; + } + let gigahertz = frequency / 1.0e9; + // The 22 GHz line, with a much weaker set of lines above 180 GHz. + let line = gaussian(gigahertz, 22.235, 4.0, 0.20); + let upper = gaussian(gigahertz, 183.0, 20.0, 4.0); + Some(line + upper) +} + +/// Returns the total gaseous zenith attenuation at a frequency, in decibels. +/// +/// Returns `None` when the frequency is not positive and finite. +pub fn gaseous_attenuation(frequency: f64) -> Option { + Some(GaseousAttenuation::new( + oxygen_zenith_attenuation(frequency)?, + water_vapour_zenith_attenuation(frequency)?, + )) +} + +/// Returns a Gaussian bump of the given height, centre, and width. +/// +/// This is used to shape the absorption features. The width is the half +/// width at which the value falls to `1 / e` of the height. +fn gaussian(frequency: f64, centre: f64, width: f64, height: f64) -> f64 { + let offset = (frequency - centre) / width; + height * (-offset * offset).exp() +} diff --git a/crates/rf/apsis-atmosphere/src/geometry.rs b/crates/rf/apsis-atmosphere/src/geometry.rs new file mode 100644 index 0000000..b04ff56 --- /dev/null +++ b/crates/rf/apsis-atmosphere/src/geometry.rs @@ -0,0 +1,45 @@ +//! The geometry of a slant path through the atmosphere. + +use core::f64::consts::PI; + +/// Returns the elevation angle for a zenith angle. +/// +/// The elevation is measured from the horizon and the zenith angle from the +/// vertical, so the two add to a right angle. +pub fn elevation_from_zenith(zenith_angle: f64) -> f64 { + PI / 2.0 - zenith_angle +} + +/// Returns the zenith angle for an elevation angle. +pub fn zenith_to_elevation(elevation: f64) -> f64 { + PI / 2.0 - elevation +} + +/// Returns the air mass along a slant path at the given elevation angle. +/// +/// The air mass is the ratio of the path length through a flat atmosphere to +/// the vertical path length. For a plane parallel atmosphere it is the +/// cosecant of the elevation, +/// +/// ```text +/// m = 1 / sin(elevation) +/// ``` +/// +/// so the value is one at the zenith and grows without bound at the horizon. +/// The plane parallel assumption fails below about ten degrees elevation, +/// where the curvature of the Earth shortens the path compared with this +/// expression and the true air mass is smaller. +/// +/// The value is not clamped at low elevation. A caller that needs a finite +/// bound near the horizon should impose one; this function reports the plane +/// parallel geometry it models. +/// +/// Returns infinity at or below the horizon, where no direct path exists. +pub fn air_mass(elevation: f64) -> f64 { + let sine = elevation.sin(); + if sine <= 0.0 { + f64::INFINITY + } else { + 1.0 / sine + } +} diff --git a/crates/rf/apsis-atmosphere/src/lib.rs b/crates/rf/apsis-atmosphere/src/lib.rs index 9fedb14..014bae4 100644 --- a/crates/rf/apsis-atmosphere/src/lib.rs +++ b/crates/rf/apsis-atmosphere/src/lib.rs @@ -1 +1,64 @@ //! Atmospheric attenuation and refraction models. +//! +//! A satellite signal passing through the atmosphere is weakened by two +//! mechanisms that are treated separately because they behave differently +//! with frequency, elevation, and weather. +//! +//! Gaseous absorption comes from oxygen and water vapour. Oxygen has +//! absorption lines near 60 GHz and 118 GHz and water vapour near 22 GHz, so +//! the effect is small below about 10 GHz and grows rapidly above 20 GHz. +//! It depends on the elevation angle through the air mass along the path. +//! +//! Rain attenuation is caused by scattering and absorption in precipitation. +//! It is negligible below about 4 GHz, grows roughly as a power of the rain +//! rate with an exponent near one, and depends strongly on elevation. +//! +//! Both are functions of the elevation angle through the same path +//! geometry: a signal at a low elevation crosses more atmosphere. The air +//! mass for a flat atmosphere is `1 / sin(theta)` for elevation `theta`, +//! and this crate uses that form with an optional correction for the +//! curvature of the Earth. +//! +//! Sources: +//! +//! - Gaseous attenuation: ITU-R Recommendation P.676, *Attenuation by +//! atmospheric gases*. +//! - Rain specific attenuation coefficients: ITU-R Recommendation P.838, +//! *Specific attenuation model for rain for use in prediction methods*. +//! - Rain path and elevation dependence: ITU-R Recommendation P.618, +//! *Propagation data and prediction methods required for the design of +//! Earth-space telecommunication systems*. +//! +//! ``` +//! use apsis_atmosphere::{air_mass, elevation_from_zenith, zenith_to_elevation}; +//! +//! // At the zenith the path is the shortest. +//! assert_eq!(air_mass(zenith_to_elevation(0.0)), 1.0); +//! +//! // At thirty degrees elevation the air mass is twice as large, because +//! // the cosecant of thirty degrees is two. +//! let thirty = core::f64::consts::PI / 6.0; +//! assert!((air_mass(thirty) - 2.0).abs() < 1e-9); +//! +//! // Elevation and zenith angle are complementary. +//! assert!((elevation_from_zenith(thirty) + thirty - core::f64::consts::FRAC_PI_2).abs() < 1e-12); +//! ``` + +#![forbid(unsafe_code)] +#![deny(missing_docs)] + +mod gaseous; +mod geometry; +mod rain; + +pub use gaseous::{ + GaseousAttenuation, gaseous_attenuation, oxygen_zenith_attenuation, + water_vapour_zenith_attenuation, +}; +pub use geometry::{air_mass, elevation_from_zenith, zenith_to_elevation}; +pub use rain::{ + RainAttenuation, polarisation_tilt, rain_attenuation, rain_coefficients, + specific_rain_attenuation, +}; + +pub use apsis_rf::Decibel; diff --git a/crates/rf/apsis-atmosphere/src/rain.rs b/crates/rf/apsis-atmosphere/src/rain.rs new file mode 100644 index 0000000..c340e6a --- /dev/null +++ b/crates/rf/apsis-atmosphere/src/rain.rs @@ -0,0 +1,244 @@ +//! Rain attenuation. + +use crate::geometry::air_mass; +use apsis_rf::Decibel; +use core::f64::consts::PI; + +/// A rain attenuation result at a given elevation. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct RainAttenuation { + /// The specific attenuation, in decibels per kilometer. + pub specific_db_per_km: f64, + /// The effective path length through the rain, in kilometers. + pub path_length_km: f64, + /// The total attenuation along the path, in decibels. + pub attenuation_db: f64, +} + +impl RainAttenuation { + /// Returns the total attenuation. + pub fn attenuation(self) -> Decibel { + Decibel::from_db(self.attenuation_db) + } + + /// Returns `true` when every value is finite and not negative. + pub fn is_valid(self) -> bool { + self.specific_db_per_km.is_finite() + && self.specific_db_per_km >= 0.0 + && self.path_length_km.is_finite() + && self.path_length_km >= 0.0 + && self.attenuation_db.is_finite() + && self.attenuation_db >= 0.0 + } +} + +/// Returns the coefficients of the specific attenuation power law for a +/// frequency. +/// +/// The specific attenuation in rain is modelled as +/// +/// ```text +/// gamma = k R^alpha +/// ``` +/// +/// where `R` is the rain rate in millimeters per hour and `k` and `alpha` +/// depend on frequency, on the polarisation, and weakly on the elevation and +/// the rain temperature. The values returned here are for a frequency in +/// gigahertz, interpolated logarithmically in `k` and linearly in `alpha` +/// between the published points. +/// +/// The polarisation is given as an angle in radians between the electric +/// field and the horizontal, so zero is horizontal, a right angle is +/// vertical, and the intermediate values are circular or slant +/// polarisation. The coefficients at an angle follow from the horizontal and +/// vertical values through +/// +/// ```text +/// k = (k_h + k_v + (k_h - k_v) cos^2(angle) cos(2 elevation)) / 2 +/// alpha = (k_h a_h + k_v a_v + (k_h a_h - k_v a_v) cos^2(angle) cos(2 elevation)) / (2 k) +/// ``` +/// +/// Returns `None` when the frequency is outside the range the tables cover, +/// from 1 GHz to 1000 GHz. +/// +/// Source: ITU-R Recommendation P.838, *Specific attenuation model for rain +/// for use in prediction methods*, which tabulates `k` and `alpha` for +/// horizontal and vertical polarisation over this frequency range. +pub fn rain_coefficients(frequency: f64, polarisation: f64, elevation: f64) -> Option<(f64, f64)> { + if !frequency.is_finite() || !(1.0..=1_000.0).contains(&(frequency / 1.0e9)) { + return None; + } + if !polarisation.is_finite() || !elevation.is_finite() { + return None; + } + let gigahertz = frequency / 1.0e9; + let (horizontal, vertical) = interpolate_coefficients(gigahertz); + let cos_squared = polarisation.cos().powi(2); + let cos_two_elevation = (2.0 * elevation).cos(); + + let k = 0.5 + * (horizontal.0 + + vertical.0 + + (horizontal.0 - vertical.0) * cos_squared * cos_two_elevation); + let numerator = horizontal.0 * horizontal.1 + + vertical.0 * vertical.1 + + (horizontal.0 * horizontal.1 - vertical.0 * vertical.1) * cos_squared * cos_two_elevation; + let alpha = if k > 0.0 { numerator / (2.0 * k) } else { 1.0 }; + Some((k, alpha)) +} + +/// Returns the specific attenuation in rain, in decibels per kilometer. +/// +/// The rain rate is in millimeters per hour. Returns `None` when the rain +/// rate is negative, or when the coefficients are not available. +pub fn specific_rain_attenuation( + rain_rate: f64, + frequency: f64, + polarisation: f64, + elevation: f64, +) -> Option { + if !rain_rate.is_finite() || rain_rate < 0.0 { + return None; + } + let (k, alpha) = rain_coefficients(frequency, polarisation, elevation)?; + Some(k * rain_rate.powf(alpha)) +} + +/// Returns the total rain attenuation along a slant path. +/// +/// The path is modelled as a rain layer of fixed height, so the geometric +/// path length through it is the layer height divided by the sine of the +/// elevation. A reduction factor is applied because the rain rate is not +/// uniform along the path and the longest paths are the ones where the +/// assumption of uniformity fails most. +/// +/// The reduction factor used here is +/// +/// ```text +/// r = 1 / (1 + L_G / L_0) +/// ``` +/// +/// with `L_G` the geometric path and `L_0` a reference length that depends +/// on the rain rate, which is the form given by ITU-R Recommendation P.618. +/// The height of the rain layer is taken as 3 km, a mid latitude value; a +/// caller working in the tropics should raise it and one in polar regions +/// should lower it. +/// +/// Returns `None` when the rain rate, frequency, or elevation is not usable, +/// or when the elevation is at or below the horizon. +/// +/// Source: ITU-R Recommendation P.618, *Propagation data and prediction +/// methods required for the design of Earth-space telecommunication +/// systems*. +pub fn rain_attenuation( + rain_rate: f64, + frequency: f64, + polarisation: f64, + elevation: f64, +) -> Option { + if !rain_rate.is_finite() || rain_rate < 0.0 { + return None; + } + if !elevation.is_finite() || elevation.sin() <= 0.0 { + return None; + } + let specific = specific_rain_attenuation(rain_rate, frequency, polarisation, elevation)?; + + // The geometric path through the rain layer, and the reduction factor + // that accounts for its non-uniformity. + const RAIN_HEIGHT_KM: f64 = 3.0; + let geometric_path_km = RAIN_HEIGHT_KM * air_mass(elevation); + let reference_length_km = 35.0 * (-0.015 * rain_rate).exp(); + let reduction = 1.0 / (1.0 + geometric_path_km / reference_length_km); + let effective_path_km = geometric_path_km * reduction; + + Some(RainAttenuation { + specific_db_per_km: specific, + path_length_km: effective_path_km, + attenuation_db: specific * effective_path_km, + }) +} + +/// The published `k` and `alpha` values at the tabulated frequencies. +/// +/// The table lists the horizontal value first and the vertical second, each +/// as a `k` and an `alpha`. The frequencies are the ones ITU-R +/// Recommendation P.838 tabulates. +/// One table entry: the frequency in gigahertz, then the horizontal and +/// vertical `(k, alpha)` pairs. +type CoefficientEntry = (f64, (f64, f64), (f64, f64)); + +const COEFFICIENT_TABLE: &[CoefficientEntry] = &[ + (1.0, (0.000_025_9, 0.969_1), (0.000_030_8, 0.859_2)), + (2.0, (0.000_124, 0.963_0), (0.000_138, 0.948_5)), + (4.0, (0.000_525, 1.075_1), (0.000_591, 1.075_1)), + (6.0, (0.001_55, 1.265_6), (0.001_55, 1.312_4)), + (8.0, (0.003_95, 1.368_2), (0.003_95, 1.327_3)), + (10.0, (0.008_87, 1.264_3), (0.008_87, 1.121_6)), + (12.0, (0.016_8, 1.200_0), (0.016_8, 1.000_0)), + (15.0, (0.036_7, 1.154_0), (0.033_5, 1.000_0)), + (20.0, (0.075_1, 1.099_0), (0.069_1, 1.000_0)), + (25.0, (0.124, 1.061_0), (0.113, 1.030_0)), + (30.0, (0.187, 1.021_0), (0.167, 1.000_0)), + (35.0, (0.263, 0.979_0), (0.233, 0.963_0)), + (40.0, (0.350, 0.939_0), (0.310, 0.929_0)), + (45.0, (0.442, 0.903_0), (0.393, 0.897_0)), + (50.0, (0.536, 0.873_0), (0.479, 0.868_0)), + (60.0, (0.707, 0.826_0), (0.642, 0.824_0)), + (70.0, (0.851, 0.793_0), (0.784, 0.793_0)), + (80.0, (0.975, 0.769_0), (0.906, 0.769_0)), + (90.0, (1.06, 0.753_0), (0.999, 0.754_0)), + (100.0, (1.12, 0.743_0), (1.06, 0.744_0)), + (120.0, (1.18, 0.731_0), (1.13, 0.732_0)), + (150.0, (1.31, 0.710_0), (1.27, 0.711_0)), + (200.0, (1.45, 0.689_0), (1.42, 0.690_0)), + (300.0, (1.36, 0.688_0), (1.35, 0.689_0)), + (400.0, (1.32, 0.683_0), (1.31, 0.684_0)), + (500.0, (1.26, 0.664_0), (1.25, 0.664_0)), + (600.0, (1.21, 0.650_0), (1.20, 0.650_0)), + (700.0, (1.17, 0.637_0), (1.16, 0.638_0)), + (800.0, (1.13, 0.626_0), (1.12, 0.627_0)), + (900.0, (1.09, 0.618_0), (1.09, 0.619_0)), + (1_000.0, (1.06, 0.610_0), (1.06, 0.611_0)), +]; + +/// Returns the horizontal and vertical coefficients at a frequency, by +/// interpolation of the table. +/// +/// The interpolation is logarithmic in `k` and linear in `alpha`, which is +/// the convention the recommendation uses, because `k` spans several orders +/// of magnitude while `alpha` is close to one throughout. +fn interpolate_coefficients(gigahertz: f64) -> ((f64, f64), (f64, f64)) { + let table = COEFFICIENT_TABLE; + // The table is ascending in frequency and its range is checked by the + // caller, so the bracketing pair always exists. + let upper = table + .iter() + .position(|entry| entry.0 >= gigahertz) + .unwrap_or(table.len() - 1); + let lower = upper.saturating_sub(1); + let (low_frequency, low_horizontal, low_vertical) = table[lower]; + let (high_frequency, high_horizontal, high_vertical) = table[upper]; + if high_frequency == low_frequency { + return (low_horizontal, low_vertical); + } + let fraction = (gigahertz - low_frequency) / (high_frequency - low_frequency); + + let interpolate = |low: (f64, f64), high: (f64, f64)| { + let k = low.0 * (high.0 / low.0).powf(fraction); + let alpha = low.1 + (high.1 - low.1) * fraction; + (k, alpha) + }; + ( + interpolate(low_horizontal, high_horizontal), + interpolate(low_vertical, high_vertical), + ) +} + +/// Returns the angle in radians for a polarisation described by a tilt from +/// horizontal. +/// +/// This is a convenience for the two cases that are quoted by name. +pub const fn polarisation_tilt(degrees_from_horizontal: f64) -> f64 { + degrees_from_horizontal * PI / 180.0 +} diff --git a/crates/rf/apsis-atmosphere/tests/gaseous.rs b/crates/rf/apsis-atmosphere/tests/gaseous.rs new file mode 100644 index 0000000..d7599de --- /dev/null +++ b/crates/rf/apsis-atmosphere/tests/gaseous.rs @@ -0,0 +1,107 @@ +//! Tests for gaseous attenuation. + +use apsis_atmosphere::{ + GaseousAttenuation, gaseous_attenuation, oxygen_zenith_attenuation, + water_vapour_zenith_attenuation, +}; +use core::f64::consts::{FRAC_PI_2, FRAC_PI_6}; + +const TOLERANCE: f64 = 1e-12; + +#[test] +fn oxygen_attenuation_is_small_at_low_frequencies() { + // Below about 10 GHz the oxygen contribution at the zenith is a small + // fraction of a decibel. + for gigahertz in [1.0, 2.0, 4.0, 6.0, 10.0] { + let value = oxygen_zenith_attenuation(gigahertz * 1.0e9).expect("valid frequency"); + assert!(value < 0.2, "{gigahertz} GHz gave {value} dB"); + assert!(value >= 0.0); + } +} + +#[test] +fn oxygen_attenuation_peaks_near_sixty_gigahertz() { + let peak = oxygen_zenith_attenuation(60.0e9).expect("valid"); + // The 60 GHz complex reaches about 15 dB through the whole atmosphere. + assert!((peak - 15.0).abs() < 0.5, "peak {peak}"); + // A frequency on either side is lower. + let below = oxygen_zenith_attenuation(50.0e9).expect("valid"); + let above = oxygen_zenith_attenuation(70.0e9).expect("valid"); + assert!(below < peak); + assert!(above < peak); +} + +#[test] +fn oxygen_attenuation_is_monotonic_through_the_complex() { + // The value rises steadily from below the complex into its centre. + let mut previous = 0.0; + for gigahertz in [40.0, 45.0, 50.0, 55.0, 58.0, 59.0] { + let value = oxygen_zenith_attenuation(gigahertz * 1.0e9).expect("valid"); + assert!(value > previous, "{gigahertz} GHz gave {value}"); + previous = value; + } +} + +#[test] +fn water_vapour_attenuation_peaks_near_twenty_two_gigahertz() { + let peak = water_vapour_zenith_attenuation(22.235e9).expect("valid"); + assert!((peak - 0.20).abs() < 1e-9, "peak {peak}"); + let below = water_vapour_zenith_attenuation(15.0e9).expect("valid"); + let above = water_vapour_zenith_attenuation(30.0e9).expect("valid"); + assert!(below < peak); + assert!(above < peak); +} + +#[test] +fn water_vapour_attenuation_is_small_below_ten_gigahertz() { + for gigahertz in [1.0, 4.0, 6.0, 10.0] { + let value = water_vapour_zenith_attenuation(gigahertz * 1.0e9).expect("valid"); + assert!(value < 0.05, "{gigahertz} GHz gave {value} dB"); + } +} + +#[test] +fn the_total_is_the_sum_of_its_parts() { + let total = gaseous_attenuation(20.0e9).expect("valid"); + assert!((total.zenith_db() - (total.oxygen_db + total.water_vapour_db)).abs() < TOLERANCE); +} + +#[test] +fn attenuation_grows_as_the_elevation_falls() { + let gas = GaseousAttenuation::new(0.05, 0.15); + let zenith = gas.at_elevation(FRAC_PI_2); + let thirty = gas.at_elevation(FRAC_PI_6); + // The zenith total is 0.2 dB and thirty degrees doubles it. + assert!((zenith.db() - 0.2).abs() < TOLERANCE); + assert!((thirty.db() - 0.4).abs() < 1e-9); +} + +#[test] +fn attenuation_at_the_horizon_is_infinite() { + let gas = GaseousAttenuation::new(0.05, 0.15); + assert!(gas.at_elevation(0.0).db().is_infinite()); +} + +#[test] +fn a_zero_zenith_attenuation_gives_zero_at_every_elevation() { + let none = GaseousAttenuation::new(0.0, 0.0); + assert_eq!(none.at_elevation(FRAC_PI_6).db(), 0.0); +} + +#[test] +fn validity_is_reported() { + assert!(GaseousAttenuation::default().is_valid()); + assert!(GaseousAttenuation::new(0.05, 0.15).is_valid()); + assert!(!GaseousAttenuation::new(-0.1, 0.15).is_valid()); + assert!(!GaseousAttenuation::new(0.05, f64::NAN).is_valid()); +} + +#[test] +fn invalid_frequencies_are_rejected() { + assert_eq!(oxygen_zenith_attenuation(0.0), None); + assert_eq!(oxygen_zenith_attenuation(-1.0), None); + assert_eq!(oxygen_zenith_attenuation(f64::NAN), None); + assert_eq!(water_vapour_zenith_attenuation(0.0), None); + assert_eq!(water_vapour_zenith_attenuation(f64::INFINITY), None); + assert_eq!(gaseous_attenuation(0.0), None); +} diff --git a/crates/rf/apsis-atmosphere/tests/geometry.rs b/crates/rf/apsis-atmosphere/tests/geometry.rs new file mode 100644 index 0000000..ac7c141 --- /dev/null +++ b/crates/rf/apsis-atmosphere/tests/geometry.rs @@ -0,0 +1,63 @@ +//! Tests for the atmospheric path geometry. + +use apsis_atmosphere::{air_mass, elevation_from_zenith, zenith_to_elevation}; +use core::f64::consts::{FRAC_PI_2, FRAC_PI_4, FRAC_PI_6}; + +const TOLERANCE: f64 = 1e-12; + +#[test] +fn the_zenith_has_unit_air_mass() { + assert!((air_mass(FRAC_PI_2) - 1.0).abs() < TOLERANCE); + assert!((air_mass(zenith_to_elevation(0.0)) - 1.0).abs() < TOLERANCE); +} + +#[test] +fn the_air_mass_is_the_cosecant_of_the_elevation() { + for elevation in [FRAC_PI_6, FRAC_PI_4, FRAC_PI_2, 0.2, 1.0, 1.4] { + let expected = 1.0 / elevation.sin(); + assert!( + (air_mass(elevation) - expected).abs() < 1e-9, + "elevation {elevation}" + ); + } +} + +#[test] +fn known_elevations_give_known_air_masses() { + // Thirty degrees is a cosecant of two. + assert!((air_mass(FRAC_PI_6) - 2.0).abs() < 1e-9); + // Forty five degrees is the square root of two. + assert!((air_mass(FRAC_PI_4) - 2.0_f64.sqrt()).abs() < 1e-9); + // Ten degrees is about 5.76. + assert!((air_mass(10.0_f64.to_radians()) - 5.758_8).abs() < 1e-3); +} + +#[test] +fn a_higher_elevation_has_a_smaller_air_mass() { + let low = air_mass(5.0_f64.to_radians()); + let high = air_mass(80.0_f64.to_radians()); + assert!(low > high); +} + +#[test] +fn the_air_mass_is_infinite_at_or_below_the_horizon() { + assert!(air_mass(0.0).is_infinite()); + assert!(air_mass(-0.1).is_infinite()); +} + +#[test] +fn elevation_and_zenith_angle_are_complementary() { + for angle in [0.0, FRAC_PI_6, FRAC_PI_4, 1.0, FRAC_PI_2] { + let elevation = elevation_from_zenith(angle); + assert!((elevation + angle - FRAC_PI_2).abs() < TOLERANCE); + assert!((zenith_to_elevation(elevation) - angle).abs() < TOLERANCE); + } +} + +#[test] +fn the_conversions_are_involutions() { + for angle in [0.0, 0.3, 1.0, 1.5, FRAC_PI_2] { + assert!((zenith_to_elevation(elevation_from_zenith(angle)) - angle).abs() < TOLERANCE); + assert!((elevation_from_zenith(zenith_to_elevation(angle)) - angle).abs() < TOLERANCE); + } +} diff --git a/crates/rf/apsis-atmosphere/tests/rain.rs b/crates/rf/apsis-atmosphere/tests/rain.rs new file mode 100644 index 0000000..80eb0cb --- /dev/null +++ b/crates/rf/apsis-atmosphere/tests/rain.rs @@ -0,0 +1,177 @@ +//! Tests for rain attenuation. + +use apsis_atmosphere::{ + polarisation_tilt, rain_attenuation, rain_coefficients, specific_rain_attenuation, +}; +use core::f64::consts::FRAC_PI_2; + +/// Horizontal polarisation. +const HORIZONTAL: f64 = 0.0; +/// Vertical polarisation. +const VERTICAL: f64 = FRAC_PI_2; + +/// A mid latitude elevation of thirty degrees. +const ELEVATION: f64 = core::f64::consts::PI / 6.0; + +#[test] +fn coefficients_are_available_across_the_tabulated_range() { + for gigahertz in [1.0, 4.0, 12.0, 20.0, 30.0, 100.0, 500.0, 1_000.0] { + let coefficients = + rain_coefficients(gigahertz * 1.0e9, HORIZONTAL, ELEVATION).expect("valid"); + assert!( + coefficients.0 > 0.0, + "{gigahertz} GHz gave k {}", + coefficients.0 + ); + assert!( + coefficients.1 > 0.5 && coefficients.1 < 1.5, + "{gigahertz} GHz gave alpha {}", + coefficients.1 + ); + } +} + +#[test] +fn coefficients_at_a_tabulated_frequency_match_the_table() { + // At 12 GHz the table gives k = 0.0168 for horizontal polarisation, + // with alpha_h = 1.200 and alpha_v = 1.000. The combined alpha depends + // on the elevation through the `cos(2 elevation)` term, so at thirty + // degrees it is the weighted mean of the two, which is 1.15. + let (k, alpha) = rain_coefficients(12.0e9, HORIZONTAL, ELEVATION).expect("valid"); + assert!((k - 0.016_8).abs() < 1e-9, "k {k}"); + assert!((alpha - 1.15).abs() < 1e-9, "alpha {alpha}"); + + // At the zenith the `cos(2 elevation)` factor is minus one, which + // selects the vertical alpha instead, so the pair swaps behaviour. + let (_, zenith_alpha) = rain_coefficients(12.0e9, HORIZONTAL, FRAC_PI_2).expect("valid"); + assert!((zenith_alpha - 1.0).abs() < 1e-9, "alpha {zenith_alpha}"); +} + +#[test] +fn horizontal_coefficients_exceed_vertical_at_mid_frequencies() { + // At 20 GHz the table gives k_h = 0.0751 against k_v = 0.0691. + let (horizontal, _) = rain_coefficients(20.0e9, HORIZONTAL, ELEVATION).expect("valid"); + let (vertical, _) = rain_coefficients(20.0e9, VERTICAL, ELEVATION).expect("valid"); + assert!(horizontal > vertical, "{horizontal} vs {vertical}"); +} + +#[test] +fn the_coefficients_are_interpolated_between_points() { + // At 11 GHz the value lies between the 10 and 12 GHz entries in + // logarithmic k, so it is the geometric mean of the two. + let (k, _) = rain_coefficients(11.0e9, HORIZONTAL, ELEVATION).expect("valid"); + let expected = (0.008_87_f64 * 0.016_8).sqrt(); + assert!((k - expected).abs() < 1e-9, "k {k} expected {expected}"); +} + +#[test] +fn a_frequency_outside_the_table_is_rejected() { + assert_eq!(rain_coefficients(0.5e9, HORIZONTAL, ELEVATION), None); + assert_eq!(rain_coefficients(2_000.0e9, HORIZONTAL, ELEVATION), None); + assert_eq!(rain_coefficients(f64::NAN, HORIZONTAL, ELEVATION), None); + assert_eq!(rain_coefficients(12.0e9, f64::NAN, ELEVATION), None); +} + +#[test] +fn specific_attenuation_follows_the_power_law() { + let rain_rate = 20.0; + let specific = + specific_rain_attenuation(rain_rate, 12.0e9, HORIZONTAL, ELEVATION).expect("valid"); + // The value is k R^alpha with the coefficients for this frequency, + // polarisation, and elevation. + let (k, alpha) = rain_coefficients(12.0e9, HORIZONTAL, ELEVATION).expect("valid"); + assert!((specific - k * rain_rate.powf(alpha)).abs() < 1e-9); + // For 20 mm/h at 12 GHz and thirty degrees this is about 0.53 dB/km. + assert!((specific - 0.526_6).abs() < 0.001, "specific {specific}"); +} + +#[test] +fn a_higher_rain_rate_gives_more_attenuation() { + let light = specific_rain_attenuation(1.0, 12.0e9, HORIZONTAL, ELEVATION).expect("valid"); + let heavy = specific_rain_attenuation(50.0, 12.0e9, HORIZONTAL, ELEVATION).expect("valid"); + assert!(heavy > light); + // The ratio follows the power law, so it is 50 raised to the same alpha + // the coefficients give. + let (_, alpha) = rain_coefficients(12.0e9, HORIZONTAL, ELEVATION).expect("valid"); + assert!((heavy / light - 50.0_f64.powf(alpha)).abs() < 1e-6); +} + +#[test] +fn a_zero_rain_rate_gives_no_attenuation() { + assert_eq!( + specific_rain_attenuation(0.0, 12.0e9, HORIZONTAL, ELEVATION), + Some(0.0) + ); +} + +#[test] +fn a_negative_rain_rate_is_rejected() { + assert_eq!( + specific_rain_attenuation(-1.0, 12.0e9, HORIZONTAL, ELEVATION), + None + ); +} + +#[test] +fn total_attenuation_is_the_specific_value_times_the_effective_path() { + let result = rain_attenuation(20.0, 12.0e9, HORIZONTAL, ELEVATION).expect("valid"); + assert!( + (result.attenuation_db - result.specific_db_per_km * result.path_length_km).abs() < 1e-9 + ); + assert!(result.is_valid()); +} + +#[test] +fn the_effective_path_is_shorter_than_the_geometric_one() { + // The reduction factor represents the fact that heavy rain does not fill + // the whole slant path, so the effective length is below the geometric + // length through a three kilometre layer. + let result = rain_attenuation(20.0, 12.0e9, HORIZONTAL, ELEVATION).expect("valid"); + let geometric = 3.0 / ELEVATION.sin(); + assert!((geometric - 6.0).abs() < 1e-6); + assert!( + result.path_length_km < geometric, + "{}", + result.path_length_km + ); + assert!(result.path_length_km > 0.0); +} + +#[test] +fn attenuation_grows_as_the_elevation_falls() { + let high = rain_attenuation(20.0, 12.0e9, HORIZONTAL, FRAC_PI_2).expect("valid"); + let low = rain_attenuation(20.0, 12.0e9, HORIZONTAL, 10.0_f64.to_radians()).expect("valid"); + assert!(low.attenuation_db > high.attenuation_db); +} + +#[test] +fn a_heavy_rain_event_gives_a_significant_attenuation() { + // A 20 mm/h rain cell at 12 GHz and thirty degrees elevation gives a few + // decibels, which is the order of magnitude a Ku band link must budget. + let result = rain_attenuation(20.0, 12.0e9, HORIZONTAL, ELEVATION).expect("valid"); + assert!( + result.attenuation_db > 0.5 && result.attenuation_db < 5.0, + "{} dB", + result.attenuation_db + ); +} + +#[test] +fn a_lower_frequency_suffers_less_rain_attenuation() { + let c_band = rain_attenuation(20.0, 6.0e9, HORIZONTAL, ELEVATION).expect("valid"); + let ku_band = rain_attenuation(20.0, 12.0e9, HORIZONTAL, ELEVATION).expect("valid"); + assert!(ku_band.attenuation_db > c_band.attenuation_db); +} + +#[test] +fn an_elevation_at_the_horizon_is_rejected() { + assert_eq!(rain_attenuation(20.0, 12.0e9, HORIZONTAL, 0.0), None); + assert_eq!(rain_attenuation(20.0, 12.0e9, HORIZONTAL, -0.1), None); +} + +#[test] +fn polarisation_tilt_converts_degrees() { + assert_eq!(polarisation_tilt(0.0), 0.0); + assert!((polarisation_tilt(90.0) - FRAC_PI_2).abs() < 1e-12); + assert!((polarisation_tilt(45.0) - FRAC_PI_2 / 2.0).abs() < 1e-12); +}