293 lines
8.5 KiB
Rust
293 lines
8.5 KiB
Rust
//! 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());
|
|
}
|
|
}
|