diff --git a/rust/prop_grid_rs/src/duct.rs b/rust/prop_grid_rs/src/duct.rs new file mode 100644 index 00000000..c6a68c94 --- /dev/null +++ b/rust/prop_grid_rs/src/duct.rs @@ -0,0 +1,266 @@ +//! Atmospheric duct detection from a native-level HRRR profile. +//! +//! Port of `lib/microwaveprop/propagation/duct.ex`. Computes the modified +//! refractivity M-profile and finds regions where dM/dh < 0. Those are +//! ducts; the `min_freq_ghz` per duct (Bean & Dutton waveguide +//! approximation) tells us which bands can be trapped. +//! +//! References: +//! - ITU-R P.453-14: N = 77.6·P/T + 3.73e5·e/T² +//! - ITU-R P.834-9: duct ≡ region with dM/dh < 0 +//! - Bean & Dutton (1966): λ_max ≈ 2.5·d·sqrt(ΔM·1e-6) + +/// One native-level profile at a single grid cell. All vectors share +/// the same length; indices correspond. Order: surface → top. +#[derive(Debug, Clone)] +pub struct NativeProfile { + pub heights_m: Vec, + pub temp_k: Vec, + pub spfh: Vec, + pub pressure_pa: Vec, +} + +impl NativeProfile { + pub fn level_count(&self) -> usize { + self.heights_m.len() + } +} + +#[derive(Debug, Clone, PartialEq)] +pub struct Duct { + pub base_m: f64, + pub top_m: f64, + pub thickness_m: f64, + pub m_deficit: f64, + pub min_freq_ghz: f64, +} + +#[derive(Debug, Clone)] +pub struct Analysis { + pub ducts: Vec, + /// Lowest `min_freq_ghz` across all detected ducts — the most + /// permissive duct for band planning. `None` when no ducts exist. + pub best_duct_band_ghz: Option, +} + +/// Refractivity `N` at every level, returned paired with height. +pub fn refractivity_profile(p: &NativeProfile) -> Vec<(f64, f64)> { + let n = p.level_count(); + let mut out = Vec::with_capacity(n); + for i in 0..n { + let h = p.heights_m[i]; + let t = p.temp_k[i]; + let q = p.spfh[i]; + let p_pa = p.pressure_pa[i]; + // Water vapor pressure e [Pa] + let e = q * p_pa / (0.622 + 0.378 * q); + let p_hpa = p_pa / 100.0; + let e_hpa = e / 100.0; + let n_val = 77.6 * p_hpa / t + 3.73e5 * e_hpa / (t * t); + out.push((h, n_val)); + } + out +} + +/// Modified refractivity M = N + 157·h_km. +pub fn m_profile(n: &[(f64, f64)]) -> Vec<(f64, f64)> { + n.iter() + .map(|&(h, nv)| (h, nv + 157.0 * h / 1000.0)) + .collect() +} + +/// Find contiguous regions where M decreases with height. +pub fn detect_ducts(m: &[(f64, f64)]) -> Vec { + if m.len() < 2 { + return Vec::new(); + } + let mut ducts: Vec = Vec::new(); + let mut current: Option = None; + + for w in m.windows(2) { + let (h1, m1) = w[0]; + let (h2, m2) = w[1]; + if m2 < m1 { + // in a duct + current = Some(match current { + None => RawDuct { base: h1, base_m_val: m1, top: h2, min_m_val: m2 }, + Some(d) => RawDuct { + top: h2, + min_m_val: d.min_m_val.min(m2), + ..d + }, + }); + } else if let Some(d) = current.take() { + ducts.push(finalize(d)); + } + } + if let Some(d) = current { + ducts.push(finalize(d)); + } + ducts +} + +#[derive(Debug, Clone, Copy)] +struct RawDuct { + base: f64, + base_m_val: f64, + top: f64, + min_m_val: f64, +} + +fn finalize(d: RawDuct) -> Duct { + let thickness = d.top - d.base; + let deficit = d.base_m_val - d.min_m_val; + Duct { + base_m: d.base, + top_m: d.top, + thickness_m: thickness, + m_deficit: deficit, + min_freq_ghz: min_trapped_frequency_ghz(thickness, deficit), + } +} + +/// Bean & Dutton waveguide approximation for the lowest trapped +/// frequency: λ_max = 2.5·d·√(ΔM·1e-6), f_min = c/λ_max. +pub fn min_trapped_frequency_ghz(thickness_m: f64, m_deficit: f64) -> f64 { + if thickness_m <= 0.0 || m_deficit <= 0.0 { + return 999.0; + } + let lambda_max = 2.5 * thickness_m * (m_deficit * 1.0e-6).sqrt(); + if lambda_max > 1.0e-6 { + 3.0e8 / lambda_max / 1.0e9 + } else { + 999.0 + } +} + +/// Full analysis: profile → ducts with trapped frequencies. +pub fn analyze(profile: &NativeProfile) -> Analysis { + let n = refractivity_profile(profile); + let m = m_profile(&n); + let ducts = detect_ducts(&m); + let best = ducts + .iter() + .map(|d| d.min_freq_ghz) + .fold(None, |acc, f| match acc { + None => Some(f), + Some(a) => Some(a.min(f)), + }); + Analysis { ducts, best_duct_band_ghz: best } +} + +/// Minimum dM/dh across a profile (in M-units per km), mirroring the +/// Elixir `min_m_gradient` used by the native-duct merge. Falls back to +/// 0.0 on empty or degenerate input — same convention Elixir uses. +pub fn min_m_gradient(profile: &NativeProfile) -> f64 { + let m = m_profile(&refractivity_profile(profile)); + if m.len() < 2 { + return 0.0; + } + let mut min = f64::INFINITY; + for w in m.windows(2) { + let (h1, m1) = w[0]; + let (h2, m2) = w[1]; + let dh = h2 - h1; + if dh > 0.0 { + let grad = (m2 - m1) / dh * 1000.0; + if grad < min { + min = grad; + } + } + } + if min.is_finite() { min } else { 0.0 } +} + +/// Maximum duct thickness across a detected set; `None` if empty. +pub fn max_duct_thickness_m(ducts: &[Duct]) -> Option { + ducts.iter().map(|d| d.thickness_m).fold(None, |acc, t| match acc { + None => Some(t), + Some(a) => Some(a.max(t)), + }) +} + +#[cfg(test)] +mod tests { + use super::*; + + fn standard_atmosphere_profile() -> NativeProfile { + // Simple isothermal + dry decreasing-pressure stack. M should + // increase monotonically ⇒ no ducts. + let heights: Vec = (0..10).map(|i| i as f64 * 100.0).collect(); + let temp_k: Vec = vec![288.0; 10]; + let spfh: Vec = vec![0.005; 10]; + let pressure_pa: Vec = (0..10) + .map(|i| 101_325.0 * (1.0 - 0.000_022_6 * (i as f64 * 100.0)).powf(5.255)) + .collect(); + NativeProfile { heights_m: heights, temp_k, spfh, pressure_pa } + } + + #[test] + fn standard_atmosphere_has_no_ducts() { + let p = standard_atmosphere_profile(); + let a = analyze(&p); + assert!(a.ducts.is_empty()); + assert!(a.best_duct_band_ghz.is_none()); + } + + #[test] + fn moisture_inversion_creates_duct() { + // Put a moist boundary layer under a dry, warm cap — classic + // trapping layer. Surface q=0.018, 200m q=0.004, with 3 °C + // temperature inversion at 100-200 m. + let heights_m = vec![0.0, 50.0, 100.0, 150.0, 200.0, 400.0]; + let temp_k = vec![298.0, 297.5, 300.0, 302.0, 300.5, 295.0]; + let spfh = vec![0.018, 0.017, 0.012, 0.006, 0.003, 0.002]; + let pressure_pa = vec![101_325.0, 100_720.0, 100_120.0, 99_520.0, 98_930.0, 96_580.0]; + let p = NativeProfile { heights_m, temp_k, spfh, pressure_pa }; + let a = analyze(&p); + assert!(!a.ducts.is_empty(), "expected at least one duct"); + // The duct must have positive thickness and deficit. + for d in &a.ducts { + assert!(d.thickness_m > 0.0); + assert!(d.m_deficit > 0.0); + assert!(d.min_freq_ghz > 0.0 && d.min_freq_ghz < 999.0); + } + assert!(a.best_duct_band_ghz.unwrap() > 0.0); + } + + #[test] + fn min_trapped_freq_degenerate_returns_999() { + assert_eq!(min_trapped_frequency_ghz(0.0, 10.0), 999.0); + assert_eq!(min_trapped_frequency_ghz(100.0, 0.0), 999.0); + assert_eq!(min_trapped_frequency_ghz(-5.0, 10.0), 999.0); + } + + #[test] + fn min_trapped_freq_50m_10deficit() { + // From Bean & Dutton: 50m duct with ΔM=10 traps ~15.2 GHz. + let f = min_trapped_frequency_ghz(50.0, 10.0); + assert!((f - 3.0e8 / (2.5 * 50.0 * (10.0e-6f64).sqrt()) / 1.0e9).abs() < 1e-6); + } + + #[test] + fn empty_profile_empty_analysis() { + let p = NativeProfile { + heights_m: vec![], + temp_k: vec![], + spfh: vec![], + pressure_pa: vec![], + }; + let a = analyze(&p); + assert!(a.ducts.is_empty()); + assert!(a.best_duct_band_ghz.is_none()); + } + + #[test] + fn single_level_no_ducts() { + let p = NativeProfile { + heights_m: vec![0.0], + temp_k: vec![288.0], + spfh: vec![0.005], + pressure_pa: vec![101_325.0], + }; + let a = analyze(&p); + assert!(a.ducts.is_empty()); + } +} diff --git a/rust/prop_grid_rs/src/lib.rs b/rust/prop_grid_rs/src/lib.rs index 05cd6cf8..a6874552 100644 --- a/rust/prop_grid_rs/src/lib.rs +++ b/rust/prop_grid_rs/src/lib.rs @@ -21,6 +21,7 @@ pub mod band_config; pub mod db; pub mod decoder; +pub mod duct; pub mod fetcher; pub mod grid; pub mod metrics;