rf: replace gaseous and rain models with ITU-R P.676-12 and P.838-3

crates/apsis-rs/Cargo.toml

@@ -157,3 +157,11 @@ apsis-scenario = { workspace = true, optional = true }
 apsis-event = { workspace = true, optional = true }
 apsis-output = { workspace = true, optional = true }
 apsis-sim = { workspace = true, optional = true }
+
+[[example]]
+name = "dump_gas"
+required-features = ["rf"]
+
+[[example]]
+name = "dump_rain"
+required-features = ["rf"]

crates/apsis-rs/examples/dump_gas.rs

@@ -0,0 +1,42 @@
+//! Writes the specific and zenith gaseous attenuation computed by this
+//! workspace, for cross-checking against ITU-R P.676-12.
+//!
+//! Usage:
+//!
+//! ```text
+//! cargo run -p apsis-rs --features rf --example dump_gas \
+//!     > validation/data/expected/rust_gas.csv
+//! ```
+
+use apsis_atmosphere::{
+    Atmosphere, equivalent_heights, specific_attenuation_dry_air, specific_attenuation_water_vapour,
+};
+
+/// The frequencies sampled, in gigahertz. They span the oxygen complex near
+/// 60 GHz, the water vapour lines at 22.235 and 183.31 GHz, and the windows
+/// between them.
+const FREQUENCIES_GHZ: &[f64] = &[
+    1.0, 2.0, 4.0, 6.0, 10.0, 12.0, 15.0, 20.0, 22.235, 30.0, 40.0, 50.0, 60.0, 70.0, 100.0, 118.0,
+    183.0, 300.0,
+];
+
+fn main() {
+    println!("frequency_ghz,dry_db_per_km,wet_db_per_km,dry_height_km,wet_height_km,zenith_db");
+    let atmosphere = Atmosphere::REFERENCE;
+    for &frequency in FREQUENCIES_GHZ {
+        let Some(dry) = specific_attenuation_dry_air(frequency, atmosphere) else {
+            continue;
+        };
+        let Some(wet) = specific_attenuation_water_vapour(frequency, atmosphere) else {
+            continue;
+        };
+        let Some(heights) = equivalent_heights(frequency, atmosphere) else {
+            continue;
+        };
+        let zenith = dry * heights.dry_km + wet * heights.wet_km;
+        println!(
+            "{frequency},{dry:.12},{wet:.12},{:.9},{:.9},{zenith:.12}",
+            heights.dry_km, heights.wet_km
+        );
+    }
+}

crates/apsis-rs/examples/dump_rain.rs

@@ -0,0 +1,36 @@
+//! Writes the P.838-3 rain coefficients computed by this workspace, for
+//! cross-checking against ITU-R P.838-3.
+//!
+//! Usage:
+//!
+//! ```text
+//! cargo run -p apsis-rs --features rf --example dump_rain \
+//!     > validation/data/expected/rust_rain.csv
+//! ```
+
+use apsis_atmosphere::rain_coefficients;
+
+/// The frequencies sampled, in gigahertz.
+const FREQUENCIES_GHZ: &[f64] = &[
+    1.0, 2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 15.0, 20.0, 22.0, 25.0, 30.0, 40.0, 50.0, 60.0, 70.0,
+    100.0,
+];
+
+fn main() {
+    println!("frequency_ghz,elevation_deg,polarisation_deg,k,alpha");
+    for &frequency in FREQUENCIES_GHZ {
+        // A horizontal path, which removes the elevation dependence and
+        // isolates the coefficient curves.
+        for &(elevation_degrees, polarisation_degrees) in
+            &[(90.0_f64, 0.0_f64), (90.0, 90.0), (30.0, 0.0), (30.0, 45.0)]
+        {
+            let elevation = elevation_degrees.to_radians();
+            let polarisation = polarisation_degrees.to_radians();
+            let Some((k, alpha)) = rain_coefficients(frequency * 1.0e9, polarisation, elevation)
+            else {
+                continue;
+            };
+            println!("{frequency},{elevation_degrees},{polarisation_degrees},{k:.12},{alpha:.12}");
+        }
+    }
+}

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

@@ -0,0 +1,170 @@
+//! Slant path attenuation through the atmosphere.
+//!
+//! The attenuation along a slant path is not simply the specific attenuation
+//! divided by the sine of the elevation, because the atmosphere thins with
+//! height and the path is curved. ITU-R Recommendation P.676 reduces the
+//! problem to two equivalent heights, one for dry air and one for water
+//! vapour, so that the zenith attenuation is the specific attenuation at the
+//! surface multiplied by that height.
+//!
+//! The equivalent height of dry air depends on frequency through the shape
+//! of the oxygen complex, and the equivalent height of water vapour depends
+//! on frequency and on the vapour density. Both depend on the surface
+//! pressure and temperature.
+//!
+//! Source: ITU-R Recommendation P.676-12, equations 30 through 38.
+
+use crate::gaseous::{Atmosphere, specific_attenuation_dry_air, specific_attenuation_water_vapour};
+
+/// The coefficients of the dry air equivalent height, table 3 of P.676-12.
+///
+/// Each entry is the coefficient and the line frequency in gigahertz.
+const DRY_HEIGHT_TERMS: &[(f64, f64)] = &[
+    (0.159_7, 118.750_334),
+    (0.106_6, 368.498_246),
+    (0.132_5, 424.763_020),
+    (0.124_2, 487.249_273),
+    (0.093_8, 715.392_902),
+    (0.144_8, 773.839_490),
+    (0.137_4, 834.145_546),
+];
+
+/// The coefficients of the water vapour equivalent height, table 4 of
+/// P.676-12.
+///
+/// Each entry is the line frequency in gigahertz, a coefficient, and a
+/// width-like parameter.
+const WET_HEIGHT_TERMS: &[(f64, f64, f64)] = &[
+    (22.235_08, 1.52, 2.56),
+    (183.310_087, 7.62, 10.2),
+    (325.152_888, 1.56, 2.7),
+    (380.197_353, 4.15, 5.7),
+    (439.150_807, 0.2, 0.91),
+    (448.001_085, 1.63, 2.46),
+    (474.689_092, 0.76, 2.22),
+    (488.490_108, 0.26, 2.49),
+    (556.935_985, 7.81, 10.0),
+    (620.700_87, 1.25, 2.35),
+    (752.033_113, 16.2, 20.0),
+    (916.171_582, 1.47, 2.58),
+    (970.315_022, 1.36, 2.44),
+    (987.926_764, 1.6, 1.86),
+];
+
+/// The equivalent height of dry air and of water vapour, in kilometers.
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub struct EquivalentHeights {
+    /// The equivalent height of dry air, in kilometers.
+    pub dry_km: f64,
+    /// The equivalent height of water vapour, in kilometers.
+    pub wet_km: f64,
+}
+
+/// Returns the equivalent heights of dry air and water vapour.
+///
+/// The heights are the lengths that, multiplied by the surface specific
+/// attenuation, give the zenith attenuation. The dry air height is limited
+/// to a pressure dependent bound below 70 GHz, where the oxygen complex is
+/// narrow enough that the height would otherwise be overestimated.
+///
+/// The frequency is in gigahertz.
+///
+/// Returns `None` when the frequency or atmosphere is not usable.
+///
+/// Source: ITU-R Recommendation P.676-12, equations 30 through 38.
+pub fn equivalent_heights(frequency_ghz: f64, atmosphere: Atmosphere) -> Option<EquivalentHeights> {
+    if !frequency_ghz.is_finite() || frequency_ghz <= 0.0 || !atmosphere.is_valid() {
+        return None;
+    }
+    let temperature = atmosphere.temperature_k;
+    let rho = atmosphere.water_vapour_density;
+    let e = atmosphere.water_vapour_partial_pressure();
+    let rp = (atmosphere.pressure_hpa + e) / 1_013.25;
+
+    // Equation 31, the resonance term centred just below 60 GHz.
+    let term1 = 5.104_0 / (1.0 + 0.066 * rp.powf(-2.3))
+        * (-((frequency_ghz - 59.7) / (2.87 + 12.4 * (-7.9 * rp).exp())).powi(2)).exp();
+
+    // Equation 32, the residual contribution of the submillimetre lines.
+    let term2: f64 = DRY_HEIGHT_TERMS
+        .iter()
+        .map(|&(coefficient, line_frequency)| {
+            coefficient * (2.12 * rp).exp()
+                / ((frequency_ghz - line_frequency).powi(2) + 0.025 * (2.2 * rp).exp())
+        })
+        .sum();
+
+    // Equation 33, the far wing correction.
+    let term3 = 0.011_4 * frequency_ghz / (1.0 + 0.14 * rp.powf(-2.6))
+        * (15.02 * frequency_ghz.powi(2) - 1_353.0 * frequency_ghz + 5.333e4)
+        / (frequency_ghz.powi(3) - 151.3 * frequency_ghz.powi(2) + 9_629.0 * frequency_ghz
+            - 6_803.0);
+
+    let a = 0.783_2 + 0.007_09 * (temperature - 273.15);
+
+    // Equation 30, with the pressure dependent cap that applies below
+    // 70 GHz.
+    let dry = 6.1 * a / (1.0 + 0.17 * rp.powf(-1.1)) * (1.0 + term1 + term2 + term3);
+    let dry = if frequency_ghz < 70.0 {
+        dry.min(10.7 * rp.powf(0.3))
+    } else {
+        dry
+    };
+
+    // Equations 36 through 38 and 35b.
+    let a_wet = 1.929_8 - 0.041_66 * (temperature - 273.15) + 0.051_7 * rho;
+    let b_wet = 1.167_4 - 0.006_22 * (temperature - 273.15) + 0.006_3 * rho;
+    let sigma = 1.013 / (1.0 + (-8.6 * (rp - 0.57)).exp());
+
+    let wet = a_wet
+        + b_wet
+            * WET_HEIGHT_TERMS
+                .iter()
+                .map(|&(line_frequency, coefficient, width)| {
+                    coefficient * sigma / ((frequency_ghz - line_frequency).powi(2) + width * sigma)
+                })
+                .sum::<f64>();
+
+    Some(EquivalentHeights {
+        dry_km: dry,
+        wet_km: wet,
+    })
+}
+
+/// Returns the zenith attenuation of the atmosphere, in decibels.
+///
+/// This is the specific attenuation at the surface multiplied by the
+/// equivalent height, summed over the two constituents.
+///
+/// The frequency is in gigahertz.
+///
+/// Returns `None` when the frequency or atmosphere is not usable.
+pub fn zenith_attenuation(frequency_ghz: f64, atmosphere: Atmosphere) -> Option<f64> {
+    let heights = equivalent_heights(frequency_ghz, atmosphere)?;
+    let dry = specific_attenuation_dry_air(frequency_ghz, atmosphere)?;
+    let wet = specific_attenuation_water_vapour(frequency_ghz, atmosphere)?;
+    Some(dry * heights.dry_km + wet * heights.wet_km)
+}
+
+/// Returns the attenuation along a slant path at an elevation angle, in
+/// decibels.
+///
+/// The zenith attenuation is divided by the sine of the elevation. The
+/// plane parallel approximation this uses is the one the recommendation
+/// applies for the approximate method; it is accurate above about five
+/// degrees of elevation.
+///
+/// The elevation is in radians.
+///
+/// Returns `None` when the frequency or atmosphere is not usable, or when
+/// the elevation is at or below the horizon.
+pub fn slant_path_attenuation(
+    frequency_ghz: f64,
+    elevation: f64,
+    atmosphere: Atmosphere,
+) -> Option<f64> {
+    if !elevation.is_finite() || elevation.sin() <= 0.0 {
+        return None;
+    }
+    Some(zenith_attenuation(frequency_ghz, atmosphere)? / elevation.sin())
+}

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

@@ -4,61 +4,65 @@
 //! mechanisms that are treated separately because they behave differently
 //! with frequency, elevation, and weather.
 //!
