feat(prop-grid-rs): port Propagation.Duct to Rust

Pure port of lib/microwaveprop/propagation/duct.ex — the ITU-R P.834
duct detector that computes the modified refractivity M-profile and
finds regions where dM/dh < 0, then uses Bean & Dutton to derive the
minimum trapped frequency per duct.

Part of Stream A (f00 port): compute_duct_metrics and the per-cell
native-level duct merge both build on analyze(). Landing it now with
standalone unit tests lets later tasks (fetch_native_duct_grid, the
wgrib2 subprocess for native levels) treat it as a known-good
dependency instead of a concurrent port.

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<f64>,
+    pub temp_k: Vec<f64>,
+    pub spfh: Vec<f64>,
+    pub pressure_pa: Vec<f64>,
+}
+
+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<Duct>,
+    /// 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<f64>,
+}
+
+/// 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<Duct> {
+    if m.len() < 2 {
+        return Vec::new();
+    }
+    let mut ducts: Vec<Duct> = Vec::new();
+    let mut current: Option<RawDuct> = 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<f64> {
+    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<f64> = (0..10).map(|i| i as f64 * 100.0).collect();
+        let temp_k: Vec<f64> = vec![288.0; 10];
+        let spfh: Vec<f64> = vec![0.005; 10];
+        let pressure_pa: Vec<f64> = (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());
+    }
+}

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;