rf: add the exact P.676 slant path with the P.835 atmosphere

crates/apsis-rs/Cargo.toml

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

crates/apsis-rs/examples/dump_gas_exact.rs

@@ -0,0 +1,40 @@
+//! Writes the exact slant path 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_exact \
+//!     > validation/data/expected/rust_gas_exact.csv
+//! ```
+
+use apsis_atmosphere::{Atmosphere, slant_path_attenuation_exact};
+
+/// The scenarios sampled, as frequency in gigahertz and elevation in
+/// degrees.
+const SCENARIOS: &[(f64, f64)] = &[
+    (12.0, 90.0),
+    (12.0, 30.0),
+    (12.0, 10.0),
+    (12.0, 5.0),
+    (12.0, 1.0),
+    (20.0, 30.0),
+    (20.0, 5.0),
+    (30.0, 10.0),
+    (60.0, 30.0),
+    (1.0, 30.0),
+    (100.0, 5.0),
+];
+
+fn main() {
+    println!("frequency_ghz,elevation_deg,attenuation_db");
+    let atmosphere = Atmosphere::REFERENCE;
+    for &(frequency, elevation_degrees) in SCENARIOS {
+        let Some(attenuation) =
+            slant_path_attenuation_exact(frequency, elevation_degrees.to_radians(), atmosphere)
+        else {
+            continue;
+        };
+        println!("{frequency},{elevation_degrees},{attenuation:.12}");
+    }
+}

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

@@ -168,3 +168,135 @@ pub fn slant_path_attenuation(
     }
     Some(zenith_attenuation(frequency_ghz, atmosphere)? / elevation.sin())
 }