-//! Gaseous absorption comes from oxygen and water vapour. Oxygen has
+//! Gaseous absorption comes from the absorption lines of oxygen and water
-//! absorption lines near 60 GHz and 118 GHz and water vapour near 22 GHz, so
+//! vapour. Oxygen has a dense group of lines between 50 and 70 GHz and an
-//! the effect is small below about 10 GHz and grows rapidly above 20 GHz.
+//! isolated line at 118.75 GHz, and water vapour has a strong line at
-//! It depends on the elevation angle through the air mass along the path.
+//! 22.235 GHz and a much stronger one at 183.31 GHz. The attenuation is the
+//! line-by-line sum over those lines, which makes it negligible below about
+//! 10 GHz and severe in the absorption bands.
 //!
-//! Rain attenuation is caused by scattering and absorption in precipitation.
+//! Rain attenuation is caused by scattering and absorption in
-//! It is negligible below about 4 GHz, grows roughly as a power of the rain
+//! precipitation. It is negligible below about 4 GHz, grows as a power of
-//! rate with an exponent near one, and depends strongly on elevation.
+//! the rain rate, and depends strongly on elevation and polarisation.
 //!
 //! Both are functions of the elevation angle through the same path
-//! geometry: a signal at a low elevation crosses more atmosphere. The air
+//! geometry: a signal at a low elevation crosses more atmosphere.
-//! 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
+//! - Gaseous specific attenuation and equivalent heights: ITU-R
-//!   atmospheric gases*.
+//!   Recommendation P.676-12, *Attenuation by atmospheric gases*.
-//! - Rain specific attenuation coefficients: ITU-R Recommendation P.838,
+//! - Rain specific attenuation coefficients: ITU-R Recommendation P.838-3,
 //!   *Specific attenuation model for rain for use in prediction methods*.
-//! - Rain path and elevation dependence: ITU-R Recommendation P.618,
+//! - Rain path and elevation dependence: ITU-R Recommendation P.618-13,
 //!   *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};
+//! use apsis_atmosphere::{Atmosphere, slant_path_attenuation, zenith_attenuation};
 //!
-//! // At the zenith the path is the shortest.
+//! // At the reference atmosphere the 60 GHz oxygen complex is heavily
-//! assert_eq!(air_mass(zenith_to_elevation(0.0)), 1.0);
+//! // attenuated while a Ku band link is barely affected.
+//! let ku = zenith_attenuation(12.0, Atmosphere::REFERENCE).expect("valid");
+//! let oxygen = zenith_attenuation(60.0, Atmosphere::REFERENCE).expect("valid");
+//! assert!(ku < 0.1);
+//! assert!(oxygen > 100.0);
 //!
-//! // At thirty degrees elevation the air mass is twice as large, because
+//! // A lower elevation crosses more atmosphere.
-//! // the cosecant of thirty degrees is two.
+//! let zenith = zenith_attenuation(12.0, Atmosphere::REFERENCE).expect("valid");
-//! let thirty = core::f64::consts::PI / 6.0;
+//! let low = slant_path_attenuation(12.0, 0.2, Atmosphere::REFERENCE).expect("valid");
-//! assert!((air_mass(thirty) - 2.0).abs() < 1e-9);
+//! assert!(low > zenith);
-//!
-//! // 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 gaseous_path;
 mod geometry;
 mod rain;
 
 pub use gaseous::{
-    GaseousAttenuation, gaseous_attenuation, oxygen_zenith_attenuation,
+    Atmosphere, specific_attenuation, specific_attenuation_dry_air,
-    water_vapour_zenith_attenuation,
+    specific_attenuation_water_vapour,
+};
+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 rain::{
-    RainAttenuation, polarisation_tilt, rain_attenuation, rain_coefficients,
+    DEFAULT_RAIN_HEIGHT_KM, RainAttenuation, polarisation_tilt, rain_attenuation,
-    specific_rain_attenuation,
+    rain_coefficients, specific_rain_attenuation,
 };
 
 pub use apsis_rf::Decibel;

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

@@ -1,8 +1,177 @@
 //! Rain attenuation.
+//!
+//! The specific attenuation of rain is modelled as a power law in the rain
+//! rate, `gamma = k R^alpha`, where `k` and `alpha` depend on the frequency,
+//! on the polarisation, and on the elevation. The coefficients are given by
+//! ITU-R Recommendation P.838 as closed form curve fits rather than as a
+//! table, which is what the current revision of that recommendation uses.
+//!
+//! The attenuation along a slant path follows ITU-R Recommendation P.618.
+//! The path length through the rain is limited by the height of the rain
+//! layer above the station, and two reduction factors account for the fact
+//! that the rain rate is not uniform along the path: a horizontal reduction
+//! and a vertical adjustment. Both depend on the specific attenuation, so
+//! the calculation is not a simple product of a rate and a length.
+//!
+//! Sources:
+//!
+//! - Specific attenuation coefficients: ITU-R Recommendation P.838-3,
+//!   *Specific attenuation model for rain for use in prediction methods*,
+//!   equations 2 and 3 and tables 1 and 2.
+//! - Slant path method: ITU-R Recommendation P.618-13, *Propagation data and
+//!   prediction methods required for the design of Earth-space
+//!   telecommunication systems*, equations 1 through 10.
 
-use crate::geometry::air_mass;
 use apsis_rf::Decibel;
-use core::f64::consts::PI;
+
+/// The effective radius of the Earth used by the path geometry, in
+/// kilometers.
+///
+/// The value accounts for the refraction of the ray in the lower atmosphere,
+/// which makes the Earth appear larger than it is.
+const EFFECTIVE_EARTH_RADIUS_KM: f64 = 8_500.0;
+
+/// The default height of the rain layer above sea level, in kilometers.
+///
+/// Recommendation P.618 takes this from the zero degree isotherm height,
+/// which depends on latitude and season and is tabulated separately. This
+/// value is representative of mid latitudes; a caller working in the tropics
+/// should raise it and one in polar regions should lower it.
+pub const DEFAULT_RAIN_HEIGHT_KM: f64 = 3.0;
+
+/// One term of a P.838-3 curve fit.
+#[derive(Debug, Clone, Copy)]
+struct CurveTerm {
+    /// The amplitude coefficient.
+    amplitude: f64,
+    /// The centre of the term, on a base ten logarithmic frequency axis.
+    centre: f64,
+    /// The width of the term, on a base ten logarithmic frequency axis.
+    width: f64,
+}
+
+/// The coefficients of the `k` curve fit for horizontal polarisation.
+const K_HORIZONTAL_TERMS: &[CurveTerm] = &[
+    CurveTerm {
+        amplitude: -5.339_80,
+        centre: -0.100_08,
+        width: 1.130_98,
+    },
+    CurveTerm {
+        amplitude: -0.353_51,
+        centre: 1.269_70,
+        width: 0.454_00,
+    },
+    CurveTerm {
+        amplitude: -0.237_89,
+        centre: 0.860_36,
+        width: 0.153_54,
+    },
+    CurveTerm {
+        amplitude: -0.941_58,
+        centre: 0.645_52,
+        width: 0.168_17,
+    },
+];
+
+/// The coefficients of the `k` curve fit for vertical polarisation.
+const K_VERTICAL_TERMS: &[CurveTerm] = &[
+    CurveTerm {
+        amplitude: -3.805_95,
+        centre: 0.569_34,
+        width: 0.810_61,
+    },
+    CurveTerm {
+        amplitude: -3.449_65,
+        centre: -0.229_11,
+        width: 0.510_59,
+    },
+    CurveTerm {
+        amplitude: -0.399_02,
+        centre: 0.730_42,
+        width: 0.118_99,
+    },
+    CurveTerm {
+        amplitude: 0.501_67,
+        centre: 1.073_19,
+        width: 0.271_95,
+    },
+];
+
+/// The coefficients of the `alpha` curve fit for horizontal polarisation.
+const ALPHA_HORIZONTAL_TERMS: &[CurveTerm] = &[
+    CurveTerm {
+        amplitude: -0.143_18,
+        centre: 1.824_42,
+        width: -0.551_87,
+    },
+    CurveTerm {
+        amplitude: 0.295_91,
+        centre: 0.775_64,
+        width: 0.198_22,
+    },
+    CurveTerm {
+        amplitude: 0.321_77,
+        centre: 0.637_73,
+        width: 0.131_64,
+    },
+    CurveTerm {
+        amplitude: -5.376_10,
+        centre: -0.962_30,
+        width: 1.478_28,
+    },
+    CurveTerm {
+        amplitude: 16.172_1,
+        centre: -3.299_80,
+        width: 3.439_90,
+    },
+];
+
+/// The coefficients of the `alpha` curve fit for vertical polarisation.
+const ALPHA_VERTICAL_TERMS: &[CurveTerm] = &[
+    CurveTerm {
+        amplitude: -0.077_71,
+        centre: 2.338_40,
+        width: -0.762_84,
+    },
+    CurveTerm {
+        amplitude: 0.567_27,
+        centre: 0.955_45,
+        width: 0.540_39,
+    },
+    CurveTerm {
+        amplitude: -0.202_38,
+        centre: 1.145_20,
+        width: 0.268_09,
+    },
+    CurveTerm {
+        amplitude: -48.299_1,
+        centre: 0.791_669,
+        width: 0.116_226,
+    },
+    CurveTerm {
+        amplitude: 48.583_3,
+        centre: 0.791_459,
+        width: 0.116_479,
+    },
+];
+
+/// The linear coefficients of the `k` curve fits.
+const K_SLOPE_HORIZONTAL: f64 = -0.189_61;
+/// The constant of the horizontal `k` curve fit.
+const K_OFFSET_HORIZONTAL: f64 = 0.711_47;
+/// The linear coefficient of the vertical `k` curve fit.
+const K_SLOPE_VERTICAL: f64 = -0.163_98;
+/// The constant of the vertical `k` curve fit.
+const K_OFFSET_VERTICAL: f64 = 0.632_97;
+/// The linear coefficient of the horizontal `alpha` curve fit.
+const ALPHA_SLOPE_HORIZONTAL: f64 = 0.678_49;
+/// The constant of the horizontal `alpha` curve fit.
+const ALPHA_OFFSET_HORIZONTAL: f64 = -1.955_37;
+/// The linear coefficient of the vertical `alpha` curve fit.
+const ALPHA_SLOPE_VERTICAL: f64 = -0.053_739;
+/// The constant of the vertical `alpha` curve fit.
+const ALPHA_OFFSET_VERTICAL: f64 = 0.834_33;
 
 /// A rain attenuation result at a given elevation.
 #[derive(Debug, Clone, Copy, PartialEq)]
@@ -13,6 +182,10 @@ pub struct RainAttenuation {
     pub path_length_km: f64,
     /// The total attenuation along the path, in decibels.
     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.
+    pub vertical_adjustment: f64,
 }
 
 impl RainAttenuation {
@@ -32,65 +205,101 @@ impl RainAttenuation {
     }
 }
 
-/// Returns the coefficients of the specific attenuation power law for a
+/// Evaluates one P.838-3 curve term.
-/// frequency.
 ///
-/// The specific attenuation in rain is modelled as
+/// The term is a Gaussian on a base ten logarithmic frequency axis, so the
+/// frequency must be in gigahertz and positive.
+fn curve_term(term: CurveTerm, frequency_ghz: f64) -> f64 {
+    let logarithm = frequency_ghz.log10();
+    term.amplitude * (-((logarithm - term.centre) / term.width).powi(2)).exp()
+}
+
+/// Returns the `k` for a polarisation, before the elevation correction.
+fn base_k(terms: &[CurveTerm], slope: f64, offset: f64, frequency_ghz: f64) -> f64 {
+    let gaussian: f64 = terms
+        .iter()
+        .map(|term| curve_term(*term, frequency_ghz))
+        .sum();
+    10f64.powf(gaussian + slope * frequency_ghz.log10() + offset)
+}
+
+/// Returns the `alpha` for a polarisation, before the elevation correction.
+fn base_alpha(terms: &[CurveTerm], slope: f64, offset: f64, frequency_ghz: f64) -> f64 {
+    let gaussian: f64 = terms
+        .iter()
+        .map(|term| curve_term(*term, frequency_ghz))
+        .sum();
+    gaussian + slope * frequency_ghz.log10() + offset
+}
+
+/// Returns the P.838-3 coefficients of the specific attenuation power law.
 ///
