ApsisRS
/
apsis-rs
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

rf: add atmospheric gaseous and rain attenuation models

This commit is contained in:
ApsisRS contributors 2026-09-26 14:42:47 +08:00
parent b3d3071778
commit 6629969150
9 changed files with 828 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
Cargo.lock generated
View File

@ -17,6 +17,10 @@ dependencies = [
[[package]]
name = "apsis-atmosphere"
version = "0.1.0"
dependencies = [
"apsis-core",
"apsis-rf",
]
[[package]]
name = "apsis-conflict"

2
crates/rf/apsis-atmosphere/Cargo.toml
View File

@ -11,3 +11,5 @@ publish = false
workspace = true
[dependencies]
apsis-core.workspace = true
apsis-rf.workspace = true

123
crates/rf/apsis-atmosphere/src/gaseous.rs Normal file
View File

@ -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<f64> {
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<f64> {
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<GaseousAttenuation> {
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()
}

45
crates/rf/apsis-atmosphere/src/geometry.rs Normal file
View File

@ -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
}
}

63
crates/rf/apsis-atmosphere/src/lib.rs
View File

@ -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;

244
crates/rf/apsis-atmosphere/src/rain.rs Normal file
View File

@ -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<f64> {
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<RainAttenuation> {
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
}

107
crates/rf/apsis-atmosphere/tests/gaseous.rs Normal file
View File

@ -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);
}

63
crates/rf/apsis-atmosphere/tests/geometry.rs Normal file
View File

@ -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);
}
}

177
crates/rf/apsis-atmosphere/tests/rain.rs Normal file
View File

@ -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);
}