From 2f342db8770e399c39c5518e82c20aa0d70465e4 Mon Sep 17 00:00:00 2001 From: ApsisRS contributors Date: Sat, 26 Sep 2026 15:13:48 +0800 Subject: [PATCH] rf: replace gaseous and rain models with ITU-R P.676-12 and P.838-3 --- crates/apsis-rs/Cargo.toml | 8 + crates/apsis-rs/examples/dump_gas.rs | 42 + crates/apsis-rs/examples/dump_rain.rs | 36 + crates/rf/apsis-atmosphere/src/gaseous.rs | 1004 +++++++++++++++-- .../rf/apsis-atmosphere/src/gaseous_path.rs | 170 +++ crates/rf/apsis-atmosphere/src/lib.rs | 62 +- crates/rf/apsis-atmosphere/src/rain.rs | 483 +++++--- crates/rf/apsis-atmosphere/tests/gaseous.rs | 288 +++-- crates/rf/apsis-atmosphere/tests/rain.rs | 390 +++++-- validation/pixi.lock | 177 +++ validation/pixi.toml | 1 + validation/python/atmosphere/check_gas.py | 136 +++ validation/python/atmosphere/check_rain.py | 93 ++ 13 files changed, 2466 insertions(+), 424 deletions(-) create mode 100644 crates/apsis-rs/examples/dump_gas.rs create mode 100644 crates/apsis-rs/examples/dump_rain.rs create mode 100644 crates/rf/apsis-atmosphere/src/gaseous_path.rs create mode 100644 validation/python/atmosphere/check_gas.py create mode 100644 validation/python/atmosphere/check_rain.py diff --git a/crates/apsis-rs/Cargo.toml b/crates/apsis-rs/Cargo.toml index 80512e6..7b26d0c 100644 --- a/crates/apsis-rs/Cargo.toml +++ b/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"] diff --git a/crates/apsis-rs/examples/dump_gas.rs b/crates/apsis-rs/examples/dump_gas.rs new file mode 100644 index 0000000..45be3d1 --- /dev/null +++ b/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 + ); + } +} diff --git a/crates/apsis-rs/examples/dump_rain.rs b/crates/apsis-rs/examples/dump_rain.rs new file mode 100644 index 0000000..61ecdb1 --- /dev/null +++ b/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}"); + } + } +} diff --git a/crates/rf/apsis-atmosphere/src/gaseous.rs b/crates/rf/apsis-atmosphere/src/gaseous.rs index 8ac4aeb..ccb22d6 100644 --- a/crates/rf/apsis-atmosphere/src/gaseous.rs +++ b/crates/rf/apsis-atmosphere/src/gaseous.rs @@ -1,123 +1,947 @@ -//! Gaseous attenuation. +//! Specific attenuation of atmospheric gases. +//! +//! The attenuation caused by oxygen and water vapour is computed by the +//! line-by-line summation of ITU-R Recommendation P.676. Each absorption +//! line contributes a resonance whose strength depends on the temperature +//! and on the pressure through the line width and the pressure shift. +//! +//! The model is parameterised by four atmospheric quantities: +//! +//! - The dry air pressure, in hectopascals. +//! - The water vapour density, in grams per cubic meter. +//! - The temperature, in kelvin. +//! - The frequency, in gigahertz. +//! +//! The pressure is the dry air pressure, not the total barometric pressure. +//! The partial pressure of water vapour follows from the density and the +//! temperature, and the two are used separately in the line widths. +//! +//! Source: ITU-R Recommendation P.676-12, *Attenuation by atmospheric +//! gases*, sections 1 and 2. The line tables are the ones tabulated in +//! tables 1 and 2 of that recommendation, and the expressions for the line +//! widths and the strengths are equations 1 through 8. -use crate::geometry::air_mass; -use apsis_rf::Decibel; +use core::f64::consts::E; -/// The zenith attenuation of a gaseous absorber, split into its two parts. +/// One oxygen absorption line. /// -/// Oxygen and water vapour are carried separately because they respond very -/// differently to frequency and to the weather: the water vapour term -/// depends on the humidity and the oxygen term does not. -#[derive(Debug, Clone, Copy, PartialEq, Default)] -pub struct GaseousAttenuation { - /// The oxygen contribution at the zenith, in decibels. - pub oxygen_db: f64, - /// The water vapour contribution at the zenith, in decibels. - pub water_vapour_db: f64, +/// The coefficients are the ones tabulated in table 1 of ITU-R P.676-12. The +/// first field is the line frequency in gigahertz and the rest are the +/// dimensionless coefficients of the line. +#[derive(Debug, Clone, Copy)] +struct OxygenLine { + /// The line frequency, in gigahertz. + f0: f64, + /// Line strength coefficient. + a1: f64, + /// Temperature dependence of the line strength. + a2: f64, + /// Line width coefficient. + a3: f64, + /// Pressure dependence of the line width. + a4: f64, + /// Pressure shift coefficient. + a5: f64, + /// Mixing ratio dependence of the pressure shift. + a6: f64, } -impl GaseousAttenuation { - /// Creates a gaseous attenuation from its two zenith components. - pub const fn new(oxygen_db: f64, water_vapour_db: f64) -> Self { +/// One water vapour absorption line. +/// +/// The coefficients are the ones tabulated in table 2 of ITU-R P.676-12. +#[derive(Debug, Clone, Copy)] +struct WaterLine { + /// The line frequency, in gigahertz. + f0: f64, + /// Line strength coefficient. + b1: f64, + /// Temperature dependence of the line strength. + b2: f64, + /// Line width coefficient. + b3: f64, + /// Pressure dependence of the line width. + b4: f64, + /// Water vapour dependence of the line width. + b5: f64, + /// Temperature dependence of the water vapour term. + b6: f64, +} + +/// The oxygen absorption lines of table 1 of ITU-R P.676-12. +/// +/// The first thirty seven entries are the fine structure of the 60 GHz +/// complex, followed by the isolated line at 118.75 GHz and six lines in the +/// submillimetre region. +const OXYGEN_LINES: &[OxygenLine] = &[ + OxygenLine { + f0: 50.474_214, + a1: 0.975_000, + a2: 9.651_000, + a3: 6.690_000, + a4: 0.000_000, + a5: 2.566_000, + a6: 6.850_000, + }, + OxygenLine { + f0: 50.987_745, + a1: 2.529_000, + a2: 8.653_000, + a3: 7.170_000, + a4: 0.000_000, + a5: 2.246_000, + a6: 6.800_000, + }, + OxygenLine { + f0: 51.503_360, + a1: 6.193_000, + a2: 7.709_000, + a3: 7.640_000, + a4: 0.000_000, + a5: 1.947_000, + a6: 6.729_000, + }, + OxygenLine { + f0: 52.021_429, + a1: 14.320_000, + a2: 6.819_000, + a3: 8.110_000, + a4: 0.000_000, + a5: 1.667_000, + a6: 6.640_000, + }, + OxygenLine { + f0: 52.542_418, + a1: 31.240_000, + a2: 5.983_000, + a3: 8.580_000, + a4: 0.000_000, + a5: 1.388_000, + a6: 6.526_000, + }, + OxygenLine { + f0: 53.066_934, + a1: 64.290_000, + a2: 5.201_000, + a3: 9.060_000, + a4: 0.000_000, + a5: 1.349_000, + a6: 6.206_000, + }, + OxygenLine { + f0: 53.595_775, + a1: 124.600_000, + a2: 4.474_000, + a3: 9.550_000, + a4: 0.000_000, + a5: 2.227_000, + a6: 5.085_000, + }, + OxygenLine { + f0: 54.130_025, + a1: 227.300_000, + a2: 3.800_000, + a3: 9.960_000, + a4: 0.000_000, + a5: 3.170_000, + a6: 3.750_000, + }, + OxygenLine { + f0: 54.671_180, + a1: 389.700_000, + a2: 3.182_000, + a3: 10.370_000, + a4: 0.000_000, + a5: 3.558_000, + a6: 2.654_000, + }, + OxygenLine { + f0: 55.221_384, + a1: 627.100_000, + a2: 2.618_000, + a3: 10.890_000, + a4: 0.000_000, + a5: 2.560_000, + a6: 2.952_000, + }, + OxygenLine { + f0: 55.783_815, + a1: 945.300_000, + a2: 2.109_000, + a3: 11.340_000, + a4: 0.000_000, + a5: -1.172_000, + a6: 6.135_000, + }, + OxygenLine { + f0: 56.264_774, + a1: 543.400_000, + a2: 0.014_000, + a3: 17.030_000, + a4: 0.000_000, + a5: 3.525_000, + a6: -0.978_000, + }, + OxygenLine { + f0: 56.363_399, + a1: 1_331.800_000, + a2: 1.654_000, + a3: 11.890_000, + a4: 0.000_000, + a5: -2.378_000, + a6: 6.547_000, + }, + OxygenLine { + f0: 56.968_211, + a1: 1_746.600_000, + a2: 1.255_000, + a3: 12.230_000, + a4: 0.000_000, + a5: -3.545_000, + a6: 6.451_000, + }, + OxygenLine { + f0: 57.612_486, + a1: 2_120.100_000, + a2: 0.910_000, + a3: 12.620_000, + a4: 0.000_000, + a5: -5.416_000, + a6: 6.056_000, + }, + OxygenLine { + f0: 58.323_877, + a1: 2_363.700_000, + a2: 0.621_000, + a3: 12.950_000, + a4: 0.000_000, + a5: -1.932_000, + a6: 0.436_000, + }, + OxygenLine { + f0: 58.446_588, + a1: 1_442.100_000, + a2: 0.083_000, + a3: 14.910_000, + a4: 0.000_000, + a5: 6.768_000, + a6: -1.273_000, + }, + OxygenLine { + f0: 59.164_204, + a1: 2_379.900_000, + a2: 0.387_000, + a3: 13.530_000, + a4: 0.000_000, + a5: -6.561_000, + a6: 2.309_000, + }, + OxygenLine { + f0: 59.590_983, + a1: 2_090.700_000, + a2: 0.207_000, + a3: 14.080_000, + a4: 0.000_000, + a5: 6.957_000, + a6: -0.776_000, + }, + OxygenLine { + f0: 60.306_056, + a1: 2_103.400_000, + a2: 0.207_000, + a3: 14.150_000, + a4: 0.000_000, + a5: -6.395_000, + a6: 0.699_000, + }, + OxygenLine { + f0: 60.434_778, + a1: 2_438.000_000, + a2: 0.386_000, + a3: 13.390_000, + a4: 0.000_000, + a5: 6.342_000, + a6: -2.825_000, + }, + OxygenLine { + f0: 61.150_562, + a1: 2_479.500_000, + a2: 0.621_000, + a3: 12.920_000, + a4: 0.000_000, + a5: 1.014_000, + a6: -0.584_000, + }, + OxygenLine { + f0: 61.800_158, + a1: 2_275.900_000, + a2: 0.910_000, + a3: 12.630_000, + a4: 0.000_000, + a5: 5.014_000, + a6: -6.619_000, + }, + OxygenLine { + f0: 62.411_220, + a1: 1_915.400_000, + a2: 1.255_000, + a3: 12.170_000, + a4: 0.000_000, + a5: 3.029_000, + a6: -6.759_000, + }, + OxygenLine { + f0: 62.486_253, + a1: 1_503.000_000, + a2: 0.083_000, + a3: 15.130_000, + a4: 0.000_000, + a5: -4.499_000, + a6: 0.844_000, + }, + OxygenLine { + f0: 62.997_984, + a1: 1_490.200_000, + a2: 1.654_000, + a3: 11.740_000, + a4: 0.000_000, + a5: 1.856_000, + a6: -6.675_000, + }, + OxygenLine { + f0: 63.568_526, + a1: 1_078.000_000, + a2: 2.108_000, + a3: 11.340_000, + a4: 0.000_000, + a5: 0.658_000, + a6: -6.139_000, + }, + OxygenLine { + f0: 64.127_775, + a1: 728.700_000, + a2: 2.617_000, + a3: 10.880_000, + a4: 0.000_000, + a5: -3.036_000, + a6: -2.895_000, + }, + OxygenLine { + f0: 64.678_910, + a1: 461.300_000, + a2: 3.181_000, + a3: 10.380_000, + a4: 0.000_000, + a5: -3.968_000, + a6: -2.590_000, + }, + OxygenLine { + f0: 65.224_078, + a1: 274.000_000, + a2: 3.800_000, + a3: 9.960_000, + a4: 0.000_000, + a5: -3.528_000, + a6: -3.680_000, + }, + OxygenLine { + f0: 65.764_779, + a1: 153.000_000, + a2: 4.473_000, + a3: 9.550_000, + a4: 0.000_000, + a5: -2.548_000, + a6: -5.002_000, + }, + OxygenLine { + f0: 66.302_096, + a1: 80.400_000, + a2: 5.200_000, + a3: 9.060_000, + a4: 0.000_000, + a5: -1.660_000, + a6: -6.091_000, + }, + OxygenLine { + f0: 66.836_834, + a1: 39.800_000, + a2: 5.982_000, + a3: 8.580_000, + a4: 0.000_000, + a5: -1.680_000, + a6: -6.393_000, + }, + OxygenLine { + f0: 67.369_601, + a1: 18.560_000, + a2: 6.818_000, + a3: 8.110_000, + a4: 0.000_000, + a5: -1.956_000, + a6: -6.475_000, + }, + OxygenLine { + f0: 67.900_868, + a1: 8.172_000, + a2: 7.708_000, + a3: 7.640_000, + a4: 0.000_000, + a5: -2.216_000, + a6: -6.545_000, + }, + OxygenLine { + f0: 68.431_006, + a1: 3.397_000, + a2: 8.652_000, + a3: 7.170_000, + a4: 0.000_000, + a5: -2.492_000, + a6: -6.600_000, + }, + OxygenLine { + f0: 68.960_312, + a1: 1.334_000, + a2: 9.650_000, + a3: 6.690_000, + a4: 0.000_000, + a5: -2.773_000, + a6: -6.650_000, + }, + OxygenLine { + f0: 118.750_334, + a1: 940.300_000, + a2: 0.010_000, + a3: 16.640_000, + a4: 0.000_000, + a5: -0.439_000, + a6: 0.079_000, + }, + OxygenLine { + f0: 368.498_246, + a1: 67.400_000, + a2: 0.048_000, + a3: 16.400_000, + a4: 0.000_000, + a5: 0.000_000, + a6: 0.000_000, + }, + OxygenLine { + f0: 424.763_020, + a1: 637.700_000, + a2: 0.044_000, + a3: 16.400_000, + a4: 0.000_000, + a5: 0.000_000, + a6: 0.000_000, + }, + OxygenLine { + f0: 487.249_273, + a1: 237.400_000, + a2: 0.049_000, + a3: 16.000_000, + a4: 0.000_000, + a5: 0.000_000, + a6: 0.000_000, + }, + OxygenLine { + f0: 715.392_902, + a1: 98.100_000, + a2: 0.145_000, + a3: 16.000_000, + a4: 0.000_000, + a5: 0.000_000, + a6: 0.000_000, + }, + OxygenLine { + f0: 773.839_490, + a1: 572.300_000, + a2: 0.141_000, + a3: 16.200_000, + a4: 0.000_000, + a5: 0.000_000, + a6: 0.000_000, + }, + OxygenLine { + f0: 834.145_546, + a1: 183.100_000, + a2: 0.145_000, + a3: 14.700_000, + a4: 0.000_000, + a5: 0.000_000, + a6: 0.000_000, + }, +]; + +/// The water vapour absorption lines of table 2 of ITU-R P.676-12. +const WATER_LINES: &[WaterLine] = &[ + WaterLine { + f0: 22.235_080, + b1: 0.107_900, + b2: 2.144_000, + b3: 26.380_000, + b4: 0.760_000, + b5: 5.087_000, + b6: 1.000_000, + }, + WaterLine { + f0: 67.803_960, + b1: 0.001_100, + b2: 8.732_000, + b3: 28.580_000, + b4: 0.690_000, + b5: 4.930_000, + b6: 0.820_000, + }, + WaterLine { + f0: 119.995_940, + b1: 0.000_700, + b2: 8.353_000, + b3: 29.480_000, + b4: 0.700_000, + b5: 4.780_000, + b6: 0.790_000, + }, + WaterLine { + f0: 183.310_087, + b1: 2.273_000, + b2: 0.668_000, + b3: 29.060_000, + b4: 0.770_000, + b5: 5.022_000, + b6: 0.850_000, + }, + WaterLine { + f0: 321.225_630, + b1: 0.047_000, + b2: 6.179_000, + b3: 24.040_000, + b4: 0.670_000, + b5: 4.398_000, + b6: 0.540_000, + }, + WaterLine { + f0: 325.152_888, + b1: 1.514_000, + b2: 1.541_000, + b3: 28.230_000, + b4: 0.640_000, + b5: 4.893_000, + b6: 0.740_000, + }, + WaterLine { + f0: 336.227_764, + b1: 0.001_000, + b2: 9.825_000, + b3: 26.930_000, + b4: 0.690_000, + b5: 4.740_000, + b6: 0.610_000, + }, + WaterLine { + f0: 380.197_353, + b1: 11.670_000, + b2: 1.048_000, + b3: 28.110_000, + b4: 0.540_000, + b5: 5.063_000, + b6: 0.890_000, + }, + WaterLine { + f0: 390.134_508, + b1: 0.004_500, + b2: 7.347_000, + b3: 21.520_000, + b4: 0.630_000, + b5: 4.810_000, + b6: 0.550_000, + }, + WaterLine { + f0: 437.346_667, + b1: 0.063_200, + b2: 5.048_000, + b3: 18.450_000, + b4: 0.600_000, + b5: 4.230_000, + b6: 0.480_000, + }, + WaterLine { + f0: 439.150_807, + b1: 0.909_800, + b2: 3.595_000, + b3: 20.070_000, + b4: 0.630_000, + b5: 4.483_000, + b6: 0.520_000, + }, + WaterLine { + f0: 443.018_343, + b1: 0.192_000, + b2: 5.048_000, + b3: 15.550_000, + b4: 0.600_000, + b5: 5.083_000, + b6: 0.500_000, + }, + WaterLine { + f0: 448.001_085, + b1: 10.410_000, + b2: 1.405_000, + b3: 25.640_000, + b4: 0.660_000, + b5: 5.028_000, + b6: 0.670_000, + }, + WaterLine { + f0: 470.888_999, + b1: 0.325_400, + b2: 3.597_000, + b3: 21.340_000, + b4: 0.660_000, + b5: 4.506_000, + b6: 0.650_000, + }, + WaterLine { + f0: 474.689_092, + b1: 1.260_000, + b2: 2.379_000, + b3: 23.200_000, + b4: 0.650_000, + b5: 4.804_000, + b6: 0.640_000, + }, + WaterLine { + f0: 488.490_108, + b1: 0.252_900, + b2: 2.852_000, + b3: 25.860_000, + b4: 0.690_000, + b5: 5.201_000, + b6: 0.720_000, + }, + WaterLine { + f0: 503.568_532, + b1: 0.037_200, + b2: 6.731_000, + b3: 16.120_000, + b4: 0.610_000, + b5: 3.980_000, + b6: 0.430_000, + }, + WaterLine { + f0: 504.482_692, + b1: 0.012_400, + b2: 6.731_000, + b3: 16.120_000, + b4: 0.610_000, + b5: 4.010_000, + b6: 0.450_000, + }, + WaterLine { + f0: 547.676_440, + b1: 0.978_500, + b2: 0.158_000, + b3: 26.000_000, + b4: 0.700_000, + b5: 4.500_000, + b6: 1.000_000, + }, + WaterLine { + f0: 552.020_960, + b1: 0.184_000, + b2: 0.158_000, + b3: 26.000_000, + b4: 0.700_000, + b5: 4.500_000, + b6: 1.000_000, + }, + WaterLine { + f0: 556.935_985, + b1: 497.000_000, + b2: 0.159_000, + b3: 30.860_000, + b4: 0.690_000, + b5: 4.552_000, + b6: 1.000_000, + }, + WaterLine { + f0: 620.700_807, + b1: 5.015_000, + b2: 2.391_000, + b3: 24.380_000, + b4: 0.710_000, + b5: 4.856_000, + b6: 0.680_000, + }, + WaterLine { + f0: 645.766_085, + b1: 0.006_700, + b2: 8.633_000, + b3: 18.000_000, + b4: 0.600_000, + b5: 4.000_000, + b6: 0.500_000, + }, + WaterLine { + f0: 658.005_280, + b1: 0.273_200, + b2: 7.816_000, + b3: 32.100_000, + b4: 0.690_000, + b5: 4.140_000, + b6: 1.000_000, + }, + WaterLine { + f0: 752.033_113, + b1: 243.400_000, + b2: 0.396_000, + b3: 30.860_000, + b4: 0.680_000, + b5: 4.352_000, + b6: 0.840_000, + }, + WaterLine { + f0: 841.051_732, + b1: 0.013_400, + b2: 8.177_000, + b3: 15.900_000, + b4: 0.330_000, + b5: 5.760_000, + b6: 0.450_000, + }, + WaterLine { + f0: 859.965_698, + b1: 0.132_500, + b2: 8.055_000, + b3: 30.600_000, + b4: 0.680_000, + b5: 4.090_000, + b6: 0.840_000, + }, + WaterLine { + f0: 899.303_175, + b1: 0.054_700, + b2: 7.914_000, + b3: 29.850_000, + b4: 0.680_000, + b5: 4.530_000, + b6: 0.900_000, + }, + WaterLine { + f0: 902.611_085, + b1: 0.038_600, + b2: 8.429_000, + b3: 28.650_000, + b4: 0.700_000, + b5: 5.100_000, + b6: 0.950_000, + }, + WaterLine { + f0: 906.205_957, + b1: 0.183_600, + b2: 5.110_000, + b3: 24.080_000, + b4: 0.700_000, + b5: 4.700_000, + b6: 0.530_000, + }, + WaterLine { + f0: 916.171_582, + b1: 8.400_000, + b2: 1.441_000, + b3: 26.730_000, + b4: 0.700_000, + b5: 5.150_000, + b6: 0.780_000, + }, + WaterLine { + f0: 923.112_692, + b1: 0.007_900, + b2: 10.293_000, + b3: 29.000_000, + b4: 0.700_000, + b5: 5.000_000, + b6: 0.800_000, + }, + WaterLine { + f0: 970.315_022, + b1: 9.009_000, + b2: 1.919_000, + b3: 25.500_000, + b4: 0.640_000, + b5: 4.940_000, + b6: 0.670_000, + }, + WaterLine { + f0: 987.926_764, + b1: 134.600_000, + b2: 0.257_000, + b3: 29.850_000, + b4: 0.680_000, + b5: 4.550_000, + b6: 0.900_000, + }, + WaterLine { + f0: 1_780.000_000, + b1: 17_506.000_000, + b2: 0.952_000, + b3: 196.300_000, + b4: 2.000_000, + b5: 24.150_000, + b6: 5.000_000, + }, +]; + +/// The atmospheric state that the gaseous attenuation depends on. +/// +/// The pressure is the dry air pressure, as the recommendation specifies. The +/// water vapour density is separate because the partial pressure of water +/// vapour is derived from it. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct Atmosphere { + /// The dry air pressure, in hectopascals. + pub pressure_hpa: f64, + /// The water vapour density, in grams per cubic meter. + pub water_vapour_density: f64, + /// The temperature, in kelvin. + pub temperature_k: f64, +} + +impl Atmosphere { + /// The ITU-R reference atmosphere at mean sea level. + /// + /// The values are the ones the recommendation uses for its reference + /// curves: 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. + pub const REFERENCE: Self = Self { + pressure_hpa: 1_013.25, + water_vapour_density: 7.5, + temperature_k: 288.15, + }; + + /// Creates an atmosphere from its three parameters. + pub const fn new(pressure_hpa: f64, water_vapour_density: f64, temperature_k: f64) -> Self { Self { - oxygen_db, - water_vapour_db, + pressure_hpa, + water_vapour_density, + temperature_k, } } - /// Returns the total zenith attenuation, in decibels. - pub fn zenith_db(self) -> f64 { - self.oxygen_db + self.water_vapour_db - } - - /// Returns the total attenuation at an elevation angle, in decibels. + /// Returns the partial pressure of water vapour, in hectopascals. /// - /// The zenith value is scaled by the air mass of the slant path. The - /// result is infinite at or below the horizon. - pub fn at_elevation(self, elevation: f64) -> Decibel { - Decibel::from_db(self.zenith_db() * air_mass(elevation)) + /// The relation is `e = rho T / 216.7`, which follows from the ideal gas + /// law with the molar masses of water and dry air. + pub fn water_vapour_partial_pressure(self) -> f64 { + self.water_vapour_density * self.temperature_k / 216.7 } - /// Returns `true` when both components are finite and not negative. + /// Returns `true` when the state describes a physical atmosphere. pub fn is_valid(self) -> bool { - self.oxygen_db.is_finite() - && self.oxygen_db >= 0.0 - && self.water_vapour_db.is_finite() - && self.water_vapour_db >= 0.0 + self.pressure_hpa.is_finite() + && self.pressure_hpa > 0.0 + && self.water_vapour_density.is_finite() + && self.water_vapour_density >= 0.0 + && self.temperature_k.is_finite() + && self.temperature_k > 0.0 } } -/// Returns the zenith attenuation of oxygen at a frequency, in decibels. +/// Returns the specific attenuation of dry air, in decibels per kilometer. /// -/// The model is a gross simplification of the oxygen absorption complex, -/// which is a set of lines near 60 GHz and a single line near 118 GHz. It -/// captures the shape that matters for a link budget: negligible attenuation -/// at the common satellite bands below 20 GHz, a rapid rise into the 60 GHz -/// complex, and a fall on the far side of it. +/// This is the oxygen contribution, computed by the line-by-line sum of +/// equation 1 of ITU-R P.676-12, plus the non-resonant Debye term that +/// accounts for the far wings of all the lines together. /// -/// The value is the attenuation through the whole atmosphere at the zenith, -/// around 4 dB in the 60 GHz complex. Outside the strong absorption region -/// it is a small fraction of a decibel. +/// The frequency is in gigahertz. /// -/// Returns `None` when the frequency is not positive and finite. -/// -/// Source: the shape follows ITU-R Recommendation P.676, *Attenuation by -/// atmospheric gases*, which tabulates the line by line result this function -/// approximates. -pub fn oxygen_zenith_attenuation(frequency: f64) -> Option { - if !frequency.is_finite() || frequency <= 0.0 { +/// Returns `None` when the frequency or atmosphere is not usable. +pub fn specific_attenuation_dry_air(frequency_ghz: f64, atmosphere: Atmosphere) -> Option { + if !frequency_ghz.is_finite() || frequency_ghz <= 0.0 || !atmosphere.is_valid() { return None; } - let gigahertz = frequency / 1.0e9; - // A narrow complex centred near 60 GHz, plus a weak line at 118 GHz. - let complex = gaussian(gigahertz, 60.0, 8.0, 15.0); - let line = gaussian(gigahertz, 118.0, 12.0, 2.0); - // A baseline that grows slowly with frequency from the far wings of the - // lines and from non-resonant absorption. - let baseline = 0.006 * gigahertz.powf(1.5).min(1.5); - Some(complex + line + baseline) + let theta = 300.0 / atmosphere.temperature_k; + let e = atmosphere.water_vapour_partial_pressure(); + let p = atmosphere.pressure_hpa; + + let mut sum = 0.0; + for line in OXYGEN_LINES { + // Equations 5 and 6: the line width and the pressure shift. + let width = line.a3 * 1.0e-4 * (p * theta.powf(0.8 - line.a4) + 1.1 * e * theta); + let width = (width * width + 2.25e-6).sqrt(); + let shift = (line.a5 + line.a6 * theta) * 1.0e-4 * (p + e) * theta.powf(0.8); + + // Equation 4: the line shape factor. + let shape = frequency_ghz / line.f0 + * ((width - shift * (line.f0 - frequency_ghz)) + / ((line.f0 - frequency_ghz).powi(2) + width * width) + + (width - shift * (line.f0 + frequency_ghz)) + / ((line.f0 + frequency_ghz).powi(2) + width * width)); + + // Equation 3: the line strength. + let strength = line.a1 * 1.0e-7 * p * theta.powi(3) * E.powf(line.a2 * (1.0 - theta)); + sum += strength * shape; + } + + // The non-resonant Debye spectrum, equation 8. + let d = 5.6e-4 * (p + e) * theta.powf(0.8); + let debye = frequency_ghz + * p + * theta.powi(2) + * (6.14e-5 / (d * (1.0 + (frequency_ghz / d).powi(2))) + + 1.4e-12 * p * theta.powf(1.5) / (1.0 + 1.9e-5 * frequency_ghz.powf(1.5))); + + Some(0.182 * frequency_ghz * (sum + debye)) } -/// Returns the zenith attenuation of water vapour at a frequency, in -/// decibels. +/// Returns the specific attenuation of water vapour, in decibels per +/// kilometer. /// -/// Water vapour has a single absorption line at about 22.235 GHz, which is -/// the dominant feature for satellite links, and a much weaker set of lines -/// above 100 GHz. This function models the 22 GHz line in the shape that -/// ITU-R Recommendation P.676 tabulates for a standard atmosphere. +/// This is the water vapour contribution, computed by the line-by-line sum of +/// equation 1 of ITU-R P.676-12. /// -/// The value is the zenith attenuation through the whole atmosphere, about -/// 0.2 dB at the line for a typical water vapour column. +/// The frequency is in gigahertz. /// -/// Returns `None` when the frequency is not positive and finite. -/// -/// Source: shape from ITU-R Recommendation P.676. -pub fn water_vapour_zenith_attenuation(frequency: f64) -> Option { - if !frequency.is_finite() || frequency <= 0.0 { +/// Returns `None` when the frequency or atmosphere is not usable. +pub fn specific_attenuation_water_vapour( + frequency_ghz: f64, + atmosphere: Atmosphere, +) -> Option { + if !frequency_ghz.is_finite() || frequency_ghz <= 0.0 || !atmosphere.is_valid() { return None; } - let gigahertz = frequency / 1.0e9; - // The 22 GHz line, with a much weaker set of lines above 180 GHz. - let line = gaussian(gigahertz, 22.235, 4.0, 0.20); - let upper = gaussian(gigahertz, 183.0, 20.0, 4.0); - Some(line + upper) + let theta = 300.0 / atmosphere.temperature_k; + let e = atmosphere.water_vapour_partial_pressure(); + let p = atmosphere.pressure_hpa; + + let mut sum = 0.0; + for line in WATER_LINES { + // Equation 7: the line width, which has a different form from the + // oxygen case. + let width = + line.b3 * 1.0e-4 * (p * theta.powf(line.b4) + line.b5 * e * theta.powf(line.b6)); + let width = 0.535 * width + + (0.217 * width * width + 2.131_6e-12 * line.f0 * line.f0 / theta).sqrt(); + + let shape = frequency_ghz / line.f0 + * (width / ((line.f0 - frequency_ghz).powi(2) + width * width) + + width / ((line.f0 + frequency_ghz).powi(2) + width * width)); + + // Equation 3 for the water vapour lines. + let strength = line.b1 * 1.0e-1 * e * theta.powf(3.5) * E.powf(line.b2 * (1.0 - theta)); + sum += strength * shape; + } + + Some(0.182 * frequency_ghz * sum) } -/// Returns the total gaseous zenith attenuation at a frequency, in decibels. +/// Returns the total specific attenuation of the atmosphere, in decibels per +/// kilometer. /// -/// Returns `None` when the frequency is not positive and finite. -pub fn gaseous_attenuation(frequency: f64) -> Option { - Some(GaseousAttenuation::new( - oxygen_zenith_attenuation(frequency)?, - water_vapour_zenith_attenuation(frequency)?, - )) -} - -/// Returns a Gaussian bump of the given height, centre, and width. +/// This is the sum of the dry air and water vapour contributions. /// -/// This is used to shape the absorption features. The width is the half -/// width at which the value falls to `1 / e` of the height. -fn gaussian(frequency: f64, centre: f64, width: f64, height: f64) -> f64 { - let offset = (frequency - centre) / width; - height * (-offset * offset).exp() +/// Returns `None` when the frequency or atmosphere is not usable. +pub fn specific_attenuation(frequency_ghz: f64, atmosphere: Atmosphere) -> Option { + Some( + specific_attenuation_dry_air(frequency_ghz, atmosphere)? + + specific_attenuation_water_vapour(frequency_ghz, atmosphere)?, + ) } diff --git a/crates/rf/apsis-atmosphere/src/gaseous_path.rs b/crates/rf/apsis-atmosphere/src/gaseous_path.rs new file mode 100644 index 0000000..ce9a978 --- /dev/null +++ b/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 { + 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::(); + + 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 { + 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 { + if !elevation.is_finite() || elevation.sin() <= 0.0 { + return None; + } + Some(zenith_attenuation(frequency_ghz, atmosphere)? / elevation.sin()) +} diff --git a/crates/rf/apsis-atmosphere/src/lib.rs b/crates/rf/apsis-atmosphere/src/lib.rs index 014bae4..15c4f3c 100644 --- a/crates/rf/apsis-atmosphere/src/lib.rs +++ b/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 -//! absorption lines near 60 GHz and 118 GHz and water vapour near 22 GHz, so -//! the effect is small below about 10 GHz and grows rapidly above 20 GHz. -//! It depends on the elevation angle through the air mass along the path. +//! Gaseous absorption comes from the absorption lines of oxygen and water +//! vapour. Oxygen has a dense group of lines between 50 and 70 GHz and an +//! isolated line at 118.75 GHz, and water vapour has a strong line at +//! 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. -//! It is negligible below about 4 GHz, grows roughly as a power of the rain -//! rate with an exponent near one, and depends strongly on elevation. +//! Rain attenuation is caused by scattering and absorption in +//! precipitation. It is negligible below about 4 GHz, grows as a power of +//! 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 -//! 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. +//! geometry: a signal at a low elevation crosses more atmosphere. //! //! Sources: //! -//! - Gaseous attenuation: ITU-R Recommendation P.676, *Attenuation by -//! atmospheric gases*. -//! - Rain specific attenuation coefficients: ITU-R Recommendation P.838, +//! - Gaseous specific attenuation and equivalent heights: ITU-R +//! Recommendation P.676-12, *Attenuation by atmospheric gases*. +//! - 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. -//! assert_eq!(air_mass(zenith_to_elevation(0.0)), 1.0); +//! // At the reference atmosphere the 60 GHz oxygen complex is heavily +//! // 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 -//! // the cosecant of thirty degrees is two. -//! let thirty = core::f64::consts::PI / 6.0; -//! assert!((air_mass(thirty) - 2.0).abs() < 1e-9); -//! -//! // Elevation and zenith angle are complementary. -//! assert!((elevation_from_zenith(thirty) + thirty - core::f64::consts::FRAC_PI_2).abs() < 1e-12); +//! // A lower elevation crosses more atmosphere. +//! let zenith = zenith_attenuation(12.0, Atmosphere::REFERENCE).expect("valid"); +//! let low = slant_path_attenuation(12.0, 0.2, Atmosphere::REFERENCE).expect("valid"); +//! assert!(low > zenith); //! ``` #![forbid(unsafe_code)] #![deny(missing_docs)] mod gaseous; +mod gaseous_path; mod geometry; mod rain; pub use gaseous::{ - GaseousAttenuation, gaseous_attenuation, oxygen_zenith_attenuation, - water_vapour_zenith_attenuation, + Atmosphere, specific_attenuation, specific_attenuation_dry_air, + 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, - specific_rain_attenuation, + DEFAULT_RAIN_HEIGHT_KM, RainAttenuation, polarisation_tilt, rain_attenuation, + rain_coefficients, specific_rain_attenuation, }; pub use apsis_rf::Decibel; diff --git a/crates/rf/apsis-atmosphere/src/rain.rs b/crates/rf/apsis-atmosphere/src/rain.rs index c340e6a..6193d8f 100644 --- a/crates/rf/apsis-atmosphere/src/rain.rs +++ b/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 -/// frequency. +/// Evaluates one P.838-3 curve term. /// -/// 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 -/// gamma = k R^alpha -/// ``` +/// The polarisation is an angle in radians between the electric field and +/// 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` -/// depend on frequency, on the polarisation, and weakly on the elevation and -/// the rain temperature. The values returned here are for a frequency in -/// gigahertz, interpolated logarithmically in `k` and linearly in `alpha` -/// between the published points. +/// The result is a pair of the coefficient `k` and the exponent `alpha` of +/// the relation `gamma = k R^alpha`, with the rain rate in millimeters per +/// hour and the result in decibels per kilometer. /// -/// The polarisation is given as an angle in radians between the electric -/// field and the horizontal, so zero is horizontal, a right angle is -/// vertical, and the intermediate values are circular or slant -/// polarisation. The coefficients at an angle follow from the horizontal and -/// vertical values through +/// Returns `None` when the frequency, polarisation, or elevation is not +/// usable. The recommendation is stated for frequencies from 1 GHz to +/// 1000 GHz. /// -/// ```text -/// k = (k_h + k_v + (k_h - k_v) cos^2(angle) cos(2 elevation)) / 2 -/// alpha = (k_h a_h + k_v a_v + (k_h a_h - k_v a_v) cos^2(angle) cos(2 elevation)) / (2 k) -/// ``` -/// -/// Returns `None` when the frequency is outside the range the tables cover, -/// from 1 GHz to 1000 GHz. -/// -/// Source: ITU-R Recommendation P.838, *Specific attenuation model for rain -/// for use in prediction methods*, which tabulates `k` and `alpha` for -/// horizontal and vertical polarisation over this frequency range. +/// Source: ITU-R Recommendation P.838-3, equations 2 and 3. 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 - * (horizontal.0 - + vertical.0 - + (horizontal.0 - vertical.0) * cos_squared * cos_two_elevation); - let numerator = horizontal.0 * horizontal.1 - + vertical.0 * vertical.1 - + (horizontal.0 * horizontal.1 - vertical.0 * vertical.1) * cos_squared * cos_two_elevation; - let alpha = if k > 0.0 { numerator / (2.0 * k) } else { 1.0 }; + let k_horizontal = base_k( + K_HORIZONTAL_TERMS, + K_SLOPE_HORIZONTAL, + K_OFFSET_HORIZONTAL, + frequency_ghz, + ); + let k_vertical = base_k( + 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 -/// rate is negative, or when the coefficients are not available. +/// The rain rate is in millimeters per hour. +/// +/// Returns `None` when the rain rate is negative, or when the coefficients +/// are not available. pub fn specific_rain_attenuation( rain_rate: f64, frequency: f64, @@ -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 -/// path length through it is the layer height divided by the sine of the -/// elevation. A reduction factor is applied because the rain rate is not -/// uniform along the path and the longest paths are the ones where the -/// assumption of uniformity fails most. +/// The rain rate is the value exceeded for the percentage of time of +/// interest, in millimeters per hour. The station height above sea level and +/// the rain layer height are in kilometers; pass [`DEFAULT_RAIN_HEIGHT_KM`] +/// for the latter when no local value is available. /// -/// The reduction factor used here is -/// -/// ```text -/// r = 1 / (1 + L_G / L_0) -/// ``` -/// -/// with `L_G` the geometric path and `L_0` a reference length that depends -/// on the rain rate, which is the form given by ITU-R Recommendation P.618. -/// The height of the rain layer is taken as 3 km, a mid latitude value; a -/// caller working in the tropics should raise it and one in polar regions -/// should lower it. +/// 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 +/// a correction for the curvature of the Earth. Two reduction factors then +/// 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 +/// latitude through the angle at which the path crosses the rain layer. /// /// 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 -/// methods required for the design of Earth-space telecommunication -/// systems*. +/// Source: ITU-R Recommendation P.618-13, equations 1 through 10. pub fn rain_attenuation( rain_rate: f64, frequency: f64, polarisation: f64, elevation: f64, + station_height_km: f64, + rain_height_km: f64, + latitude: f64, ) -> Option { 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 - // that accounts for its non-uniformity. - const RAIN_HEIGHT_KM: f64 = 3.0; - let geometric_path_km = RAIN_HEIGHT_KM * air_mass(elevation); - let reference_length_km = 35.0 * (-0.015 * rain_rate).exp(); - let reduction = 1.0 / (1.0 + geometric_path_km / reference_length_km); - let effective_path_km = geometric_path_km * reduction; + // Equation 2: the slant path length, with the Earth curvature term. + let height_difference = rain_height_km - station_height_km; + if height_difference <= 0.0 { + // The station is at or above the rain layer, so the path is empty. + return Some(RainAttenuation { + specific_db_per_km: specific, + 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, - attenuation_db: specific * effective_path_km, + path_length_km: effective_length, + 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 } diff --git a/crates/rf/apsis-atmosphere/tests/gaseous.rs b/crates/rf/apsis-atmosphere/tests/gaseous.rs index d7599de..feee915 100644 --- a/crates/rf/apsis-atmosphere/tests/gaseous.rs +++ b/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, - water_vapour_zenith_attenuation, + Atmosphere, equivalent_heights, specific_attenuation, specific_attenuation_dry_air, + 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() { - // Below about 10 GHz the oxygen contribution at the zenith is a small - // fraction of a decibel. - for gigahertz in [1.0, 2.0, 4.0, 6.0, 10.0] { - let value = oxygen_zenith_attenuation(gigahertz * 1.0e9).expect("valid frequency"); - assert!(value < 0.2, "{gigahertz} GHz gave {value} dB"); - assert!(value >= 0.0); +fn the_reference_atmosphere_matches_the_standard() { + let atmosphere = reference(); + assert_eq!(atmosphere.pressure_hpa, 1_013.25); + assert_eq!(atmosphere.water_vapour_density, 7.5); + assert_eq!(atmosphere.temperature_k, 288.15); + assert!(atmosphere.is_valid()); +} + +#[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() { - let peak = oxygen_zenith_attenuation(60.0e9).expect("valid"); - // The 60 GHz complex reaches about 15 dB through the whole atmosphere. - assert!((peak - 15.0).abs() < 0.5, "peak {peak}"); - // A frequency on either side is lower. - let below = oxygen_zenith_attenuation(50.0e9).expect("valid"); - let above = oxygen_zenith_attenuation(70.0e9).expect("valid"); - assert!(below < peak); - assert!(above < peak); -} - -#[test] -fn oxygen_attenuation_is_monotonic_through_the_complex() { - // The value rises steadily from below the complex into its centre. - let mut previous = 0.0; - for gigahertz in [40.0, 45.0, 50.0, 55.0, 58.0, 59.0] { - let value = oxygen_zenith_attenuation(gigahertz * 1.0e9).expect("valid"); - assert!(value > previous, "{gigahertz} GHz gave {value}"); - previous = value; +fn water_vapour_attenuation_matches_the_recommendation() { + let cases = [ + (1.0, 0.000_050_904_617), + (10.0, 0.005_974_125_245), + (22.235, 0.178_977_992_373), + (60.0, 0.154_841_840_636), + (183.0, 27.665_008_314_166), + (300.0, 5.221_329_041_105), + ]; + for (frequency, expected) in cases { + let actual = specific_attenuation_water_vapour(frequency, reference()).expect("valid"); + assert!( + (actual - expected).abs() < TOLERANCE, + "{frequency} GHz: got {actual}, expected {expected}" + ); } } #[test] -fn water_vapour_attenuation_peaks_near_twenty_two_gigahertz() { - let peak = water_vapour_zenith_attenuation(22.235e9).expect("valid"); - assert!((peak - 0.20).abs() < 1e-9, "peak {peak}"); - let below = water_vapour_zenith_attenuation(15.0e9).expect("valid"); - let above = water_vapour_zenith_attenuation(30.0e9).expect("valid"); - assert!(below < peak); - assert!(above < peak); +fn the_two_contributions_add() { + let frequency = 22.235; + let dry = specific_attenuation_dry_air(frequency, reference()).expect("valid"); + let wet = specific_attenuation_water_vapour(frequency, reference()).expect("valid"); + let total = specific_attenuation(frequency, reference()).expect("valid"); + assert!((total - (dry + wet)).abs() < 1e-12); } #[test] -fn water_vapour_attenuation_is_small_below_ten_gigahertz() { - for gigahertz in [1.0, 4.0, 6.0, 10.0] { - let value = water_vapour_zenith_attenuation(gigahertz * 1.0e9).expect("valid"); - assert!(value < 0.05, "{gigahertz} GHz gave {value} dB"); +fn the_oxygen_complex_dominates_near_sixty_gigahertz() { + // The 60 GHz complex is a set of more than thirty lines, so the dry air + // term is many orders of magnitude above the water vapour term there. + 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() { - let total = gaseous_attenuation(20.0e9).expect("valid"); - assert!((total.zenith_db() - (total.oxygen_db + total.water_vapour_db)).abs() < TOLERANCE); +fn dry_air_attenuation_peaks_inside_the_oxygen_complex() { + // The strongest line is near 60.4 GHz, and the value there must exceed + // 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() { - let gas = GaseousAttenuation::new(0.05, 0.15); - let zenith = gas.at_elevation(FRAC_PI_2); - let thirty = gas.at_elevation(FRAC_PI_6); - // The zenith total is 0.2 dB and thirty degrees doubles it. - assert!((zenith.db() - 0.2).abs() < TOLERANCE); - assert!((thirty.db() - 0.4).abs() < 1e-9); +fn water_vapour_attenuation_grows_with_the_vapour_density() { + let dry_air = Atmosphere::new(1_013.25, 0.0, 288.15); + let wet_air = Atmosphere::new(1_013.25, 15.0, 288.15); + let dry = specific_attenuation_water_vapour(22.235, dry_air).expect("valid"); + let wet = specific_attenuation_water_vapour(22.235, wet_air).expect("valid"); + assert_eq!(dry, 0.0); + assert!(wet > 0.0); } #[test] -fn attenuation_at_the_horizon_is_infinite() { - let gas = GaseousAttenuation::new(0.05, 0.15); - assert!(gas.at_elevation(0.0).db().is_infinite()); +fn equivalent_heights_match_the_recommendation() { + let cases = [ + (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() { - let none = GaseousAttenuation::new(0.0, 0.0); - assert_eq!(none.at_elevation(FRAC_PI_6).db(), 0.0); +fn the_dry_height_is_capped_below_seventy_gigahertz() { + // Equation 30 limits the dry air height to a pressure dependent bound at + // 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() { - assert!(GaseousAttenuation::default().is_valid()); - assert!(GaseousAttenuation::new(0.05, 0.15).is_valid()); - assert!(!GaseousAttenuation::new(-0.1, 0.15).is_valid()); - assert!(!GaseousAttenuation::new(0.05, f64::NAN).is_valid()); +fn zenith_attenuation_matches_the_recommendation() { + let cases = [ + (1.0, 0.029_998_832_95), + (12.0, 0.058_856_728_55), + (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() { - assert_eq!(oxygen_zenith_attenuation(0.0), None); - assert_eq!(oxygen_zenith_attenuation(-1.0), None); - assert_eq!(oxygen_zenith_attenuation(f64::NAN), None); - assert_eq!(water_vapour_zenith_attenuation(0.0), None); - assert_eq!(water_vapour_zenith_attenuation(f64::INFINITY), None); - assert_eq!(gaseous_attenuation(0.0), None); +fn invalid_inputs_are_rejected() { + assert_eq!(specific_attenuation(0.0, reference()), None); + assert_eq!(specific_attenuation(-1.0, reference()), None); + assert_eq!(specific_attenuation(f64::NAN, reference()), None); + assert_eq!( + specific_attenuation_dry_air(12.0, Atmosphere::new(0.0, 7.5, 288.15)), + 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); } diff --git a/crates/rf/apsis-atmosphere/tests/rain.rs b/crates/rf/apsis-atmosphere/tests/rain.rs index 80eb0cb..62ff8f8 100644 --- a/crates/rf/apsis-atmosphere/tests/rain.rs +++ b/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; + +/// The relative tolerance on the coefficient values. +const TOLERANCE: f64 = 1e-9; + +fn assert_relative(actual: f64, expected: f64, context: &str) { + let relative = (actual - expected).abs() / expected.abs(); + assert!( + relative < TOLERANCE, + "{context}: got {actual}, expected {expected}, relative {relative}" + ); +} #[test] -fn coefficients_are_available_across_the_tabulated_range() { - for gigahertz in [1.0, 4.0, 12.0, 20.0, 30.0, 100.0, 500.0, 1_000.0] { - let coefficients = - rain_coefficients(gigahertz * 1.0e9, HORIZONTAL, ELEVATION).expect("valid"); - assert!( - coefficients.0 > 0.0, - "{gigahertz} GHz gave k {}", - coefficients.0 - ); - assert!( - coefficients.1 > 0.5 && coefficients.1 < 1.5, - "{gigahertz} GHz gave alpha {}", - coefficients.1 - ); +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() { - // At 12 GHz the table gives k = 0.0168 for horizontal polarisation, - // with alpha_h = 1.200 and alpha_v = 1.000. The combined alpha depends - // on the elevation through the `cos(2 elevation)` term, so at thirty - // degrees it is the weighted mean of the two, which is 1.15. - let (k, alpha) = rain_coefficients(12.0e9, HORIZONTAL, ELEVATION).expect("valid"); - assert!((k - 0.016_8).abs() < 1e-9, "k {k}"); - assert!((alpha - 1.15).abs() < 1e-9, "alpha {alpha}"); - - // At the zenith the `cos(2 elevation)` factor is minus one, which - // selects the vertical alpha instead, so the pair swaps behaviour. - let (_, zenith_alpha) = rain_coefficients(12.0e9, HORIZONTAL, FRAC_PI_2).expect("valid"); - assert!((zenith_alpha - 1.0).abs() < 1e-9, "alpha {zenith_alpha}"); +fn coefficients_match_the_recommendation_at_thirty_degrees() { + let cases = [ + (12.0, 0.023_944_109_789_7, 1.174_670_745_2), + (20.0, 0.092_201_236_241_3, 1.047_387_560_19), + (30.0, 0.238_905_952_257, 0.944_198_048_818), + ]; + for (frequency, k_expected, alpha_expected) in cases { + let (k, alpha) = + rain_coefficients(frequency * 1.0e9, HORIZONTAL, THIRTY_DEGREES).expect("valid"); + assert_relative(k, k_expected, &format!("{frequency} GHz k")); + assert_relative(alpha, alpha_expected, &format!("{frequency} GHz alpha")); + } } #[test] -fn horizontal_coefficients_exceed_vertical_at_mid_frequencies() { - // At 20 GHz the table gives k_h = 0.0751 against k_v = 0.0691. - let (horizontal, _) = rain_coefficients(20.0e9, HORIZONTAL, ELEVATION).expect("valid"); - let (vertical, _) = rain_coefficients(20.0e9, VERTICAL, ELEVATION).expect("valid"); - assert!(horizontal > vertical, "{horizontal} vs {vertical}"); +fn horizontal_and_vertical_differ_off_the_zenith() { + // Away from the zenith the polarisation matters. At 12 GHz and low + // elevation the vertical coefficient is the larger of the two, and the + // two converge as the elevation rises and the geometry becomes + // 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() { - // At 11 GHz the value lies between the 10 and 12 GHz entries in - // logarithmic k, so it is the geometric mean of the two. - let (k, _) = rain_coefficients(11.0e9, HORIZONTAL, ELEVATION).expect("valid"); - let expected = (0.008_87_f64 * 0.016_8).sqrt(); - assert!((k - expected).abs() < 1e-9, "k {k} expected {expected}"); +fn circular_polarisation_lies_between_the_linear_cases() { + // A tilt of forty five degrees is circular polarisation, whose + // coefficients are the mean of the two linear cases. + let horizontal = rain_coefficients(12.0e9, HORIZONTAL, THIRTY_DEGREES).expect("valid"); + let vertical = rain_coefficients(12.0e9, VERTICAL, THIRTY_DEGREES).expect("valid"); + 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() { - assert_eq!(rain_coefficients(0.5e9, HORIZONTAL, ELEVATION), None); - assert_eq!(rain_coefficients(2_000.0e9, HORIZONTAL, ELEVATION), None); - assert_eq!(rain_coefficients(f64::NAN, HORIZONTAL, ELEVATION), None); - assert_eq!(rain_coefficients(12.0e9, f64::NAN, ELEVATION), None); +fn the_coefficients_are_smooth_across_the_band() { + // The curve fits are continuous, so a small change in frequency gives a + // small change in the coefficients. + let mut previous = rain_coefficients(10.0e9, HORIZONTAL, ZENITH).expect("valid"); + 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"); - // The value is k R^alpha with the coefficients for this frequency, - // polarisation, and elevation. - let (k, alpha) = rain_coefficients(12.0e9, HORIZONTAL, ELEVATION).expect("valid"); - assert!((specific - k * rain_rate.powf(alpha)).abs() < 1e-9); - // For 20 mm/h at 12 GHz and thirty degrees this is about 0.53 dB/km. - assert!((specific - 0.526_6).abs() < 0.001, "specific {specific}"); + specific_rain_attenuation(rain_rate, frequency, HORIZONTAL, THIRTY_DEGREES).expect("valid"); + let (k, alpha) = rain_coefficients(frequency, HORIZONTAL, THIRTY_DEGREES).expect("valid"); + assert_relative(specific, k * rain_rate.powf(alpha), "specific attenuation"); + // At 12 GHz, thirty degrees elevation, and 20 mm/h this is 0.808 dB/km. + assert!( + (specific - 0.808_128_772_7).abs() < 1e-6, + "specific {specific}" + ); } #[test] fn a_higher_rain_rate_gives_more_attenuation() { - let light = specific_rain_attenuation(1.0, 12.0e9, HORIZONTAL, ELEVATION).expect("valid"); - let heavy = specific_rain_attenuation(50.0, 12.0e9, HORIZONTAL, ELEVATION).expect("valid"); + let light = specific_rain_attenuation(1.0, 12.0e9, HORIZONTAL, THIRTY_DEGREES).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 - // the coefficients give. - let (_, alpha) = rain_coefficients(12.0e9, HORIZONTAL, ELEVATION).expect("valid"); - assert!((heavy / light - 50.0_f64.powf(alpha)).abs() < 1e-6); + let (_, alpha) = rain_coefficients(12.0e9, HORIZONTAL, THIRTY_DEGREES).expect("valid"); + assert_relative(heavy / light, 50.0_f64.powf(alpha), "rate ratio"); } #[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"); - assert!( - (result.attenuation_db - result.specific_db_per_km * result.path_length_km).abs() < 1e-9 + let result = rain_attenuation( + 20.0, + 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() { - // The reduction factor represents the fact that heavy rain does not fill - // the whole slant path, so the effective length is below the geometric - // length through a three kilometre layer. - let result = rain_attenuation(20.0, 12.0e9, HORIZONTAL, ELEVATION).expect("valid"); - let geometric = 3.0 / ELEVATION.sin(); - assert!((geometric - 6.0).abs() < 1e-6); +fn the_reduction_factors_are_physical() { + // The horizontal reduction factor accounts for rain not filling the + // whole horizontal extent of the path, so it is at most one. The + // vertical adjustment is the ratio of two path estimates and can exceed + // one when the path crosses the rain layer at a shallow angle, which is + // the case Recommendation P.618 allows for. + 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.path_length_km + result.horizontal_reduction > 0.0 && result.horizontal_reduction <= 1.0, + "horizontal reduction {}", + 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 low = rain_attenuation(20.0, 12.0e9, HORIZONTAL, 10.0_f64.to_radians()).expect("valid"); + let high = rain_attenuation( + 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 ku_band = rain_attenuation(20.0, 12.0e9, HORIZONTAL, ELEVATION).expect("valid"); + let c_band = rain_attenuation( + 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!(rain_attenuation(20.0, 12.0e9, HORIZONTAL, -0.1), None); + assert_eq!( + 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!((polarisation_tilt(45.0) - FRAC_PI_2 / 2.0).abs() < 1e-12); + assert_relative(polarisation_tilt(90.0), FRAC_PI_2, "ninety degrees"); + assert_relative( + polarisation_tilt(45.0), + FRAC_PI_2 / 2.0, + "forty five degrees", + ); } diff --git a/validation/pixi.lock b/validation/pixi.lock index 673f380..34e8bfb 100644 --- a/validation/pixi.lock +++ b/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 diff --git a/validation/pixi.toml b/validation/pixi.toml index e93a083..e718682 100644 --- a/validation/pixi.toml +++ b/validation/pixi.toml @@ -15,5 +15,6 @@ numpy = "*" sgp4 = "*" astropy = "*" skyfield = "*" +itur = ">=0.4.0, <0.5" [tasks] diff --git a/validation/python/atmosphere/check_gas.py b/validation/python/atmosphere/check_gas.py new file mode 100644 index 0000000..8a53b14 --- /dev/null +++ b/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)) diff --git a/validation/python/atmosphere/check_rain.py b/validation/python/atmosphere/check_rain.py new file mode 100644 index 0000000..f554ce4 --- /dev/null +++ b/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))