prop/rust/prop_grid_rs/src/sounding_params.rs

//! 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<f64>,
}

/// 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<Level>) -> Option<f64> {
    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<f64> = 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}");
    }
}