+
+/// The number of layer boundaries used by the exact slant path integration.
+///
+/// Recommendation P.676 defines the boundaries by
+/// `delta_i = 1e-4 exp(i / 100)` for `i` from zero to 921, which gives 922
+/// boundaries. The layers between them are thin near the surface and grow
+/// with height, so the integration resolves the dense lower atmosphere where
+/// the absorption is concentrated.
+const EXACT_BOUNDARY_COUNT: usize = 922;
+
+/// The Earth radius used by the ray trace, in kilometers.
+///
+/// The recommendation uses the mean radius, not the geopotential radius of
+/// the standard atmosphere, because this value enters the geometry of the
+/// ray rather than the definition of the height.
+const RAY_TRACE_EARTH_RADIUS_KM: f64 = 6_371.0;
+
+/// Returns the attenuation along a slant path by the exact method of
+/// Recommendation P.676.
+///
+/// The exact method divides the atmosphere into a stack of layers whose
+/// thickness grows exponentially with height, evaluates the specific
+/// attenuation in each, and sums the contribution of each layer weighted by
+/// the ray path length through it. Because the ray is refracted, the path
+/// length is found by ray tracing with Snell's law applied to the ratio of
+/// the refractive indices of adjacent layers, which the approximate method
+/// ignores.
+///
+/// This is the form to use at low elevation, where the plane parallel
+/// approximation that [`slant_path_attenuation`] uses is inaccurate. It is
+/// far more expensive, so the approximate form remains available for the
+/// high elevation case where it is adequate.
+///
+/// The frequency is in gigahertz, the elevation in radians, and the surface
+/// atmosphere supplies the pressure, the water vapour density, and the
+/// temperature at ground level. The surface water vapour density scales the
+/// exponential vapour profile, so the atmosphere defines the whole stack.
+///
+/// Returns `None` when the frequency, elevation, or atmosphere is not
+/// usable.
+///
+/// Source: ITU-R Recommendation P.676-12, equations 13 through 19, with the
+/// layered atmosphere of ITU-R Recommendation P.835-6.
+pub fn slant_path_attenuation_exact(
+    frequency_ghz: f64,
+    elevation: f64,
+    atmosphere: Atmosphere,
+) -> Option<f64> {
+    if !frequency_ghz.is_finite() || frequency_ghz <= 0.0 || !atmosphere.is_valid() {
+        return None;
+    }
+    if !elevation.is_finite() || elevation <= 0.0 || elevation > core::f64::consts::FRAC_PI_2 {
+        return None;
+    }
+
+    // Build the stack of layer boundaries. Each entry holds the geometric
+    // height, the thickness, the profile values, and the refractive index.
+    let thickness_step = (1.0_f64 / 100.0).exp();
+    let mut heights = [0.0_f64; EXACT_BOUNDARY_COUNT];
+    let mut thicknesses = [0.0_f64; EXACT_BOUNDARY_COUNT];
+    let mut temperatures = [0.0_f64; EXACT_BOUNDARY_COUNT];
+    let mut pressures = [0.0_f64; EXACT_BOUNDARY_COUNT];
+    let mut densities = [0.0_f64; EXACT_BOUNDARY_COUNT];
+    let mut indices = [1.0_f64; EXACT_BOUNDARY_COUNT];
+
+    for index in 0..EXACT_BOUNDARY_COUNT {
+        let exponent = index as f64 / 100.0;
+        thicknesses[index] = 1.0e-4 * exponent.exp();
+        // Equation 15: the geometric height of the boundary.
+        heights[index] = 1.0e-4 * (exponent.exp() - 1.0) / (thickness_step - 1.0);
+
+        let temperature = crate::p835::standard_temperature(heights[index]);
+        let pressure = crate::p835::standard_pressure(heights[index]);
+        let density = crate::p835::standard_water_vapour_density(
+            heights[index],
+            atmosphere.water_vapour_density,
+        );
+        temperatures[index] = temperature;
+        pressures[index] = pressure;
+        densities[index] = density;
+
+        // The refractive index takes the total pressure, as the
+        // recommendation specifies, not the dry air pressure.
+        let vapour_pressure = density * temperature / 216.7;
+        indices[index] =
+            crate::p835::radio_refractive_index(pressure, vapour_pressure, temperature);
+    }
+
+    // The ray starts below the local horizon by the complement of the
+    // elevation, and is bent upward as it climbs.
+    let mut angle = core::f64::consts::FRAC_PI_2 - elevation;
+    let mut total = 0.0;
+
+    for index in 0..EXACT_BOUNDARY_COUNT {
+        let layer_atmosphere =
+            Atmosphere::new(pressures[index], densities[index], temperatures[index]);
+        let specific = crate::gaseous::specific_attenuation(frequency_ghz, layer_atmosphere)?;
+
+        // Equation 17: the path length through the layer.
+        let radius = RAY_TRACE_EARTH_RADIUS_KM + heights[index];
+        let delta = thicknesses[index];
+        let length = -radius * angle.cos()
+            + 0.5
+                * (4.0 * radius * radius * angle.cos().powi(2)
+                    + 8.0 * radius * delta
+                    + 4.0 * delta * delta)
+                    .sqrt();
+        total += length * specific;
+
+        // Equation 18a: the angle at which the ray leaves the layer.
+        let cosine_argument = ((-length * length - 2.0 * radius * delta - delta * delta)
+            / (2.0 * length * radius + 2.0 * length * delta))
+            .clamp(-1.0, 1.0);
+        let leave_angle = core::f64::consts::PI - cosine_argument.acos();
+
+        // Equation 19a: Snell's law at the boundary, using the ratio of the
+        // refractive index here to the one above. The topmost boundary has
+        // no layer above it, so its ratio is one.
+        let next = (index + 1).min(EXACT_BOUNDARY_COUNT - 1);
+        let ratio = if next == index {
+            1.0
+        } else {
+            indices[index] / indices[next]
+        };
+        angle = (leave_angle.sin() * ratio).clamp(-1.0, 1.0).asin();
+        if !angle.is_finite() {
+            break;
+        }
+    }
+
+    Some(total)
+}

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

@@ -50,6 +50,7 @@
 mod gaseous;
 mod gaseous_path;
 mod geometry;
+mod p835;
 mod p839;
 mod p839_table;
 mod rain;