-/// ```text
+/// The polarisation is an angle in radians between the electric field and
-/// gamma = k R^alpha
+/// the horizontal, so zero is horizontal and a right angle is vertical. The
-/// ```
+/// elevation is in radians.
 ///
-/// where `R` is the rain rate in millimeters per hour and `k` and `alpha`
+/// The result is a pair of the coefficient `k` and the exponent `alpha` of
-/// depend on frequency, on the polarisation, and weakly on the elevation and
+/// the relation `gamma = k R^alpha`, with the rain rate in millimeters per
-/// the rain temperature. The values returned here are for a frequency in
+/// hour and the result in decibels per kilometer.
-/// 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
+/// Returns `None` when the frequency, polarisation, or elevation is not
-/// field and the horizontal, so zero is horizontal, a right angle is
+/// usable. The recommendation is stated for frequencies from 1 GHz to
-/// vertical, and the intermediate values are circular or slant
+/// 1000 GHz.
-/// polarisation. The coefficients at an angle follow from the horizontal and
-/// vertical values through
 ///
-/// ```text
+/// Source: ITU-R Recommendation P.838-3, equations 2 and 3.
-/// 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)) {
+    let frequency_ghz = frequency / 1.0e9;
+    if !frequency_ghz.is_finite() || !(1.0..=1_000.0).contains(&frequency_ghz) {
         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
+    let k_horizontal = base_k(
-        * (horizontal.0
+        K_HORIZONTAL_TERMS,
-            + vertical.0
+        K_SLOPE_HORIZONTAL,
-            + (horizontal.0 - vertical.0) * cos_squared * cos_two_elevation);
+        K_OFFSET_HORIZONTAL,
-    let numerator = horizontal.0 * horizontal.1
+        frequency_ghz,
-        + vertical.0 * vertical.1
+    );
-        + (horizontal.0 * horizontal.1 - vertical.0 * vertical.1) * cos_squared * cos_two_elevation;
+    let k_vertical = base_k(
-    let alpha = if k > 0.0 { numerator / (2.0 * k) } else { 1.0 };
+        K_VERTICAL_TERMS,
+        K_SLOPE_VERTICAL,
+        K_OFFSET_VERTICAL,
+        frequency_ghz,
+    );
+    let alpha_horizontal = base_alpha(
+        ALPHA_HORIZONTAL_TERMS,
+        ALPHA_SLOPE_HORIZONTAL,
+        ALPHA_OFFSET_HORIZONTAL,
+        frequency_ghz,
+    );
+    let alpha_vertical = base_alpha(
+        ALPHA_VERTICAL_TERMS,
+        ALPHA_SLOPE_VERTICAL,
+        ALPHA_OFFSET_VERTICAL,
+        frequency_ghz,
+    );
+
+    // The elevation and polarisation enter through the square of the cosine
+    // of the elevation angle and the cosine of twice the polarisation tilt.
+    // The elevation term vanishes at the zenith, where the distinction
+    // between horizontal and vertical polarisation disappears.
+    let geometry_factor = elevation.cos().powi(2) * (2.0 * polarisation).cos();
+    let k = 0.5 * (k_horizontal + k_vertical + (k_horizontal - k_vertical) * geometry_factor);
+    let alpha = (k_horizontal * alpha_horizontal
+        + k_vertical * alpha_vertical
+        + (k_horizontal * alpha_horizontal - k_vertical * alpha_vertical) * geometry_factor)
+        / (2.0 * k);
     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
+/// The rain rate is in millimeters per hour.
-/// rate is negative, or when the coefficients are not available.
+///
+/// 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,
@@ -106,35 +315,32 @@ pub fn specific_rain_attenuation(
 
 /// Returns the total rain attenuation along a slant path.
 ///
-/// The path is modelled as a rain layer of fixed height, so the geometric
+/// The rain rate is the value exceeded for the percentage of time of
-/// path length through it is the layer height divided by the sine of the
+/// interest, in millimeters per hour. The station height above sea level and
-/// elevation. A reduction factor is applied because the rain rate is not
+/// the rain layer height are in kilometers; pass [`DEFAULT_RAIN_HEIGHT_KM`]
-/// uniform along the path and the longest paths are the ones where the
+/// for the latter when no local value is available.
-/// assumption of uniformity fails most.
 ///
-/// The reduction factor used here is
+/// The method is the one of Recommendation P.618. The slant path length
-///
+/// follows from the rain height, the station height, and the elevation, with
-/// ```text
+/// a correction for the curvature of the Earth. Two reduction factors then
-/// r = 1 / (1 + L_G / L_0)
+/// account for the non-uniformity of the rain: a horizontal factor that
-/// ```
+/// depends on the specific attenuation, the frequency, and the horizontal
-///
+/// projection of the path, and a vertical factor that depends on the
-/// with `L_G` the geometric path and `L_0` a reference length that depends
+/// latitude through the angle at which the path crosses the rain layer.
-/// 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.
+/// when the station is above the rain layer, or when the elevation is at or
+/// below the horizon.
 ///
-/// Source: ITU-R Recommendation P.618, *Propagation data and prediction
+/// Source: ITU-R Recommendation P.618-13, equations 1 through 10.
-/// methods required for the design of Earth-space telecommunication
-/// systems*.
 pub fn rain_attenuation(
     rain_rate: f64,
     frequency: f64,
     polarisation: f64,
     elevation: f64,
+    station_height_km: f64,
+    rain_height_km: f64,
+    latitude: f64,
 ) -> Option<RainAttenuation> {
     if !rain_rate.is_finite() || rain_rate < 0.0 {
         return None;
@@ -142,103 +348,80 @@ pub fn rain_attenuation(
     if !elevation.is_finite() || elevation.sin() <= 0.0 {
         return None;
     }
+    if !station_height_km.is_finite() || !rain_height_km.is_finite() {
+        return None;
+    }
+    if !latitude.is_finite() {
+        return None;
+    }
     let specific = specific_rain_attenuation(rain_rate, frequency, polarisation, elevation)?;
 
-    // The geometric path through the rain layer, and the reduction factor
+    // Equation 2: the slant path length, with the Earth curvature term.
-    // that accounts for its non-uniformity.
+    let height_difference = rain_height_km - station_height_km;
-    const RAIN_HEIGHT_KM: f64 = 3.0;
+    if height_difference <= 0.0 {
-    let geometric_path_km = RAIN_HEIGHT_KM * air_mass(elevation);
+        // The station is at or above the rain layer, so the path is empty.
-    let reference_length_km = 35.0 * (-0.015 * rain_rate).exp();
+        return Some(RainAttenuation {
-    let reduction = 1.0 / (1.0 + geometric_path_km / reference_length_km);
+            specific_db_per_km: specific,
-    let effective_path_km = geometric_path_km * reduction;
+            path_length_km: 0.0,
+            attenuation_db: 0.0,
+            horizontal_reduction: 1.0,
+            vertical_adjustment: 1.0,
+        });
+    }
+    let sine = elevation.sin();
+    let slant_length = 2.0 * height_difference
+        / ((sine * sine + 2.0 * height_difference / EFFECTIVE_EARTH_RADIUS_KM).sqrt() + sine);
 
+    // Equation 3: the horizontal projection of the slant path.
+    let horizontal_length = (slant_length * elevation.cos()).abs();
+
+    // Equation 6: the horizontal reduction factor.
+    let frequency_ghz = frequency / 1.0e9;
+    let horizontal_reduction = 1.0
+        / (1.0 + 0.78 * (horizontal_length * specific / frequency_ghz).sqrt()
+            - 0.38 * (1.0 - (-2.0 * horizontal_length).exp()));
+
+    // Equation 7: the angle at which the path crosses the rain layer.
+    let eta = (height_difference / (horizontal_length * horizontal_reduction)).atan();
+
+    // Equation 8: the longer of the two path estimates.
+    let path_length = if eta > elevation {
+        horizontal_length * horizontal_reduction / elevation.cos()
+    } else {
+        height_difference / sine
+    };
+
+    // Equation 9: the vertical adjustment factor.
+    let latitude_magnitude = latitude.abs();
+    let xi = if latitude_magnitude < 36.0 {
+        36.0 - latitude_magnitude
+    } else {
+        0.0
+    };
+    let vertical_adjustment = 1.0
+        / (1.0
+            + sine.sqrt()
+                * (31.0
+                    * (1.0 - (-(elevation / (1.0 + xi))).exp())
+                    * (path_length * specific).sqrt()
+                    / frequency_ghz.powi(2)
+                    - 0.45));
+
+    // Equation 10: the effective path length and the attenuation.
+    let effective_length = path_length * vertical_adjustment;
     Some(RainAttenuation {
         specific_db_per_km: specific,
-        path_length_km: effective_path_km,
+        path_length_km: effective_length,
-        attenuation_db: specific * effective_path_km,
+        attenuation_db: specific * effective_length,
+        horizontal_reduction,
+        vertical_adjustment,
     })
 }
 
-/// 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
+    degrees_from_horizontal * core::f64::consts::PI / 180.0
 }

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

@@ -1,107 +1,257 @@
 //! Tests for gaseous attenuation.
+//!
+//! The expected values are those of ITU-R Recommendation P.676-12 for the
+//! reference atmosphere: a dry air pressure of 1013.25 hPa, a water vapour
+//! density of 7.5 g/m^3, and a temperature of 288.15 K. They were produced
+//! with a line-by-line implementation of the recommendation and are recorded
+//! here so that a change to the line tables or the formulae shows up as a
+//! failing assertion.
 
 use apsis_atmosphere::{
-    GaseousAttenuation, gaseous_attenuation, oxygen_zenith_attenuation,
+    Atmosphere, equivalent_heights, specific_attenuation, specific_attenuation_dry_air,
-    water_vapour_zenith_attenuation,
+    specific_attenuation_water_vapour,
 };
-use core::f64::consts::{FRAC_PI_2, FRAC_PI_6};
 
-const TOLERANCE: f64 = 1e-12;
+/// The tolerance on a specific attenuation, in decibels per kilometer.
+///
+/// The implementation reproduces the reference to about twelve significant
+/// digits, so this allows for the rounding of the recorded values.
+const TOLERANCE: f64 = 1e-6;
+
+/// The reference atmosphere of the recommendation.
+fn reference() -> Atmosphere {
+    Atmosphere::REFERENCE
+}
 
 #[test]
-fn oxygen_attenuation_is_small_at_low_frequencies() {
+fn the_reference_atmosphere_matches_the_standard() {
-    // Below about 10 GHz the oxygen contribution at the zenith is a small
+    let atmosphere = reference();
-    // fraction of a decibel.
+    assert_eq!(atmosphere.pressure_hpa, 1_013.25);
-    for gigahertz in [1.0, 2.0, 4.0, 6.0, 10.0] {
+    assert_eq!(atmosphere.water_vapour_density, 7.5);
-        let value = oxygen_zenith_attenuation(gigahertz * 1.0e9).expect("valid frequency");
+    assert_eq!(atmosphere.temperature_k, 288.15);
-        assert!(value < 0.2, "{gigahertz} GHz gave {value} dB");
+    assert!(atmosphere.is_valid());
-        assert!(value >= 0.0);
+}
+
+#[test]
+fn the_water_vapour_partial_pressure_follows_the_ideal_gas_law() {
+    let atmosphere = reference();
+    let expected = 7.5 * 288.15 / 216.7;
+    assert!((atmosphere.water_vapour_partial_pressure() - expected).abs() < 1e-12);
+}
+
+#[test]
+fn dry_air_attenuation_matches_the_recommendation() {
+    // The recorded values are from the reference implementation of P.676-12.
+    let cases = [
+        (1.0, 0.005_388_658_168),
+        (10.0, 0.008_224_416_703),
+        (22.235, 0.013_292_678_183),
+        (50.0, 0.277_268_638_531),
+        (60.0, 14.623_474_796_486),
+        (70.0, 0.304_105_097_147),
+        (118.0, 1.134_866_201_871),
+        (300.0, 0.025_759_576_282),
+    ];
+    for (frequency, expected) in cases {
+        let actual = specific_attenuation_dry_air(frequency, reference()).expect("valid");
+        assert!(
+            (actual - expected).abs() < TOLERANCE,
+            "{frequency} GHz: got {actual}, expected {expected}"
+        );
     }
 }
 
 #[test]
-fn oxygen_attenuation_peaks_near_sixty_gigahertz() {
+fn water_vapour_attenuation_matches_the_recommendation() {
-    let peak = oxygen_zenith_attenuation(60.0e9).expect("valid");
+    let cases = [
-    // The 60 GHz complex reaches about 15 dB through the whole atmosphere.
+        (1.0, 0.000_050_904_617),
-    assert!((peak - 15.0).abs() < 0.5, "peak {peak}");
+        (10.0, 0.005_974_125_245),
-    // A frequency on either side is lower.
+        (22.235, 0.178_977_992_373),
-    let below = oxygen_zenith_attenuation(50.0e9).expect("valid");
+        (60.0, 0.154_841_840_636),
-    let above = oxygen_zenith_attenuation(70.0e9).expect("valid");
+        (183.0, 27.665_008_314_166),
-    assert!(below < peak);
+        (300.0, 5.221_329_041_105),
-    assert!(above < peak);
+    ];
-}
+    for (frequency, expected) in cases {
-
+        let actual = specific_attenuation_water_vapour(frequency, reference()).expect("valid");
-#[test]
+        assert!(
-fn oxygen_attenuation_is_monotonic_through_the_complex() {
+            (actual - expected).abs() < TOLERANCE,
-    // The value rises steadily from below the complex into its centre.
+            "{frequency} GHz: got {actual}, expected {expected}"
-    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() {
+fn the_two_contributions_add() {
-    let peak = water_vapour_zenith_attenuation(22.235e9).expect("valid");
+    let frequency = 22.235;
-    assert!((peak - 0.20).abs() < 1e-9, "peak {peak}");
+    let dry = specific_attenuation_dry_air(frequency, reference()).expect("valid");
-    let below = water_vapour_zenith_attenuation(15.0e9).expect("valid");
+    let wet = specific_attenuation_water_vapour(frequency, reference()).expect("valid");
-    let above = water_vapour_zenith_attenuation(30.0e9).expect("valid");
+    let total = specific_attenuation(frequency, reference()).expect("valid");
-    assert!(below < peak);
+    assert!((total - (dry + wet)).abs() < 1e-12);
-    assert!(above < peak);
 }
 
 #[test]
-fn water_vapour_attenuation_is_small_below_ten_gigahertz() {
+fn the_oxygen_complex_dominates_near_sixty_gigahertz() {
-    for gigahertz in [1.0, 4.0, 6.0, 10.0] {
+    // The 60 GHz complex is a set of more than thirty lines, so the dry air
-        let value = water_vapour_zenith_attenuation(gigahertz * 1.0e9).expect("valid");
+    // term is many orders of magnitude above the water vapour term there.
-        assert!(value < 0.05, "{gigahertz} GHz gave {value} dB");
+    let frequency = 60.0;
+    let dry = specific_attenuation_dry_air(frequency, reference()).expect("valid");
+    let wet = specific_attenuation_water_vapour(frequency, reference()).expect("valid");
+    // The dry air value is about 14.6 dB/km against 0.15 dB/km for water
+    // vapour, a ratio near one hundred.
+    assert!(dry > 50.0 * wet, "dry {dry} wet {wet}");
+}
+
+#[test]
+fn the_water_vapour_line_dominates_at_one_hundred_and_eighty_three_gigahertz() {
+    let frequency = 183.0;
+    let dry = specific_attenuation_dry_air(frequency, reference()).expect("valid");
+    let wet = specific_attenuation_water_vapour(frequency, reference()).expect("valid");
+    assert!(wet > dry);
+    // The line reaches over twenty decibels per kilometer.
+    assert!(wet > 20.0, "wet {wet}");
+}
+
+#[test]
+fn attenuation_is_small_in_the_microwave_windows() {
+    // Below 15 GHz the total specific attenuation is a hundredth of a
+    // decibel per kilometer or less, which is why those bands are used.
+    for frequency in [1.0, 2.0, 4.0, 6.0, 10.0, 12.0, 15.0] {
+        let total = specific_attenuation(frequency, reference()).expect("valid");
+        assert!(total < 0.05, "{frequency} GHz gave {total} dB/km");
     }
 }
 
 #[test]
-fn the_total_is_the_sum_of_its_parts() {
+fn dry_air_attenuation_peaks_inside_the_oxygen_complex() {
-    let total = gaseous_attenuation(20.0e9).expect("valid");
+    // The strongest line is near 60.4 GHz, and the value there must exceed
-    assert!((total.zenith_db() - (total.oxygen_db + total.water_vapour_db)).abs() < TOLERANCE);
+    // the values on either side.
+    let peak = specific_attenuation_dry_air(60.434_778, reference()).expect("valid");
+    for frequency in [50.0, 55.0, 65.0, 70.0] {
+        let value = specific_attenuation_dry_air(frequency, reference()).expect("valid");
+        assert!(value < peak, "{frequency} GHz gave {value} against {peak}");
+    }
 }
 
 #[test]
-fn attenuation_grows_as_the_elevation_falls() {
+fn water_vapour_attenuation_grows_with_the_vapour_density() {
-    let gas = GaseousAttenuation::new(0.05, 0.15);
+    let dry_air = Atmosphere::new(1_013.25, 0.0, 288.15);
-    let zenith = gas.at_elevation(FRAC_PI_2);
+    let wet_air = Atmosphere::new(1_013.25, 15.0, 288.15);
-    let thirty = gas.at_elevation(FRAC_PI_6);
+    let dry = specific_attenuation_water_vapour(22.235, dry_air).expect("valid");
-    // The zenith total is 0.2 dB and thirty degrees doubles it.
+    let wet = specific_attenuation_water_vapour(22.235, wet_air).expect("valid");
-    assert!((zenith.db() - 0.2).abs() < TOLERANCE);
+    assert_eq!(dry, 0.0);
-    assert!((thirty.db() - 0.4).abs() < 1e-9);
+    assert!(wet > 0.0);
 }
 
 #[test]
-fn attenuation_at_the_horizon_is_infinite() {
+fn equivalent_heights_match_the_recommendation() {
-    let gas = GaseousAttenuation::new(0.05, 0.15);
+    let cases = [
-    assert!(gas.at_elevation(0.0).db().is_infinite());
+        (1.0, 5.551_003_440, 1.696_760_106),
+        (22.235, 4.867_310_088, 2.358_922_817),
+        (60.0, 10.731_486_113, 1.694_560_689),
+        (118.0, 13.446_411_089, 1.695_050_659),
+        (183.0, 5.581_437_611, 2.522_911_251),
+    ];
+    for (frequency, dry_expected, wet_expected) in cases {
+        let heights = equivalent_heights(frequency, reference()).expect("valid");
+        assert!(
+            (heights.dry_km - dry_expected).abs() < 1e-6,
+            "{frequency} GHz dry height: got {}, expected {dry_expected}",
+            heights.dry_km
+        );
+        assert!(
+            (heights.wet_km - wet_expected).abs() < 1e-6,
+            "{frequency} GHz wet height: got {}, expected {wet_expected}",
+            heights.wet_km
+        );
+    }
 }
 
 #[test]
-fn a_zero_zenith_attenuation_gives_zero_at_every_elevation() {
+fn the_dry_height_is_capped_below_seventy_gigahertz() {
-    let none = GaseousAttenuation::new(0.0, 0.0);
+    // Equation 30 limits the dry air height to a pressure dependent bound at
-    assert_eq!(none.at_elevation(FRAC_PI_6).db(), 0.0);
+    // the lower frequencies. The bound is `10.7 rp^0.3`, where `rp` is the
+    // ratio of the total barometric pressure, dry air plus water vapour, to
+    // the reference pressure.
+    let atmosphere = reference();
+    let total = atmosphere.pressure_hpa + atmosphere.water_vapour_partial_pressure();
+    let rp = total / 1_013.25;
+    let bound = 10.7 * rp.powf(0.3);
+
+    // At 60 GHz the uncapped value is above the bound, so the bound applies.
+    let heights = equivalent_heights(60.0, atmosphere).expect("valid");
+    assert!((heights.dry_km - bound).abs() < 1e-9, "{}", heights.dry_km);
+
+    // Above 70 GHz no cap is applied, and the value exceeds it.
+    let high = equivalent_heights(118.0, atmosphere).expect("valid");
+    assert!(high.dry_km > bound);
 }
 
 #[test]
-fn validity_is_reported() {
+fn zenith_attenuation_matches_the_recommendation() {
-    assert!(GaseousAttenuation::default().is_valid());
+    let cases = [
-    assert!(GaseousAttenuation::new(0.05, 0.15).is_valid());
+        (1.0, 0.029_998_832_95),
-    assert!(!GaseousAttenuation::new(-0.1, 0.15).is_valid());
+        (12.0, 0.058_856_728_55),
-    assert!(!GaseousAttenuation::new(0.05, f64::NAN).is_valid());
+        (22.235, 0.486_894_856_6),
+        (60.0, 157.194_005_6),
+        (183.0, 69.867_434_26),
+    ];
+    for (frequency, expected) in cases {
+        let dry = specific_attenuation_dry_air(frequency, reference()).expect("valid");
+        let wet = specific_attenuation_water_vapour(frequency, reference()).expect("valid");
+        let heights = equivalent_heights(frequency, reference()).expect("valid");
+        let zenith = dry * heights.dry_km + wet * heights.wet_km;
+        assert!(
+            (zenith - expected).abs() < 1e-6,
+            "{frequency} GHz: got {zenith}, expected {expected}"
+        );
+    }
 }
 
 #[test]
-fn invalid_frequencies_are_rejected() {
+fn invalid_inputs_are_rejected() {
-    assert_eq!(oxygen_zenith_attenuation(0.0), None);
+    assert_eq!(specific_attenuation(0.0, reference()), None);
-    assert_eq!(oxygen_zenith_attenuation(-1.0), None);
+    assert_eq!(specific_attenuation(-1.0, reference()), None);
-    assert_eq!(oxygen_zenith_attenuation(f64::NAN), None);
+    assert_eq!(specific_attenuation(f64::NAN, reference()), None);
-    assert_eq!(water_vapour_zenith_attenuation(0.0), None);
+    assert_eq!(
-    assert_eq!(water_vapour_zenith_attenuation(f64::INFINITY), None);
+        specific_attenuation_dry_air(12.0, Atmosphere::new(0.0, 7.5, 288.15)),
-    assert_eq!(gaseous_attenuation(0.0), None);
+        None
+    );
+    assert_eq!(
+        specific_attenuation_water_vapour(12.0, Atmosphere::new(1_013.25, -1.0, 288.15)),
+        None
+    );
+    assert_eq!(equivalent_heights(0.0, reference()), None);
+    assert_eq!(
+        equivalent_heights(12.0, Atmosphere::new(1_013.25, 7.5, 0.0)),
+        None
+    );
+}
+
+#[test]
+fn a_denser_atmosphere_attenuates_more() {
+    let sea_level = Atmosphere::new(1_013.25, 7.5, 288.15);
+    let high_altitude = Atmosphere::new(700.0, 3.0, 280.0);
+    let low = specific_attenuation_dry_air(60.0, sea_level).expect("valid");
+    let high = specific_attenuation_dry_air(60.0, high_altitude).expect("valid");
+    assert!(low > high);
+}
+
+#[test]
+fn slant_path_attenuation_scales_with_the_elevation() {
+    use apsis_atmosphere::{slant_path_attenuation, zenith_attenuation};
+    let frequency = 12.0;
+    let zenith = zenith_attenuation(frequency, reference()).expect("valid");
+    // At thirty degrees elevation the path is twice the zenith path.
+    let thirty =
+        slant_path_attenuation(frequency, core::f64::consts::PI / 6.0, reference()).expect("valid");
+    assert!((thirty - 2.0 * zenith).abs() < 1e-9);
+    assert!(thirty > zenith);
+}
+
+#[test]
+fn slant_path_attenuation_rejects_an_elevation_at_the_horizon() {
+    use apsis_atmosphere::slant_path_attenuation;
+    assert_eq!(slant_path_attenuation(12.0, 0.0, reference()), None);
+    assert_eq!(slant_path_attenuation(12.0, -0.1, reference()), None);
+    assert_eq!(slant_path_attenuation(12.0, f64::NAN, reference()), None);
 }

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

@@ -1,105 +1,152 @@
 //! Tests for rain attenuation.
+//!
+//! The coefficient values are those of ITU-R Recommendation P.838-3,
+//! produced with a reference implementation of the recommendation and
+//! recorded here. The path results exercise the P.618-13 method.
 
 use apsis_atmosphere::{
-    polarisation_tilt, rain_attenuation, rain_coefficients, specific_rain_attenuation,
+    DEFAULT_RAIN_HEIGHT_KM, polarisation_tilt, rain_attenuation, rain_coefficients,
+    specific_rain_attenuation,
 };
-use core::f64::consts::FRAC_PI_2;
+use core::f64::consts::{FRAC_PI_2, PI};
 
 /// Horizontal polarisation.
 const HORIZONTAL: f64 = 0.0;
 /// Vertical polarisation.
 const VERTICAL: f64 = FRAC_PI_2;
-
+/// The zenith, where the two polarisations coincide.
+const ZENITH: f64 = FRAC_PI_2;
 /// A mid latitude elevation of thirty degrees.
-const ELEVATION: f64 = core::f64::consts::PI / 6.0;
+const THIRTY_DEGREES: f64 = PI / 6.0;
 
-#[test]
+/// The relative tolerance on the coefficient values.
-fn coefficients_are_available_across_the_tabulated_range() {
+const TOLERANCE: f64 = 1e-9;
-    for gigahertz in [1.0, 4.0, 12.0, 20.0, 30.0, 100.0, 500.0, 1_000.0] {
+
-        let coefficients =
+fn assert_relative(actual: f64, expected: f64, context: &str) {
-            rain_coefficients(gigahertz * 1.0e9, HORIZONTAL, ELEVATION).expect("valid");
+    let relative = (actual - expected).abs() / expected.abs();
     assert!(
-            coefficients.0 > 0.0,
+        relative < TOLERANCE,
-            "{gigahertz} GHz gave k {}",
+        "{context}: got {actual}, expected {expected}, relative {relative}"
-            coefficients.0
-        );
-        assert!(
-            coefficients.1 > 0.5 && coefficients.1 < 1.5,
-            "{gigahertz} GHz gave alpha {}",
-            coefficients.1
     );
 }
+
+#[test]
+fn coefficients_match_the_recommendation_at_the_zenith() {
+    // At the zenith the elevation term vanishes, so horizontal and vertical
+    // polarisation give the same coefficients.
+    let cases = [
+        (12.0, 0.024_203_061_156_1, 1.151_599_196_3),
+        (20.0, 0.093_876_937_766_6, 1.019_877_631_17),
+        (30.0, 0.234_699_253_968, 0.931_114_875_787),
+    ];
+    for (frequency, k_expected, alpha_expected) in cases {
+        for polarisation in [HORIZONTAL, VERTICAL] {
+            let (k, alpha) =
+                rain_coefficients(frequency * 1.0e9, polarisation, ZENITH).expect("valid");
+            assert_relative(k, k_expected, &format!("{frequency} GHz k"));
+            assert_relative(alpha, alpha_expected, &format!("{frequency} GHz alpha"));
+        }
+    }
 }
 
 #[test]
-fn coefficients_at_a_tabulated_frequency_match_the_table() {
+fn coefficients_match_the_recommendation_at_thirty_degrees() {
-    // At 12 GHz the table gives k = 0.0168 for horizontal polarisation,
+    let cases = [
-    // with alpha_h = 1.200 and alpha_v = 1.000. The combined alpha depends
+        (12.0, 0.023_944_109_789_7, 1.174_670_745_2),
-    // on the elevation through the `cos(2 elevation)` term, so at thirty
+        (20.0, 0.092_201_236_241_3, 1.047_387_560_19),
-    // degrees it is the weighted mean of the two, which is 1.15.
+        (30.0, 0.238_905_952_257, 0.944_198_048_818),
-    let (k, alpha) = rain_coefficients(12.0e9, HORIZONTAL, ELEVATION).expect("valid");
+    ];
-    assert!((k - 0.016_8).abs() < 1e-9, "k {k}");
+    for (frequency, k_expected, alpha_expected) in cases {
-    assert!((alpha - 1.15).abs() < 1e-9, "alpha {alpha}");
+        let (k, alpha) =
-
+            rain_coefficients(frequency * 1.0e9, HORIZONTAL, THIRTY_DEGREES).expect("valid");
-    // At the zenith the `cos(2 elevation)` factor is minus one, which
+        assert_relative(k, k_expected, &format!("{frequency} GHz k"));
-    // selects the vertical alpha instead, so the pair swaps behaviour.
+        assert_relative(alpha, alpha_expected, &format!("{frequency} GHz alpha"));
-    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() {
+fn horizontal_and_vertical_differ_off_the_zenith() {
-    // At 20 GHz the table gives k_h = 0.0751 against k_v = 0.0691.
+    // Away from the zenith the polarisation matters. At 12 GHz and low
-    let (horizontal, _) = rain_coefficients(20.0e9, HORIZONTAL, ELEVATION).expect("valid");
+    // elevation the vertical coefficient is the larger of the two, and the
-    let (vertical, _) = rain_coefficients(20.0e9, VERTICAL, ELEVATION).expect("valid");
+    // two converge as the elevation rises and the geometry becomes
-    assert!(horizontal > vertical, "{horizontal} vs {vertical}");
+    // symmetric.
+    let low_elevation = 10.0_f64.to_radians();
+    let horizontal = rain_coefficients(12.0e9, HORIZONTAL, low_elevation).expect("valid");
+    let vertical = rain_coefficients(12.0e9, VERTICAL, low_elevation).expect("valid");
+    assert_relative(horizontal.0, 0.023_868_203_786_369_1, "horizontal k");
+    assert_relative(vertical.0, 0.024_537_918_525_868_188, "vertical k");
+
+    // At the zenith they coincide.
+    let zenith_horizontal = rain_coefficients(12.0e9, HORIZONTAL, ZENITH).expect("valid");
+    let zenith_vertical = rain_coefficients(12.0e9, VERTICAL, ZENITH).expect("valid");
+    assert_relative(zenith_horizontal.0, zenith_vertical.0, "zenith k");
 }
 
 #[test]
-fn the_coefficients_are_interpolated_between_points() {
+fn circular_polarisation_lies_between_the_linear_cases() {
-    // At 11 GHz the value lies between the 10 and 12 GHz entries in
+    // A tilt of forty five degrees is circular polarisation, whose
-    // logarithmic k, so it is the geometric mean of the two.
+    // coefficients are the mean of the two linear cases.
-    let (k, _) = rain_coefficients(11.0e9, HORIZONTAL, ELEVATION).expect("valid");
+    let horizontal = rain_coefficients(12.0e9, HORIZONTAL, THIRTY_DEGREES).expect("valid");
-    let expected = (0.008_87_f64 * 0.016_8).sqrt();
+    let vertical = rain_coefficients(12.0e9, VERTICAL, THIRTY_DEGREES).expect("valid");
-    assert!((k - expected).abs() < 1e-9, "k {k} expected {expected}");
+    let circular =
+        rain_coefficients(12.0e9, polarisation_tilt(45.0), THIRTY_DEGREES).expect("valid");
+    let mean = 0.5 * (horizontal.0 + vertical.0);
+    assert_relative(circular.0, mean, "circular k");
 }
 
 #[test]
-fn a_frequency_outside_the_table_is_rejected() {
+fn the_coefficients_are_smooth_across_the_band() {
-    assert_eq!(rain_coefficients(0.5e9, HORIZONTAL, ELEVATION), None);
+    // The curve fits are continuous, so a small change in frequency gives a
-    assert_eq!(rain_coefficients(2_000.0e9, HORIZONTAL, ELEVATION), None);
+    // small change in the coefficients.
-    assert_eq!(rain_coefficients(f64::NAN, HORIZONTAL, ELEVATION), None);
+    let mut previous = rain_coefficients(10.0e9, HORIZONTAL, ZENITH).expect("valid");
-    assert_eq!(rain_coefficients(12.0e9, f64::NAN, ELEVATION), None);
+    for tenth in 101..=200 {
+        let frequency = tenth as f64 * 0.1e9;
+        let current = rain_coefficients(frequency, HORIZONTAL, ZENITH).expect("valid");
+        assert!(
+            (current.0 - previous.0).abs() / previous.0 < 0.05,
+            "jump at {frequency} GHz"
+        );
+        previous = current;
+    }
+}
+
+#[test]
+fn a_frequency_outside_the_recommendation_is_rejected() {
+    assert_eq!(rain_coefficients(0.5e9, HORIZONTAL, ZENITH), None);
+    assert_eq!(rain_coefficients(2_000.0e9, HORIZONTAL, ZENITH), None);
+    assert_eq!(rain_coefficients(f64::NAN, HORIZONTAL, ZENITH), None);
+    assert_eq!(rain_coefficients(12.0e9, f64::NAN, ZENITH), None);
+    assert_eq!(rain_coefficients(12.0e9, HORIZONTAL, f64::NAN), None);
 }
 
 #[test]
 fn specific_attenuation_follows_the_power_law() {
     let rain_rate = 20.0;
+    let frequency = 12.0e9;
     let specific =
-        specific_rain_attenuation(rain_rate, 12.0e9, HORIZONTAL, ELEVATION).expect("valid");
+        specific_rain_attenuation(rain_rate, frequency, HORIZONTAL, THIRTY_DEGREES).expect("valid");
-    // The value is k R^alpha with the coefficients for this frequency,
+    let (k, alpha) = rain_coefficients(frequency, HORIZONTAL, THIRTY_DEGREES).expect("valid");
-    // polarisation, and elevation.
+    assert_relative(specific, k * rain_rate.powf(alpha), "specific attenuation");
-    let (k, alpha) = rain_coefficients(12.0e9, HORIZONTAL, ELEVATION).expect("valid");
+    // At 12 GHz, thirty degrees elevation, and 20 mm/h this is 0.808 dB/km.
-    assert!((specific - k * rain_rate.powf(alpha)).abs() < 1e-9);
+    assert!(
-    // For 20 mm/h at 12 GHz and thirty degrees this is about 0.53 dB/km.
+        (specific - 0.808_128_772_7).abs() < 1e-6,
-    assert!((specific - 0.526_6).abs() < 0.001, "specific {specific}");
+        "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 light = specific_rain_attenuation(1.0, 12.0e9, HORIZONTAL, THIRTY_DEGREES).expect("valid");
-    let heavy = specific_rain_attenuation(50.0, 12.0e9, HORIZONTAL, ELEVATION).expect("valid");
+    let heavy = specific_rain_attenuation(50.0, 12.0e9, HORIZONTAL, THIRTY_DEGREES).expect("valid");
     assert!(heavy > light);
-    // The ratio follows the power law, so it is 50 raised to the same alpha
+    let (_, alpha) = rain_coefficients(12.0e9, HORIZONTAL, THIRTY_DEGREES).expect("valid");
-    // the coefficients give.
+    assert_relative(heavy / light, 50.0_f64.powf(alpha), "rate ratio");
-    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),
+        specific_rain_attenuation(0.0, 12.0e9, HORIZONTAL, THIRTY_DEGREES),
         Some(0.0)
     );
 }
@@ -107,40 +154,90 @@ fn a_zero_rain_rate_gives_no_attenuation() {
 #[test]
 fn a_negative_rain_rate_is_rejected() {
     assert_eq!(
-        specific_rain_attenuation(-1.0, 12.0e9, HORIZONTAL, ELEVATION),
+        specific_rain_attenuation(-1.0, 12.0e9, HORIZONTAL, THIRTY_DEGREES),
+        None
+    );
+    assert_eq!(
+        rain_attenuation(
+            -1.0,
+            12.0e9,
+            HORIZONTAL,
+            THIRTY_DEGREES,
+            0.0,
+            DEFAULT_RAIN_HEIGHT_KM,
+            45.0
+        ),
         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");
+    let result = rain_attenuation(
-    assert!(
+        20.0,
-        (result.attenuation_db - result.specific_db_per_km * result.path_length_km).abs() < 1e-9
+        12.0e9,
+        HORIZONTAL,
+        THIRTY_DEGREES,
+        0.0,
+        DEFAULT_RAIN_HEIGHT_KM,
+        45.0,
+    )
+    .expect("valid");
+    assert_relative(
+        result.attenuation_db,
+        result.specific_db_per_km * result.path_length_km,
+        "total attenuation",
     );
     assert!(result.is_valid());
 }
 
 #[test]
-fn the_effective_path_is_shorter_than_the_geometric_one() {
+fn the_reduction_factors_are_physical() {
-    // The reduction factor represents the fact that heavy rain does not fill
+    // The horizontal reduction factor accounts for rain not filling the
-    // the whole slant path, so the effective length is below the geometric
+    // whole horizontal extent of the path, so it is at most one. The
-    // length through a three kilometre layer.
+    // vertical adjustment is the ratio of two path estimates and can exceed
-    let result = rain_attenuation(20.0, 12.0e9, HORIZONTAL, ELEVATION).expect("valid");
+    // one when the path crosses the rain layer at a shallow angle, which is
-    let geometric = 3.0 / ELEVATION.sin();
+    // the case Recommendation P.618 allows for.
-    assert!((geometric - 6.0).abs() < 1e-6);
+    let result = rain_attenuation(
+        20.0,
+        12.0e9,
+        HORIZONTAL,
+        THIRTY_DEGREES,
+        0.0,
+        DEFAULT_RAIN_HEIGHT_KM,
+        45.0,
+    )
+    .expect("valid");
     assert!(
-        result.path_length_km < geometric,
+        result.horizontal_reduction > 0.0 && result.horizontal_reduction <= 1.0,
-        "{}",
+        "horizontal reduction {}",
-        result.path_length_km
+        result.horizontal_reduction
     );
-    assert!(result.path_length_km > 0.0);
+    assert!(result.vertical_adjustment > 0.0, "vertical adjustment");
 }
 
 #[test]
 fn attenuation_grows_as_the_elevation_falls() {
-    let high = rain_attenuation(20.0, 12.0e9, HORIZONTAL, FRAC_PI_2).expect("valid");
+    let high = rain_attenuation(
-    let low = rain_attenuation(20.0, 12.0e9, HORIZONTAL, 10.0_f64.to_radians()).expect("valid");
+        20.0,
+        12.0e9,
+        HORIZONTAL,
+        ZENITH,
+        0.0,
+        DEFAULT_RAIN_HEIGHT_KM,
+        45.0,
+    )
+    .expect("valid");
+    let low = rain_attenuation(
+        20.0,
+        12.0e9,
+        HORIZONTAL,
+        10.0_f64.to_radians(),
+        0.0,
+        DEFAULT_RAIN_HEIGHT_KM,
+        45.0,
+    )
+    .expect("valid");
     assert!(low.attenuation_db > high.attenuation_db);
 }
 
@@ -148,9 +245,18 @@ fn attenuation_grows_as_the_elevation_falls() {
 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");
+    let result = rain_attenuation(
+        20.0,
+        12.0e9,
+        HORIZONTAL,
+        THIRTY_DEGREES,
+        0.0,
+        DEFAULT_RAIN_HEIGHT_KM,
+        45.0,
+    )
+    .expect("valid");
     assert!(
-        result.attenuation_db > 0.5 && result.attenuation_db < 5.0,
+        result.attenuation_db > 1.0 && result.attenuation_db < 6.0,
         "{} dB",
         result.attenuation_db
     );
@@ -158,20 +264,132 @@ fn a_heavy_rain_event_gives_a_significant_attenuation() {
 
 #[test]
 fn a_lower_frequency_suffers_less_rain_attenuation() {
-    let c_band = rain_attenuation(20.0, 6.0e9, HORIZONTAL, ELEVATION).expect("valid");
+    let c_band = rain_attenuation(
-    let ku_band = rain_attenuation(20.0, 12.0e9, HORIZONTAL, ELEVATION).expect("valid");
+        20.0,
+        6.0e9,
+        HORIZONTAL,
+        THIRTY_DEGREES,
+        0.0,
+        DEFAULT_RAIN_HEIGHT_KM,
+        45.0,
+    )
+    .expect("valid");
+    let ku_band = rain_attenuation(
+        20.0,
+        12.0e9,
+        HORIZONTAL,
+        THIRTY_DEGREES,
+        0.0,
+        DEFAULT_RAIN_HEIGHT_KM,
+        45.0,
+    )
+    .expect("valid");
     assert!(ku_band.attenuation_db > c_band.attenuation_db);
 }
 
+#[test]
+fn a_station_above_the_rain_layer_has_no_rain_path() {
+    let result = rain_attenuation(
+        20.0,
+        12.0e9,
+        HORIZONTAL,
+        THIRTY_DEGREES,
+        4.0,
+        DEFAULT_RAIN_HEIGHT_KM,
+        45.0,
+    )
+    .expect("valid");
+    assert_eq!(result.attenuation_db, 0.0);
+    assert_eq!(result.path_length_km, 0.0);
+}
+
 #[test]
 fn an_elevation_at_the_horizon_is_rejected() {
-    assert_eq!(rain_attenuation(20.0, 12.0e9, HORIZONTAL, 0.0), None);
+    assert_eq!(
-    assert_eq!(rain_attenuation(20.0, 12.0e9, HORIZONTAL, -0.1), None);
+        rain_attenuation(
+            20.0,
+            12.0e9,
+            HORIZONTAL,
+            0.0,
+            0.0,
+            DEFAULT_RAIN_HEIGHT_KM,
+            45.0
+        ),
+        None
+    );
+    assert_eq!(
+        rain_attenuation(
+            20.0,
+            12.0e9,
+            HORIZONTAL,
+            -0.1,
+            0.0,
+            DEFAULT_RAIN_HEIGHT_KM,
+            45.0
+        ),
+        None
+    );
+}
+
+#[test]
+fn the_latitude_enters_through_the_vertical_adjustment() {
+    // Below thirty six degrees the latitude term is nonzero, so the same
+    // geometry gives a different result at two latitudes.
+    let temperate = rain_attenuation(
+        20.0,
+        12.0e9,
+        HORIZONTAL,
+        THIRTY_DEGREES,
+        0.0,
+        DEFAULT_RAIN_HEIGHT_KM,
+        45.0,
+    )
+    .expect("valid");
+    let tropical = rain_attenuation(
+        20.0,
+        12.0e9,
+        HORIZONTAL,
+        THIRTY_DEGREES,
+        0.0,
+        DEFAULT_RAIN_HEIGHT_KM,
+        10.0,
+    )
+    .expect("valid");
+    assert!(
+        (temperate.vertical_adjustment - tropical.vertical_adjustment).abs() > 1e-6,
+        "{} vs {}",
+        temperate.vertical_adjustment,
+        tropical.vertical_adjustment
+    );
+}
+
+#[test]
+fn non_finite_geometry_is_rejected() {
+    assert_eq!(
+        rain_attenuation(
+            20.0,
+            12.0e9,
+            HORIZONTAL,
+            THIRTY_DEGREES,
+            f64::NAN,
+            3.0,
+            45.0
+        ),
+        None
+    );
+    assert_eq!(
+        rain_attenuation(20.0, 12.0e9, HORIZONTAL, THIRTY_DEGREES, 0.0, 3.0, f64::NAN),
+        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_relative(polarisation_tilt(90.0), FRAC_PI_2, "ninety degrees");
-    assert!((polarisation_tilt(45.0) - FRAC_PI_2 / 2.0).abs() < 1e-12);
+    assert_relative(
+        polarisation_tilt(45.0),
+        FRAC_PI_2 / 2.0,
+        "forty five degrees",
+    );
 }

validation/pixi.lock

@@ -48,6 +48,8 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/noarch/ca-certificates-2026.7.22-hbd8a1cb_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/tzdata-2026c-h151e31d_0.conda
       - pypi: https://files.pythonhosted.org/packages/0b/a7/71ac2cff56fec219ed242bb11b8efb69fcc4bec75db06fb7bfe35de520e6/certifi-2026.7.22-py3-none-any.whl
+      - pypi: https://files.pythonhosted.org/packages/5f/ed/035354de8fbc1c1350f5b7f66361995f187d89675622eb63016bf7261dd5/pyproj-3.8.0-cp312-cp312-manylinux_2_28_x86_64.whl
+      - pypi: https://files.pythonhosted.org/packages/61/7b/e682678c0a6fcdd4529abc5f22324149cd7a725138465d76885a3a53c88f/itur-0.4.0-py2.py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/63/34/ba1c580383c9eada3711951fef0795c80b829a078d72188184bcab9dd527/packaging-26.3-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/65/af/aa78d1a88805456e212b65461354cd943197fb9acecc4c90fd12295123a3/numpy-2.5.3-cp312-cp312-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl
       - pypi: https://files.pythonhosted.org/packages/7f/13/a24edfb7480eb3f2d01d2181dad934f9761d9bca2c256c75ac2753a22890/astropy_iers_data-0.2026.9.21.0.56.25-py3-none-any.whl
@@ -55,6 +57,7 @@ environments:
       - pypi: https://files.pythonhosted.org/packages/8d/5a/c4cc09d4421de255ca0504485efceb3646a4e5381bfb665647e39bea5cb1/skyfield-1.55-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/9d/3f/b9d5739352badc11ca637c8f72525d519458622936bc3313ddefdc7dee96/jplephem-2.24-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/d6/75/633ad7874db149df164168090951e676cdc3541a12d43b2cbb9d20ecb721/sgp4-2.27-cp311-abi3-manylinux1_x86_64.manylinux_2_28_x86_64.manylinux_2_5_x86_64.whl
+      - pypi: https://files.pythonhosted.org/packages/df/64/ff35eb9e54894cf471ff4716abd3c81eb0a0626869217ce3e6ba4ccf17d7/scipy-1.18.1-cp312-cp312-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl
       - pypi: https://files.pythonhosted.org/packages/e5/e0/050018d855d26d3c0b4a7d1b2ed692be758ce276d8289e2a2b44ba1014a5/pyerfa-2.0.1.5-cp39-abi3-manylinux_2_17_x86_64.manylinux2014_x86_64.whl
       - pypi: https://files.pythonhosted.org/packages/fa/14/6f4427419f8a02e1d0a885fb64d8d4454c5d443627aa7cb1a3a36e16abb7/astropy-8.0.1-cp311-abi3-manylinux2014_x86_64.manylinux_2_17_x86_64.manylinux_2_28_x86_64.whl
       osx-arm64:
@@ -74,16 +77,19 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/python-3.12.14-hd05a0c4_3_cpython.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/readline-8.3-h8b90a29_1.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/tk-8.6.13-hbeba79b_4.conda
+      - pypi: https://files.pythonhosted.org/packages/05/4a/78c6285577c375e7cf27277ea8ee6961224327f1e1a0c44af5f17f23635c/scipy-1.18.1-cp312-cp312-macosx_12_0_arm64.whl
       - pypi: https://files.pythonhosted.org/packages/0b/a7/71ac2cff56fec219ed242bb11b8efb69fcc4bec75db06fb7bfe35de520e6/certifi-2026.7.22-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/11/4a/31a363370478b63c6289a34743f2ba2d3ae1bd8223e004d18ab28fb92385/pyerfa-2.0.1.5-cp39-abi3-macosx_11_0_arm64.whl
       - pypi: https://files.pythonhosted.org/packages/36/d9/05fd166d8f74b326590e78bb13148d6b359ae32be5cd39082fd8a292bf50/sgp4-2.27-cp311-abi3-macosx_11_0_arm64.whl
       - pypi: https://files.pythonhosted.org/packages/60/39/789131c1188c078dcb3a1692e72e1e050c68b88ffe72c9ccaac9bcd7a9cd/numpy-2.5.3-cp312-cp312-macosx_11_0_arm64.whl
+      - pypi: https://files.pythonhosted.org/packages/61/7b/e682678c0a6fcdd4529abc5f22324149cd7a725138465d76885a3a53c88f/itur-0.4.0-py2.py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/63/34/ba1c580383c9eada3711951fef0795c80b829a078d72188184bcab9dd527/packaging-26.3-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/7f/13/a24edfb7480eb3f2d01d2181dad934f9761d9bca2c256c75ac2753a22890/astropy_iers_data-0.2026.9.21.0.56.25-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/89/a0/6cf41a19a1f2f3feab0e9c0b74134aa2ce6849093d5517a0c550fe37a648/pyyaml-6.0.3-cp312-cp312-macosx_11_0_arm64.whl
       - pypi: https://files.pythonhosted.org/packages/8d/5a/c4cc09d4421de255ca0504485efceb3646a4e5381bfb665647e39bea5cb1/skyfield-1.55-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/9d/3f/b9d5739352badc11ca637c8f72525d519458622936bc3313ddefdc7dee96/jplephem-2.24-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/a6/34/074f367d5699a1008b24c50acc1539f05f2cfc881b1901f4b110d57eae20/astropy-8.0.1-cp311-abi3-macosx_11_0_arm64.whl
+      - pypi: https://files.pythonhosted.org/packages/c8/29/6598570c90cbfc84ddefc3ccac4aa412bf51a527d72c74cc4fe64a5e6f24/pyproj-3.8.0.tar.gz
       win-64:
       - conda: https://conda.anaconda.org/conda-forge/noarch/ca-certificates-2026.7.22-h4c7d964_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/tzdata-2026c-h151e31d_0.conda
@@ -102,13 +108,16 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/win-64/vc14_runtime-14.51.36247-habf1de7_41.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/vcomp14-14.51.36247-habf1de7_41.conda
       - pypi: https://files.pythonhosted.org/packages/0b/a7/71ac2cff56fec219ed242bb11b8efb69fcc4bec75db06fb7bfe35de520e6/certifi-2026.7.22-py3-none-any.whl
+      - pypi: https://files.pythonhosted.org/packages/39/e7/979fd14e75008623df31ba70d6bb144700f68feadcea042021c06a05bf82/scipy-1.18.1-cp312-cp312-win_amd64.whl
       - pypi: https://files.pythonhosted.org/packages/3c/a1/accf6d4f0c80c5d9ba9735d6b1550e444180599f34dec69ca01360f717ad/numpy-2.5.3-cp312-cp312-win_amd64.whl
+      - pypi: https://files.pythonhosted.org/packages/61/7b/e682678c0a6fcdd4529abc5f22324149cd7a725138465d76885a3a53c88f/itur-0.4.0-py2.py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/63/34/ba1c580383c9eada3711951fef0795c80b829a078d72188184bcab9dd527/packaging-26.3-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/7f/13/a24edfb7480eb3f2d01d2181dad934f9761d9bca2c256c75ac2753a22890/astropy_iers_data-0.2026.9.21.0.56.25-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/86/bf/899e81e4cce32febab4fb42bb97dcdf66bc135272882d1987881a4b519e9/pyyaml-6.0.3-cp312-cp312-win_amd64.whl
       - pypi: https://files.pythonhosted.org/packages/8d/5a/c4cc09d4421de255ca0504485efceb3646a4e5381bfb665647e39bea5cb1/skyfield-1.55-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/9d/3f/b9d5739352badc11ca637c8f72525d519458622936bc3313ddefdc7dee96/jplephem-2.24-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/b4/11/97233cf23ad5411ac6f13b1d6ee3888f90ace4f974d9bf9db887aa428912/pyerfa-2.0.1.5-cp39-abi3-win_amd64.whl
+      - pypi: https://files.pythonhosted.org/packages/b7/d7/2c861a429581577a441e03020d27274812d7969748f5951d972f6e34597f/pyproj-3.8.0-cp312-cp312-win_amd64.whl
       - pypi: https://files.pythonhosted.org/packages/d5/b5/6e76ddece37a130c1c4c33a0ca2ae381c7aad733ca8e3d9b451594711fc0/sgp4-2.27-cp311-abi3-win_amd64.whl
       - pypi: https://files.pythonhosted.org/packages/ec/6a/8534fd98e2d3eb51cbaa261e6c6b59e4ad969a1787d72bb1208499080dbd/astropy-8.0.1-cp311-abi3-win_amd64.whl
 packages:
@@ -894,6 +903,52 @@ packages:
     - vcomp14 >=14.51.36247
   size: 155910
   timestamp: 1785359349999
+- pypi: https://files.pythonhosted.org/packages/05/4a/78c6285577c375e7cf27277ea8ee6961224327f1e1a0c44af5f17f23635c/scipy-1.18.1-cp312-cp312-macosx_12_0_arm64.whl
+  name: scipy
+  version: 1.18.1
+  sha256: e708533e8b2ae2497d65346538a7dcc92814410b25b81432eac66de0f2af8265
+  requires_dist:
+  - numpy>=2.0.0,<2.8
+  - pytest>=8.0.0 ; extra == 'test'
+  - pytest-cov ; extra == 'test'
+  - pytest-timeout ; extra == 'test'
+  - pytest-xdist ; extra == 'test'
+  - asv ; extra == 'test'
+  - mpmath ; extra == 'test'
+  - gmpy2 ; extra == 'test'
+  - threadpoolctl ; extra == 'test'
+  - scikit-umfpack ; extra == 'test'
+  - pooch ; extra == 'test'
+  - hypothesis>=6.30 ; extra == 'test'
+  - array-api-strict>=2.3.1 ; extra == 'test'
+  - cython ; extra == 'test'
+  - meson ; extra == 'test'
+  - ninja ; sys_platform != 'emscripten' and extra == 'test'
+  - scipy-doctest>=2.0.0 ; extra == 'test'
+  - sphinx>=5.0.0,<8.2.0 ; extra == 'doc'
+  - intersphinx-registry ; extra == 'doc'
+  - pydata-sphinx-theme>=0.15.2 ; extra == 'doc'
+  - sphinx-copybutton ; extra == 'doc'
+  - sphinx-design>=0.4.0 ; extra == 'doc'
+  - matplotlib>=3.5 ; extra == 'doc'
+  - numpydoc ; extra == 'doc'
+  - jupytext ; extra == 'doc'
+  - myst-nb>=1.2.0 ; extra == 'doc'
+  - pooch ; extra == 'doc'
+  - jupyterlite-sphinx>=0.19.1 ; extra == 'doc'
+  - jupyterlite-pyodide-kernel ; extra == 'doc'
+  - linkify-it-py ; extra == 'doc'
+  - tabulate ; extra == 'doc'
+  - click<8.3.0 ; extra == 'dev'
+  - spin ; extra == 'dev'
+  - mypy==1.19.1 ; extra == 'dev'
+  - pyrefly==0.63.0 ; extra == 'dev'
+  - typing-extensions ; extra == 'dev'
+  - types-psutil ; extra == 'dev'
+  - pycodestyle ; extra == 'dev'
+  - ruff>=0.12.0 ; extra == 'dev'
+  - cython-lint>=0.12.2 ; extra == 'dev'
+  requires_python: '>=3.12'
 - pypi: https://files.pythonhosted.org/packages/0b/a7/71ac2cff56fec219ed242bb11b8efb69fcc4bec75db06fb7bfe35de520e6/certifi-2026.7.22-py3-none-any.whl
   name: certifi
   version: 2026.7.22
@@ -914,16 +969,78 @@ packages:
   version: '2.27'
   sha256: ee68f73cb2383bc950eada2ca1b10357b0165dbf5ab75185962009375bab8761
   requires_python: '>=2.7,!=3.0.*,!=3.1.*,!=3.2.*,!=3.3.*,!=3.4.*,!=3.5.*,!=3.6.*,!=3.7.*,!=3.8.*,!=3.9.*'
+- pypi: https://files.pythonhosted.org/packages/39/e7/979fd14e75008623df31ba70d6bb144700f68feadcea042021c06a05bf82/scipy-1.18.1-cp312-cp312-win_amd64.whl
+  name: scipy
+  version: 1.18.1
+  sha256: 5e4d44984abc0020154ea81b247adeddcc3ac5527b975ff798bd1ba0adc513c2
+  requires_dist:
+  - numpy>=2.0.0,<2.8
+  - pytest>=8.0.0 ; extra == 'test'
+  - pytest-cov ; extra == 'test'
+  - pytest-timeout ; extra == 'test'
+  - pytest-xdist ; extra == 'test'
+  - asv ; extra == 'test'
+  - mpmath ; extra == 'test'
+  - gmpy2 ; extra == 'test'
+  - threadpoolctl ; extra == 'test'
+  - scikit-umfpack ; extra == 'test'
+  - pooch ; extra == 'test'
+  - hypothesis>=6.30 ; extra == 'test'
+  - array-api-strict>=2.3.1 ; extra == 'test'
+  - cython ; extra == 'test'
+  - meson ; extra == 'test'
+  - ninja ; sys_platform != 'emscripten' and extra == 'test'
+  - scipy-doctest>=2.0.0 ; extra == 'test'
+  - sphinx>=5.0.0,<8.2.0 ; extra == 'doc'
+  - intersphinx-registry ; extra == 'doc'
+  - pydata-sphinx-theme>=0.15.2 ; extra == 'doc'
+  - sphinx-copybutton ; extra == 'doc'
+  - sphinx-design>=0.4.0 ; extra == 'doc'
+  - matplotlib>=3.5 ; extra == 'doc'
+  - numpydoc ; extra == 'doc'
+  - jupytext ; extra == 'doc'
+  - myst-nb>=1.2.0 ; extra == 'doc'
+  - pooch ; extra == 'doc'
+  - jupyterlite-sphinx>=0.19.1 ; extra == 'doc'
+  - jupyterlite-pyodide-kernel ; extra == 'doc'
+  - linkify-it-py ; extra == 'doc'
+  - tabulate ; extra == 'doc'
+  - click<8.3.0 ; extra == 'dev'
+  - spin ; extra == 'dev'
+  - mypy==1.19.1 ; extra == 'dev'
+  - pyrefly==0.63.0 ; extra == 'dev'
+  - typing-extensions ; extra == 'dev'
+  - types-psutil ; extra == 'dev'
+  - pycodestyle ; extra == 'dev'
+  - ruff>=0.12.0 ; extra == 'dev'
+  - cython-lint>=0.12.2 ; extra == 'dev'
+  requires_python: '>=3.12'
 - pypi: https://files.pythonhosted.org/packages/3c/a1/accf6d4f0c80c5d9ba9735d6b1550e444180599f34dec69ca01360f717ad/numpy-2.5.3-cp312-cp312-win_amd64.whl
   name: numpy
   version: 2.5.3
   sha256: 0a59a421a32580a009e8a1751345bf829631b990dc1794b80514ab722b435def
   requires_python: '>=3.12'
+- pypi: https://files.pythonhosted.org/packages/5f/ed/035354de8fbc1c1350f5b7f66361995f187d89675622eb63016bf7261dd5/pyproj-3.8.0-cp312-cp312-manylinux_2_28_x86_64.whl
+  name: pyproj
+  version: 3.8.0
+  sha256: 6b34cec6bdd66b721980c003cc2a0fce0cab7949444ad037fc2a7ecb6f3b996c
+  requires_dist:
+  - certifi
+  requires_python: '>=3.12'
 - pypi: https://files.pythonhosted.org/packages/60/39/789131c1188c078dcb3a1692e72e1e050c68b88ffe72c9ccaac9bcd7a9cd/numpy-2.5.3-cp312-cp312-macosx_11_0_arm64.whl
   name: numpy
   version: 2.5.3
   sha256: f59a878c33d6b88122d80d239bb3b845d58708750b0cb06a09aebb9b18ec696c
   requires_python: '>=3.12'
+- pypi: https://files.pythonhosted.org/packages/61/7b/e682678c0a6fcdd4529abc5f22324149cd7a725138465d76885a3a53c88f/itur-0.4.0-py2.py3-none-any.whl
+  name: itur
+  version: 0.4.0
+  sha256: d7a357172216075b9f0b8f38cd68ce975dba1b7e22db8329f013e6ef651db9b2
+  requires_dist:
+  - astropy
+  - scipy
+  - numpy
+  - pyproj
 - pypi: https://files.pythonhosted.org/packages/63/34/ba1c580383c9eada3711951fef0795c80b829a078d72188184bcab9dd527/packaging-26.3-py3-none-any.whl
   name: packaging
   version: '26.3'
@@ -1054,6 +1171,20 @@ packages:
   - pytest ; extra == 'test'
   - pytest-doctestplus>=0.7 ; extra == 'test'
   requires_python: '>=3.9'
+- pypi: https://files.pythonhosted.org/packages/b7/d7/2c861a429581577a441e03020d27274812d7969748f5951d972f6e34597f/pyproj-3.8.0-cp312-cp312-win_amd64.whl
+  name: pyproj
+  version: 3.8.0
+  sha256: 0556ce011e1530aea2084a12e39dc9c6ccd25a16c1db5457f39bf83161e7c301
+  requires_dist:
+  - certifi
+  requires_python: '>=3.12'
+- pypi: https://files.pythonhosted.org/packages/c8/29/6598570c90cbfc84ddefc3ccac4aa412bf51a527d72c74cc4fe64a5e6f24/pyproj-3.8.0.tar.gz
+  name: pyproj
+  version: 3.8.0
+  sha256: efa59725bba68bf97fa808b61302df32934acdceb6a5c92a8dd0e71dc266a876
+  requires_dist:
+  - certifi
+  requires_python: '>=3.12'
 - pypi: https://files.pythonhosted.org/packages/d5/b5/6e76ddece37a130c1c4c33a0ca2ae381c7aad733ca8e3d9b451594711fc0/sgp4-2.27-cp311-abi3-win_amd64.whl
   name: sgp4
   version: '2.27'
@@ -1064,6 +1195,52 @@ packages:
   version: '2.27'
   sha256: 4d3775313120dcb6239535fa0c94cb6d5b089fd84bbfb15220923ae38dec7296
   requires_python: '>=2.7,!=3.0.*,!=3.1.*,!=3.2.*,!=3.3.*,!=3.4.*,!=3.5.*,!=3.6.*,!=3.7.*,!=3.8.*,!=3.9.*'
+- pypi: https://files.pythonhosted.org/packages/df/64/ff35eb9e54894cf471ff4716abd3c81eb0a0626869217ce3e6ba4ccf17d7/scipy-1.18.1-cp312-cp312-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl
+  name: scipy
+  version: 1.18.1
+  sha256: f55fa87b6c612ecd6b058f167c53231b1d14e412efe361d3d6e38b3631c73218
+  requires_dist:
+  - numpy>=2.0.0,<2.8
+  - pytest>=8.0.0 ; extra == 'test'
+  - pytest-cov ; extra == 'test'
+  - pytest-timeout ; extra == 'test'
+  - pytest-xdist ; extra == 'test'
+  - asv ; extra == 'test'
+  - mpmath ; extra == 'test'
+  - gmpy2 ; extra == 'test'
+  - threadpoolctl ; extra == 'test'
+  - scikit-umfpack ; extra == 'test'
+  - pooch ; extra == 'test'
+  - hypothesis>=6.30 ; extra == 'test'
+  - array-api-strict>=2.3.1 ; extra == 'test'
+  - cython ; extra == 'test'
+  - meson ; extra == 'test'
+  - ninja ; sys_platform != 'emscripten' and extra == 'test'
+  - scipy-doctest>=2.0.0 ; extra == 'test'
+  - sphinx>=5.0.0,<8.2.0 ; extra == 'doc'
+  - intersphinx-registry ; extra == 'doc'
+  - pydata-sphinx-theme>=0.15.2 ; extra == 'doc'
+  - sphinx-copybutton ; extra == 'doc'
+  - sphinx-design>=0.4.0 ; extra == 'doc'
+  - matplotlib>=3.5 ; extra == 'doc'
+  - numpydoc ; extra == 'doc'
+  - jupytext ; extra == 'doc'
+  - myst-nb>=1.2.0 ; extra == 'doc'
+  - pooch ; extra == 'doc'
+  - jupyterlite-sphinx>=0.19.1 ; extra == 'doc'
+  - jupyterlite-pyodide-kernel ; extra == 'doc'
+  - linkify-it-py ; extra == 'doc'
+  - tabulate ; extra == 'doc'
+  - click<8.3.0 ; extra == 'dev'
+  - spin ; extra == 'dev'
+  - mypy==1.19.1 ; extra == 'dev'
+  - pyrefly==0.63.0 ; extra == 'dev'
+  - typing-extensions ; extra == 'dev'
+  - types-psutil ; extra == 'dev'
+  - pycodestyle ; extra == 'dev'
+  - ruff>=0.12.0 ; extra == 'dev'
+  - cython-lint>=0.12.2 ; extra == 'dev'
+  requires_python: '>=3.12'
 - pypi: https://files.pythonhosted.org/packages/e5/e0/050018d855d26d3c0b4a7d1b2ed692be758ce276d8289e2a2b44ba1014a5/pyerfa-2.0.1.5-cp39-abi3-manylinux_2_17_x86_64.manylinux2014_x86_64.whl
   name: pyerfa
   version: 2.0.1.5

validation/pixi.toml

@@ -15,5 +15,6 @@ numpy = "*"
 sgp4 = "*"
 astropy = "*"
 skyfield = "*"
+itur = ">=0.4.0, <0.5"
 
 [tasks]

validation/python/atmosphere/check_gas.py

@@ -0,0 +1,136 @@
+"""Cross-check the ApsisRS gaseous attenuation against ITU-R P.676-12.
+
+Reads the values written by ``cargo run -p apsis-rs --features rf --example
+dump_gas`` and recomputes each one with the ``itur`` package, which
+implements the ITU-R recommendations directly, then reports the largest
+difference in each column.
+
+The ``itur`` package is a reference implementation of the recommendations
+maintained as part of the ITU-R study programme, so it is an independent
+check on the line tables and the formulae rather than a restatement of them.
+
+The script needs no network access.
+
+Usage:
+
+    pixi run python validation/python/atmosphere/check_gas.py [path]
+
+Exit status is zero when every value agrees within tolerance and one
+otherwise.
+"""
+
+from __future__ import annotations
+
+import csv
+import sys
+from pathlib import Path
+
+import numpy as np
+import astropy.units as u
+
+from itur.models.itu676 import _ITU676_12_
+
+# The reference atmosphere of the recommendation, matching the Rust side.
+PRESSURE_HPA = 1013.25
+WATER_VAPOUR_DENSITY = 7.5
+TEMPERATURE_K = 288.15
+
+# Tolerances.
+#
+# The specific attenuation spans several orders of magnitude, and the water
+# vapour term at low frequencies is a small difference between larger
+# contributions, so a pure relative test is too strict there. The check
+# passes when either the relative difference is small or the absolute
+# difference is negligible.
+RELATIVE_TOLERANCE = 1e-8
+ABSOLUTE_TOLERANCE = 1e-10
+
+DEFAULT_INPUT = Path("validation/data/expected/rust_gas.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("frequency_ghz")]
+
+
+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_gas > " + 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 = _ITU676_12_()
+
+    columns = ["dry_db_per_km", "wet_db_per_km", "dry_height_km", "wet_height_km", "zenith_db"]
+    worst = {column: 0.0 for column in columns}
+    worst_absolute = {column: 0.0 for column in columns}
+    worst_at = {column: None for column in columns}
+
+    for row in rows:
+        frequency = float(row["frequency_ghz"])
+        dry = float(model.gamma0_exact(frequency, PRESSURE_HPA, WATER_VAPOUR_DENSITY, TEMPERATURE_K))
+        wet = float(model.gammaw_exact(frequency, PRESSURE_HPA, WATER_VAPOUR_DENSITY, TEMPERATURE_K))
+        heights = model.slant_inclined_path_equivalent_height(
+            frequency, PRESSURE_HPA, WATER_VAPOUR_DENSITY, TEMPERATURE_K
+        )
+        dry_height = float(np.asarray(heights[0]).ravel()[0])
+        wet_height = float(np.asarray(heights[1]).ravel()[0])
+        zenith = dry * dry_height + wet * wet_height
+
+        expected = {
+            "dry_db_per_km": dry,
+            "wet_db_per_km": wet,
+            "dry_height_km": dry_height,
+            "wet_height_km": wet_height,
+            "zenith_db": zenith,
+        }
+        for column, reference in expected.items():
+            actual = float(row[column])
+            absolute = abs(actual - reference)
+            relative = absolute / max(abs(reference), ABSOLUTE_TOLERANCE)
+            # The column passes on either measure, so the reported figure is
+            # the relative one and the acceptance test below also allows the
+            # absolute difference to be negligible.
+            if relative > worst[column]:
+                worst[column] = relative
+                worst_at[column] = frequency
+            if absolute > worst_absolute[column]:
+                worst_absolute[column] = absolute
+
+    print(f"frequencies compared: {len(rows)}")
+    for column in columns:
+        print(
+            f"{column:>15}: worst relative {worst[column]:.3e} "
+            f"absolute {worst_absolute[column]:.3e} at {worst_at[column]} GHz"
+        )
+    print(
+        f"{'tolerance':>15}:       {RELATIVE_TOLERANCE:.3e} relative or "
+        f"{ABSOLUTE_TOLERANCE:.3e} absolute"
+    )
+
+    passes = all(
+        worst[column] <= RELATIVE_TOLERANCE
+        or worst_absolute[column] <= ABSOLUTE_TOLERANCE
+        for column in columns
+    )
+    if passes:
+        print("result:          pass")
+        return 0
+    print("result:          FAIL")
+    return 1
+
+
+if __name__ == "__main__":
+    raise SystemExit(main(sys.argv))

validation/python/atmosphere/check_rain.py

@@ -0,0 +1,93 @@
+"""Cross-check the ApsisRS rain coefficients against ITU-R P.838-3.
+
+Reads the values written by ``cargo run -p apsis-rs --features rf --example
+dump_rain`` and recomputes each one with the ``itur`` package, which
+implements the recommendation directly.
+
+The script needs no network access.
+
+Usage:
+
+    pixi run python validation/python/atmosphere/check_rain.py [path]
+
+Exit status is zero when every value agrees within tolerance and one
+otherwise.
+"""
+
+from __future__ import annotations
+
+import csv
+import sys
+from pathlib import Path
+
+from itur.models.itu838 import _ITU838_3_
+
+# The implementation reproduces the reference to about twelve significant
+# digits. The Rust output is written with twelve decimal places, which limits
+# the comparison for the smallest coefficients, so the tolerance is set to
+# one part in ten to the eighth.
+RELATIVE_TOLERANCE = 1e-8
+
+DEFAULT_INPUT = Path("validation/data/expected/rust_rain.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("frequency_ghz")]
+
+
+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 > " + 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 = _ITU838_3_()
+    worst_k = 0.0
+    worst_alpha = 0.0
+    worst_at = None
+
+    for row in rows:
+        frequency = float(row["frequency_ghz"])
+        elevation = float(row["elevation_deg"])
+        polarisation = float(row["polarisation_deg"])
+        k, alpha = model.rain_specific_attenuation_coefficients(
+            frequency, elevation, polarisation
+        )
+        k = float(k)
+        alpha = float(alpha)
+
+        k_error = abs(float(row["k"]) - k) / k
+        alpha_error = abs(float(row["alpha"]) - alpha) / abs(alpha)
+        if k_error > worst_k or alpha_error > worst_alpha:
+            worst_at = (frequency, elevation, polarisation)
+            worst_k = max(worst_k, k_error)
+            worst_alpha = max(worst_alpha, alpha_error)
+
+    print(f"combinations compared: {len(rows)}")
+    print(f"worst relative k:      {worst_k:.3e}")
+    print(f"worst relative alpha:  {worst_alpha:.3e}")
+    print(f"worst combination:     {worst_at}")
+    print(f"tolerance:             {RELATIVE_TOLERANCE:.3e} relative")
+
+    if worst_k <= RELATIVE_TOLERANCE and worst_alpha <= RELATIVE_TOLERANCE:
+        print("result:                pass")
+        return 0
+    print("result:                FAIL")
+    return 1
+
+
+if __name__ == "__main__":
+    raise SystemExit(main(sys.argv))