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

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

This commit is contained in:
ApsisRS contributors 2026-09-26 15:13:48 +08:00
parent fb1dd47a9f
commit d5ba27e784
13 changed files with 2466 additions and 424 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
crates/apsis-rs/Cargo.toml
View File

@ -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"]

42
crates/apsis-rs/examples/dump_gas.rs Normal file
View File

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

36
crates/apsis-rs/examples/dump_rain.rs Normal file
View File

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

1002
crates/rf/apsis-atmosphere/src/gaseous.rs
View File

File diff suppressed because it is too large Load Diff

170
crates/rf/apsis-atmosphere/src/gaseous_path.rs Normal file
View File

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

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

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

483
crates/rf/apsis-atmosphere/src/rain.rs
View File

@ -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<RainAttenuation> {
if !rain_rate.is_finite() || rain_rate < 0.0 {
return None;
@ -142,103 +348,80 @@ pub fn rain_attenuation(
if !elevation.is_finite() || elevation.sin() <= 0.0 {
return None;
}
if !station_height_km.is_finite() || !rain_height_km.is_finite() {
return None;
}
if !latitude.is_finite() {
return None;
}
let specific = specific_rain_attenuation(rain_rate, frequency, polarisation, elevation)?;
// The geometric path through the rain layer, and the reduction factor
// 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
}

288
crates/rf/apsis-atmosphere/tests/gaseous.rs
View File

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

388
crates/rf/apsis-atmosphere/tests/rain.rs
View File

@ -1,105 +1,152 @@
//! Tests for rain attenuation.
//!
//! The coefficient values are those of ITU-R Recommendation P.838-3,
//! produced with a reference implementation of the recommendation and
//! recorded here. The path results exercise the P.618-13 method.
use apsis_atmosphere::{
polarisation_tilt, rain_attenuation, rain_coefficients, specific_rain_attenuation,
DEFAULT_RAIN_HEIGHT_KM, polarisation_tilt, rain_attenuation, rain_coefficients,
specific_rain_attenuation,
};
use core::f64::consts::FRAC_PI_2;
use core::f64::consts::{FRAC_PI_2, PI};
/// Horizontal polarisation.
const HORIZONTAL: f64 = 0.0;
/// Vertical polarisation.
const VERTICAL: f64 = FRAC_PI_2;
/// The zenith, where the two polarisations coincide.
const ZENITH: f64 = FRAC_PI_2;
/// A mid latitude elevation of thirty degrees.
const ELEVATION: f64 = core::f64::consts::PI / 6.0;
const THIRTY_DEGREES: f64 = PI / 6.0;
#[test]
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");
/// 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!(
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
relative < TOLERANCE,
"{context}: got {actual}, expected {expected}, relative {relative}"
);
}
#[test]
fn coefficients_match_the_recommendation_at_the_zenith() {
// At the zenith the elevation term vanishes, so horizontal and vertical
// polarisation give the same coefficients.
let cases = [
(12.0, 0.024_203_061_156_1, 1.151_599_196_3),
(20.0, 0.093_876_937_766_6, 1.019_877_631_17),
(30.0, 0.234_699_253_968, 0.931_114_875_787),
];
for (frequency, k_expected, alpha_expected) in cases {
for polarisation in [HORIZONTAL, VERTICAL] {
let (k, alpha) =
rain_coefficients(frequency * 1.0e9, polarisation, ZENITH).expect("valid");
assert_relative(k, k_expected, &format!("{frequency} GHz k"));
assert_relative(alpha, alpha_expected, &format!("{frequency} GHz alpha"));
}
}
}
#[test]
fn coefficients_at_a_tabulated_frequency_match_the_table() {
// 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",
);
}

177
validation/pixi.lock
View File

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

1
validation/pixi.toml
View File

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

136
validation/python/atmosphere/check_gas.py Normal file
View File

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

93
validation/python/atmosphere/check_rain.py Normal file
View File

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