@@ -59,9 +60,14 @@ pub use gaseous::{
     specific_attenuation_water_vapour,
 };
 pub use gaseous_path::{
-    EquivalentHeights, equivalent_heights, slant_path_attenuation, zenith_attenuation,
+    EquivalentHeights, equivalent_heights, slant_path_attenuation, slant_path_attenuation_exact,
+    zenith_attenuation,
 };
 pub use geometry::{air_mass, elevation_from_zenith, zenith_to_elevation};
+pub use p835::{
+    geopotential_height, radio_refractive_index, standard_pressure, standard_temperature,
+    standard_water_vapour_density,
+};
 pub use p839::{RAIN_HEIGHT_OFFSET_KM, grid_step_degrees, rain_height, zero_isotherm_height};
 pub use rain::{
     DEFAULT_RAIN_HEIGHT_KM, RainAttenuation, polarisation_tilt, rain_attenuation,

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

@@ -0,0 +1,148 @@
+//! The reference standard atmosphere of ITU-R Recommendation P.835.
+//!
+//! The exact slant path method of Recommendation P.676 integrates the
+//! specific attenuation through a stack of layers, so it needs the
+//! temperature, pressure, and water vapour density at each layer height.
+//! This module supplies those from the mean annual global reference
+//! atmosphere of Recommendation P.835.
+//!
+//! Height is measured as geopotential height up to the base of the
+//! thermosphere. The geopotential height is the geometric height adjusted
+//! for the variation of gravity with altitude, and it is the independent
+//! variable the temperature and pressure profiles are expressed in.
+//!
+//! Source: ITU-R Recommendation P.835-6, *Reference Standard Atmospheres*,
+//! section 1, the mean annual global reference atmosphere.
+
+/// The radius of the Earth used in the geopotential conversion, in
+/// kilometers.
+///
+/// The recommendation fixes this value for the standard atmosphere, which is
+/// slightly different from the polar radius used elsewhere. It is the value
+/// that makes the geopotential height consistent with the tabulated
+/// profiles.
+const GEOPOTENTIAL_EARTH_RADIUS_KM: f64 = 6_356.766;
+
+/// Returns the geopotential height for a geometric height, in kilometers.
+///
+/// The relation is
+///
+/// ```text
+/// h_p = R H / (R + H)
+/// ```
+///
+/// with `R` the fixed Earth radius of the standard atmosphere and `H` the
+/// geometric height. The geopotential height is always the smaller of the
+/// two.
+pub fn geopotential_height(geometric_height_km: f64) -> f64 {
+    GEOPOTENTIAL_EARTH_RADIUS_KM * geometric_height_km
+        / (GEOPOTENTIAL_EARTH_RADIUS_KM + geometric_height_km)
+}
+
+/// Returns the temperature of the reference atmosphere, in kelvin.
+///
+/// The profile is piecewise linear in geopotential height up to 84.852 km
+/// and is then defined by two further segments that extend to 100 km. Above
+/// 100 km the recommendation does not define a global profile and this
+/// function returns the value at 100 km.
+///
+/// The geometric height is in kilometers.
+pub fn standard_temperature(geometric_height_km: f64) -> f64 {
+    let height = geopotential_height(geometric_height_km);
+    if height <= 11.0 {
+        288.15 - 6.5 * height
+    } else if height <= 20.0 {
+        216.65
+    } else if height <= 32.0 {
+        216.65 + (height - 20.0)
+    } else if height <= 47.0 {
+        228.65 + 2.8 * (height - 32.0)
+    } else if height <= 51.0 {
+        270.65
+    } else if height <= 71.0 {
+        270.65 - 2.8 * (height - 51.0)
+    } else if height <= 84.852 {
+        214.65 - 2.0 * (height - 71.0)
+    } else if (86.0..=91.0).contains(&geometric_height_km) {
+        186.867_3
+    } else if geometric_height_km > 91.0 && geometric_height_km <= 100.0 {
+        let ratio = (geometric_height_km - 91.0) / 19.942_9;
+        263.190_5 - 76.323_2 * (1.0 - ratio * ratio).max(0.0).sqrt()
+    } else {
+        195.081_34
+    }
+}
+
+/// Returns the pressure of the reference atmosphere, in hectopascals.
+///
+/// The profile is piecewise, continuous at the segment boundaries, and
+/// follows from the hydrostatic relation applied to each temperature
+/// segment.
+///
+/// The geometric height is in kilometers.
+pub fn standard_pressure(geometric_height_km: f64) -> f64 {
+    let height = geopotential_height(geometric_height_km);
+    if height <= 11.0 {
+        1_013.25 * (288.15 / (288.15 - 6.5 * height)).powf(-34.163_2 / 6.5)
+    } else if height <= 20.0 {
+        226.322_6 * (-34.163_2 * (height - 11.0) / 216.65).exp()
+    } else if height <= 32.0 {
+        54.749_80 * (216.65 / (216.65 + (height - 20.0))).powf(34.163_2)
+    } else if height <= 47.0 {
+        8.680_422 * (228.65 / (228.65 + 2.8 * (height - 32.0))).powf(34.163_2 / 2.8)
+    } else if height <= 51.0 {
+        1.109_106 * (-34.163_2 * (height - 47.0) / 270.65).exp()
+    } else if height <= 71.0 {
+        0.669_416_7 * (270.65 / (270.65 - 2.8 * (height - 51.0))).powf(-34.163_2 / 2.8)
+    } else if height <= 84.852 {
+        0.039_566_49 * (214.65 / (214.65 - 2.0 * (height - 71.0))).powf(-34.163_2 / 2.0)
+    } else if (86.0..=100.0).contains(&geometric_height_km) {
+        let h = geometric_height_km;
+        (95.571_899 - 4.011_801 * h + 6.424_731e-2 * h * h - 4.789_660e-4 * h * h * h
+            + 1.340_543e-6 * h * h * h * h)
+            .exp()
+    } else {
+        1.0e-62
+    }
+}
+
+/// Returns the water vapour density of the reference atmosphere, in grams
+/// per cubic meter.
+///
+/// The profile is exponential with a scale height of two kilometers, scaled
+/// so that the density at the surface is the value given.
+///
+/// The geometric height is in kilometers.
+pub fn standard_water_vapour_density(geometric_height_km: f64, surface_density: f64) -> f64 {
+    surface_density * (-geometric_height_km / 2.0).exp()
+}
+
+/// Returns the radio refractive index of moist air, dimensionless.
+///
+/// The relation is
+///
+/// ```text
+/// N = 77.6 P_d / T + 72 e / T + 3.75e5 e / T^2
+/// n = 1 + N 1e-6
+/// ```
+///
+/// with `P_d` the dry air pressure and `e` the water vapour partial
+/// pressure, both in hectopascals, and `T` the temperature in kelvin. The
+/// first two terms are the non-dispersive refractivity and the third is the
+/// dipole term that gives water vapour its strong effect.
+///
+/// Source: ITU-R Recommendation P.453, which gives this expression for the
+/// refractive index of moist air.
+pub fn radio_refractive_index(
+    dry_pressure_hpa: f64,
+    vapour_pressure_hpa: f64,
+    temperature_k: f64,
+) -> f64 {
+    if temperature_k == 0.0 {
+        return 1.0;
+    }
+    let refractivity = 77.6 * dry_pressure_hpa / temperature_k
+        + 72.0 * vapour_pressure_hpa / temperature_k
+        + 3.75e5 * vapour_pressure_hpa / (temperature_k * temperature_k);
+    1.0 + refractivity * 1.0e-6
+}

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

@@ -255,3 +255,79 @@ fn slant_path_attenuation_rejects_an_elevation_at_the_horizon() {
     assert_eq!(slant_path_attenuation(12.0, -0.1, reference()), None);
     assert_eq!(slant_path_attenuation(12.0, f64::NAN, reference()), None);
 }
+
+#[test]
+fn the_exact_slant_path_matches_the_recommendation() {
+    use apsis_atmosphere::slant_path_attenuation_exact;
+    // These values are from a reference implementation of the exact method
+    // of ITU-R Recommendation P.676-12.
+    let cases: [(f64, f64, f64); 10] = [
+        (12.0, 90.0, 0.059926619268),
+        (12.0, 30.0, 0.119683280157),
+        (12.0, 10.0, 0.340072667215),
+        (12.0, 5.0, 0.651451990239),
+        (12.0, 1.0, 2.054182412438),
+        (20.0, 30.0, 0.504703085886),
+        (20.0, 5.0, 2.793084848959),
+        (30.0, 10.0, 1.319753394647),
+        (60.0, 30.0, 309.245091571274),
+        (100.0, 5.0, 10.109978745794),
+    ];
+    for (frequency, elevation_degrees, expected) in cases {
+        let actual =
+            slant_path_attenuation_exact(frequency, elevation_degrees.to_radians(), reference())
+                .expect("valid arguments");
+        let relative = (actual - expected).abs() / expected;
+        assert!(
+            relative < 1e-6,
+            "{frequency} GHz at {elevation_degrees} degrees: got {actual}, expected {expected}"
+        );
+    }
+}
+
+#[test]
+fn the_exact_and_approximate_methods_agree_at_high_elevation() {
+    use apsis_atmosphere::{slant_path_attenuation, slant_path_attenuation_exact};
+    // The two methods use different formulations, so they agree to a few
+    // per cent rather than exactly even at the zenith, where the plane
+    // parallel geometry is exact. The recommendation notes that the
+    // approximate method is intended only for the elevation range above
+    // about five degrees and for rough work.
+    let exact = slant_path_attenuation_exact(12.0, core::f64::consts::FRAC_PI_2, reference())
+        .expect("valid");
+    let approximate =
+        slant_path_attenuation(12.0, core::f64::consts::FRAC_PI_2, reference()).expect("valid");
+    let relative = (exact - approximate).abs() / approximate;
+    assert!(relative < 0.05, "exact {exact}, approximate {approximate}");
+}
+
+#[test]
+fn the_exact_method_exceeds_the_approximate_one_at_low_elevation() {
+    use apsis_atmosphere::{slant_path_attenuation, slant_path_attenuation_exact};
+    // At a low elevation the curved path through the thinning atmosphere is
+    // shorter than the plane parallel estimate, and the refraction bends
+    // the ray downward, so the two methods diverge.
+    let elevation = 2.0_f64.to_radians();
+    let exact = slant_path_attenuation_exact(12.0, elevation, reference()).expect("valid");
+    let approximate = slant_path_attenuation(12.0, elevation, reference()).expect("valid");
+    assert!(
+        (exact - approximate).abs() / approximate > 0.001,
+        "exact {exact}, approximate {approximate}"
+    );
+}
+
+#[test]
+fn the_exact_method_rejects_unusable_arguments() {
+    use apsis_atmosphere::slant_path_attenuation_exact;
+    assert_eq!(slant_path_attenuation_exact(0.0, 1.0, reference()), None);
+    assert_eq!(slant_path_attenuation_exact(12.0, 0.0, reference()), None);
+    assert_eq!(slant_path_attenuation_exact(12.0, -0.1, reference()), None);
+    assert_eq!(
+        slant_path_attenuation_exact(12.0, f64::NAN, reference()),
+        None
+    );
+    assert_eq!(
+        slant_path_attenuation_exact(12.0, 1.0, Atmosphere::new(0.0, 7.5, 288.15)),
+        None
+    );
+}

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

@@ -0,0 +1,153 @@
+//! Tests for the ITU-R P.835 reference standard atmosphere.
+
+use apsis_atmosphere::{
+    geopotential_height, radio_refractive_index, standard_pressure, standard_temperature,
+    standard_water_vapour_density,
+};
+
+const TOLERANCE: f64 = 1e-9;
+
+#[test]
+fn the_geopotential_height_is_below_the_geometric_height() {
+    // The conversion compresses the height because gravity weakens with
+    // altitude, so the geopotential height is always the smaller value.
+    for geometric in [0.0, 1.0, 11.0, 50.0, 100.0] {
+        let geopotential = geopotential_height(geometric);
+        assert!(geopotential <= geometric, "height {geometric}");
+    }
+    // At the surface the two coincide.
+    assert_eq!(geopotential_height(0.0), 0.0);
+}
+
+#[test]
+fn the_surface_conditions_match_the_recommendation() {
+    // The mean annual global reference atmosphere starts at 288.15 K and
+    // 1013.25 hPa.
+    assert!((standard_temperature(0.0) - 288.15).abs() < TOLERANCE);
+    assert!((standard_pressure(0.0) - 1_013.25).abs() < TOLERANCE);
+    assert!((standard_water_vapour_density(0.0, 7.5) - 7.5).abs() < TOLERANCE);
+}
+
+#[test]
+fn the_temperature_profile_has_the_tabulated_segments() {
+    // The troposphere cools at 6.5 K per kilometer of geopotential height,
+    // so the rate is expressed against that height, not the geometric one.
+    let expected = 288.15 - 6.5 * geopotential_height(1.0);
+    assert!((standard_temperature(1.0) - expected).abs() < TOLERANCE);
+
+    // The stratosphere is isothermal at 216.65 K between 11 and 20 km of
+    // geopotential height. A geometric height of 15 km has a geopotential
+    // height of about 14.97 km, which is inside that band.
+    let geopotential = geopotential_height(15.0);
+    assert!(
+        (11.0..=20.0).contains(&geopotential),
+        "height {geopotential}"
+    );
+    assert!((standard_temperature(15.0) - 216.65).abs() < TOLERANCE);
+
+    // Above 20 km geopotential the temperature rises at 1 K per kilometer.
+    // A geometric height of 30 km reaches the middle of that segment.
+    let geopotential = geopotential_height(30.0);
+    assert!(
+        (20.0..32.0).contains(&geopotential),
+        "height {geopotential}"
+    );
+    let expected = 216.65 + (geopotential - 20.0);
+    assert!((standard_temperature(30.0) - expected).abs() < TOLERANCE);
+
+    // The temperature is continuous across each segment boundary. The
+    // boundaries are in geopotential height, so they are converted to the
+    // geometric heights at which they occur before sampling either side.
+    for geopotential in [11.0_f64, 20.0, 32.0, 47.0, 51.0, 71.0] {
+        let radius = 6_356.766;
+        let geometric = radius * geopotential / (radius - geopotential);
+        let below = standard_temperature(geometric - 0.001);
+        let above = standard_temperature(geometric + 0.001);
+        assert!(
+            (below - above).abs() < 0.01,
+            "boundary {geopotential}: {below} against {above}"
+        );
+    }
+}
+
+#[test]
+fn the_temperature_is_continuous_across_every_segment() {
+    for boundary in [11.0_f64, 20.0, 32.0, 47.0, 51.0, 71.0, 84.852] {
+        // The boundaries are in geopotential height, so a geometric height
+        // just either side is used.
+        let below = standard_temperature(boundary - 0.01);
+        let above = standard_temperature(boundary + 0.01);
+        assert!(
+            (below - above).abs() < 0.2,
+            "boundary {boundary}: {below} against {above}"
+        );
+    }
+}
+
+#[test]
+fn the_pressure_falls_with_height() {
+    let mut previous = f64::INFINITY;
+    for height in [0.0, 1.0, 5.0, 11.0, 20.0, 32.0, 47.0, 51.0, 71.0, 84.852] {
+        let pressure = standard_pressure(height);
+        assert!(pressure < previous, "height {height} gave {pressure}");
+        previous = pressure;
+    }
+}
+
+#[test]
+fn the_pressure_is_continuous_across_the_segments() {
+    for boundary in [11.0_f64, 20.0, 32.0, 47.0, 51.0, 71.0, 84.852] {
+        let below = standard_pressure(boundary - 0.01);
+        let above = standard_pressure(boundary + 0.01);
+        let relative = (below - above).abs() / below;
+        assert!(
+            relative < 0.01,
+            "boundary {boundary}: {below} against {above}, relative {relative}"
+        );
+    }
+}
+
+#[test]
+fn the_pressure_at_the_tropopause_matches_the_reference() {
+    // The tropopause is at 11 km of geopotential height, which is a
+    // geometric height of about 11.02 km. The pressure there is 226.32 hPa.
+    let geometric = 11.019_068;
+    assert!((geopotential_height(geometric) - 11.0).abs() < 1e-5);
+    let pressure = standard_pressure(geometric);
+    assert!((pressure - 226.32).abs() < 0.01, "pressure {pressure}");
+}
+
+#[test]
+fn the_water_vapour_profile_is_exponential() {
+    // The scale height is two kilometers, so the density halves about every
+    // 1.386 km.
+    let surface = 7.5;
+    let at_two = standard_water_vapour_density(2.0, surface);
+    assert!((at_two - surface / core::f64::consts::E).abs() < 1e-9);
+    // The profile scales with the surface value.
+    let doubled = standard_water_vapour_density(2.0, 15.0);
+    assert!((doubled / at_two - 2.0).abs() < 1e-12);
+}
+
+#[test]
+fn the_radio_refractive_index_at_the_surface_is_about_one() {
+    // A refractivity of about three hundred N units gives an index about
+    // 1.0003.
+    let index = radio_refractive_index(1_013.25, 10.0, 288.15);
+    assert!(index > 1.0 && index < 1.001, "index {index}");
+    assert!((index - 1.000_315).abs() < 1e-5, "index {index}");
+}
+
+#[test]
+fn the_radio_refractive_index_grows_with_humidity() {
+    let dry = radio_refractive_index(1_013.25, 0.0, 288.15);
+    let moist = radio_refractive_index(1_013.25, 20.0, 288.15);
+    assert!(moist > dry);
+}
+
+#[test]
+fn a_zero_temperature_is_handled_without_dividing_by_zero() {
+    // The expression divides by the temperature, so a zero would give an
+    // infinite value. The function returns the index of vacuum instead.
+    assert_eq!(radio_refractive_index(1_013.25, 10.0, 0.0), 1.0);
+}

validation/python/atmosphere/check_gas_exact.py

@@ -0,0 +1,93 @@
+"""Cross-check the ApsisRS exact slant path against ITU-R P.676-12.
+
+Reads the values written by ``cargo run -p apsis-rs --features rf --example
+dump_gas_exact`` and recomputes each one with the ``itur`` package, which
+implements the exact method of the recommendation directly. The exact method
+integrates the specific attenuation through a layered standard atmosphere and
+traces the refracted ray, so it is the form to use at low elevation.
+
+The script needs no network access.
+
+Usage:
+
+    pixi run python validation/python/atmosphere/check_gas_exact.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.itu676 import _ITU676_12_
+
+# The reference atmosphere, matching the Rust side.
+PRESSURE_HPA = 1013.25
+WATER_VAPOUR_DENSITY = 7.5
+TEMPERATURE_K = 288.15
+
+# The implementation reproduces the reference to about nine significant
+# digits, so this allows for the rounding of the recorded values.
+RELATIVE_TOLERANCE = 1e-7
+
+DEFAULT_INPUT = Path("validation/data/expected/rust_gas_exact.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_exact > " + 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_()
+    worst = 0.0
+    worst_at = None
+
+    for row in rows:
+        frequency = float(row["frequency_ghz"])
+        elevation = float(row["elevation_deg"])
+        reference = model.gaseous_attenuation_slant_path(
+            frequency, elevation, WATER_VAPOUR_DENSITY, PRESSURE_HPA, TEMPERATURE_K,
+            mode="exact",
+        )
+        reference = float(reference.value if hasattr(reference, "value") else reference)
+
+        actual = float(row["attenuation_db"])
+        relative = abs(actual - reference) / max(abs(reference), 1e-12)
+        if relative > worst:
+            worst = relative
+            worst_at = (frequency, elevation)
+
+    print(f"scenarios compared: {len(rows)}")
+    print(f"worst relative:     {worst:.3e} at {worst_at}")
+    print(f"tolerance:          {RELATIVE_TOLERANCE:.3e} relative")
+
+    if worst <= RELATIVE_TOLERANCE:
+        print("result:             pass")
+        return 0
+    print("result:             FAIL")
+    return 1
+
+
+if __name__ == "__main__":
+    raise SystemExit(main(sys.argv))