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()
}

/// Surface refractivity N = 77.6 * p / T_K + 3.73e5 * e / T_K² for the
/// lowest-altitude level that has a dewpoint. None if no usable level.
pub fn surface_refractivity(levels: &[Level]) -> Option<f64> {
    let usable: Vec<&Level> = levels.iter().filter(|l| l.dwpc.is_some()).collect();
    if usable.is_empty() {
        return None;
    }
    // Highest pressure = lowest altitude = "surface" in this profile.
    let sfc = usable
        .iter()
        .max_by(|a, b| {
            a.pres_mb
                .partial_cmp(&b.pres_mb)
                .unwrap_or(std::cmp::Ordering::Equal)
        })
        .copied()?;
    let t_k = sfc.tmpc + 273.15;
    let e = sat_vap_pres(sfc.dwpc?);
    Some(77.6 * sfc.pres_mb / t_k + 3.73e5 * e / (t_k * t_k))
}

/// A duct (trapping layer) is present when any part of the refractivity
/// profile has dN/dh < -157 N/km. We use the *minimum* gradient as a
/// cheap proxy — matches Elixir's `ducting_detected` computed from
/// `SoundingParams.derive/1`.
pub fn ducting_detected(min_gradient: Option<f64>) -> bool {
    matches!(min_gradient, Some(g) if g < -157.0)
}

/// 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;
        // Only zero-thickness layers are skipped. The previous 10 m
        // guard discarded the very thin native-level layers HRRR uses
        // to resolve sharp surface-based inversions and ducts.
        if dh_km.abs() < 1.0e-9 {
            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}"
        );
    }
}