From be00a4f08893b9c6eae3ef55cc012b90bbaa792f Mon Sep 17 00:00:00 2001 From: ApsisRS contributors Date: Sat, 26 Sep 2026 15:43:12 +0800 Subject: [PATCH] rf: add the exact P.676 slant path with the P.835 atmosphere --- crates/apsis-rs/Cargo.toml | 4 + crates/apsis-rs/examples/dump_gas_exact.rs | 40 +++++ .../rf/apsis-atmosphere/src/gaseous_path.rs | 132 +++++++++++++++ crates/rf/apsis-atmosphere/src/lib.rs | 8 +- crates/rf/apsis-atmosphere/src/p835.rs | 148 +++++++++++++++++ crates/rf/apsis-atmosphere/tests/gaseous.rs | 76 +++++++++ crates/rf/apsis-atmosphere/tests/p835.rs | 153 ++++++++++++++++++ .../python/atmosphere/check_gas_exact.py | 93 +++++++++++ 8 files changed, 653 insertions(+), 1 deletion(-) create mode 100644 crates/apsis-rs/examples/dump_gas_exact.rs create mode 100644 crates/rf/apsis-atmosphere/src/p835.rs create mode 100644 crates/rf/apsis-atmosphere/tests/p835.rs create mode 100644 validation/python/atmosphere/check_gas_exact.py diff --git a/crates/apsis-rs/Cargo.toml b/crates/apsis-rs/Cargo.toml index 5fd5f43..7d86e17 100644 --- a/crates/apsis-rs/Cargo.toml +++ b/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"] diff --git a/crates/apsis-rs/examples/dump_gas_exact.rs b/crates/apsis-rs/examples/dump_gas_exact.rs new file mode 100644 index 0000000..eae3a5d --- /dev/null +++ b/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}"); + } +} diff --git a/crates/rf/apsis-atmosphere/src/gaseous_path.rs b/crates/rf/apsis-atmosphere/src/gaseous_path.rs index ce9a978..f2705de 100644 --- a/crates/rf/apsis-atmosphere/src/gaseous_path.rs +++ b/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 { + 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) +} diff --git a/crates/rf/apsis-atmosphere/src/lib.rs b/crates/rf/apsis-atmosphere/src/lib.rs index 03e5f10..3dd59c0 100644 --- a/crates/rf/apsis-atmosphere/src/lib.rs +++ b/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, diff --git a/crates/rf/apsis-atmosphere/src/p835.rs b/crates/rf/apsis-atmosphere/src/p835.rs new file mode 100644 index 0000000..2e7f173 --- /dev/null +++ b/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 +} diff --git a/crates/rf/apsis-atmosphere/tests/gaseous.rs b/crates/rf/apsis-atmosphere/tests/gaseous.rs index feee915..216335d 100644 --- a/crates/rf/apsis-atmosphere/tests/gaseous.rs +++ b/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 + ); +} diff --git a/crates/rf/apsis-atmosphere/tests/p835.rs b/crates/rf/apsis-atmosphere/tests/p835.rs new file mode 100644 index 0000000..a8c03e8 --- /dev/null +++ b/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); +} diff --git a/validation/python/atmosphere/check_gas_exact.py b/validation/python/atmosphere/check_gas_exact.py new file mode 100644 index 0000000..c7e0440 --- /dev/null +++ b/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))