//! 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());
    }
}