//! Minimum-refractivity-gradient derivation from pressure-level //! profiles. 1:1 port of the subset of //! `lib/microwaveprop/weather/sounding_params.ex` needed by the //! propagation grid scorer: `min_refractivity_gradient` per cell. //! //! The full Elixir module derives CAPE/LI/K-index/ducts — those are //! used for contact enrichment and f00's native-duct merge, neither of //! which is in the Rust f01..f18 scope. #[derive(Debug, Clone)] pub struct Level { pub pres_mb: f64, pub hght_m: f64, pub tmpc: f64, pub dwpc: Option, } /// Buck equation for saturation vapor pressure (hPa) at temperature in °C. pub fn sat_vap_pres(t_c: f64) -> f64 { 6.1121 * ((18.678 - t_c / 234.5) * (t_c / (257.14 + t_c))).exp() } /// Minimum refractivity gradient (dN/dh, N-units per km) from a pressure /// profile. Returns `None` if fewer than 3 usable levels or adjacent /// pairs too close in height. Levels are internally sorted descending by /// pressure (surface first) to match the Elixir orientation. pub fn min_refractivity_gradient(mut levels: Vec) -> Option { levels.retain(|l| l.dwpc.is_some()); if levels.len() < 3 { return None; } levels.sort_by(|a, b| { b.pres_mb .partial_cmp(&a.pres_mb) .unwrap_or(std::cmp::Ordering::Equal) }); let sfc_hght = levels[0].hght_m; let refract: Vec<(f64, f64)> = levels .iter() .map(|l| { let t_k = l.tmpc + 273.15; let e = sat_vap_pres(l.dwpc.expect("filtered above")); let n = 77.6 * l.pres_mb / t_k + 3.73e5 * e / (t_k * t_k); let h_agl = l.hght_m - sfc_hght; (h_agl, n) }) .collect(); let mut min_grad: Option = None; for pair in refract.windows(2) { let (h0, n0) = pair[0]; let (h1, n1) = pair[1]; let dh_km = (h1 - h0) / 1000.0; if dh_km.abs() < 0.01 { continue; } let g = (n1 - n0) / dh_km; min_grad = Some(match min_grad { Some(prev) if prev <= g => prev, _ => g, }); } min_grad } #[cfg(test)] mod tests { use super::*; #[test] fn buck_sat_vap_at_freezing() { // At 0 °C, Buck gives ~6.112 hPa. let e = sat_vap_pres(0.0); assert!((e - 6.112).abs() < 0.02, "{e}"); } #[test] fn nil_for_too_few_levels() { let levels = vec![Level { pres_mb: 1000.0, hght_m: 0.0, tmpc: 25.0, dwpc: Some(20.0), }]; assert!(min_refractivity_gradient(levels).is_none()); } #[test] fn typical_tropospheric_gradient_is_negative() { // Synthetic profile: refractivity decreases with altitude → negative // dN/dh. Realistic numbers from a standard atmosphere. let levels = vec![ Level { pres_mb: 1000.0, hght_m: 100.0, tmpc: 25.0, dwpc: Some(20.0), }, Level { pres_mb: 925.0, hght_m: 800.0, tmpc: 20.0, dwpc: Some(14.0), }, Level { pres_mb: 850.0, hght_m: 1500.0, tmpc: 15.0, dwpc: Some(8.0), }, Level { pres_mb: 700.0, hght_m: 3000.0, tmpc: 5.0, dwpc: Some(-5.0), }, ]; let g = min_refractivity_gradient(levels).unwrap(); assert!(g < 0.0, "gradient should be negative: {g}"); assert!(g > -500.0, "no duct expected: {g}"); } #[test] fn surface_temperature_inversion_yields_stronger_gradient() { // Surface inversion: cold moist near the ground, warm dry above — // strong negative refractivity gradient. let inversion = vec![ Level { pres_mb: 1000.0, hght_m: 100.0, tmpc: 18.0, dwpc: Some(17.0), }, Level { pres_mb: 990.0, hght_m: 300.0, tmpc: 22.0, dwpc: Some(5.0), }, Level { pres_mb: 950.0, hght_m: 700.0, tmpc: 20.0, dwpc: Some(3.0), }, Level { pres_mb: 850.0, hght_m: 1500.0, tmpc: 15.0, dwpc: Some(0.0), }, ]; let strong = min_refractivity_gradient(inversion).unwrap(); let neutral = vec![ Level { pres_mb: 1000.0, hght_m: 100.0, tmpc: 25.0, dwpc: Some(18.0), }, Level { pres_mb: 925.0, hght_m: 800.0, tmpc: 20.0, dwpc: Some(14.0), }, Level { pres_mb: 850.0, hght_m: 1500.0, tmpc: 15.0, dwpc: Some(8.0), }, ]; let baseline = min_refractivity_gradient(neutral).unwrap(); assert!( strong < baseline, "inversion {strong} vs baseline {baseline}" ); } }