Phase 4: Ray-traced duct geometry
Duct module (Propagation.Duct): - refractivity_profile/1: ITU-R P.453 N at each native level - m_profile/1: modified refractivity M = N + 157*h(km) - detect_ducts/1: find contiguous regions where dM/dh < 0, returning base/top height, thickness, and M-deficit per duct - min_trapped_frequency_ghz/1: waveguide approximation (Bean & Dutton) for the minimum frequency a duct of given geometry can trap - analyze/1: full pipeline from native profile to duct list + best trapped frequency across all ducts Derive task updated to also compute ducts JSONB and best_duct_band_ghz alongside the Phase 2 turbulence fields. Backtest features: duct_thickness, best_duct_freq, duct_usable_10ghz, duct_usable_24ghz, duct_usable_47ghz. Real-data validation: 2022-08-20 12Z TX profile shows 0 ducts (M increases monotonically) — correct for a well-mixed boundary layer on a turbulent August afternoon (Ri=0.16).
This commit is contained in:
parent
e2daf5a81d
commit
82bf248ab7
4 changed files with 408 additions and 20 deletions
|
|
@ -159,6 +159,62 @@ defmodule Microwaveprop.Backtest.Features do
|
|||
end
|
||||
end
|
||||
|
||||
@doc """
|
||||
Max duct thickness (m) from the nearest native profile.
|
||||
|
||||
Larger ducts trap lower frequencies and support longer-range propagation.
|
||||
Returns nil if no ducts detected or no native profile available.
|
||||
"""
|
||||
@spec duct_thickness(float, float, DateTime.t()) :: float | nil
|
||||
def duct_thickness(lat, lon, valid_time) do
|
||||
with %HrrrNativeProfile{ducts: ducts} when is_list(ducts) and ducts != [] <-
|
||||
find_nearest_native(lat, lon, valid_time) do
|
||||
ducts |> Enum.map(& &1["thickness_m"]) |> Enum.max(fn -> nil end)
|
||||
else
|
||||
_ -> nil
|
||||
end
|
||||
end
|
||||
|
||||
@doc """
|
||||
Lowest trapped frequency (GHz) across all ducts in the nearest native profile.
|
||||
|
||||
A lower value means stronger ducting that can trap more bands.
|
||||
"""
|
||||
@spec best_duct_freq(float, float, DateTime.t()) :: float | nil
|
||||
def best_duct_freq(lat, lon, valid_time) do
|
||||
with %HrrrNativeProfile{best_duct_band_ghz: f} when is_float(f) <-
|
||||
find_nearest_native(lat, lon, valid_time) do
|
||||
f
|
||||
else
|
||||
_ -> nil
|
||||
end
|
||||
end
|
||||
|
||||
@doc """
|
||||
Boolean (1.0 or 0.0): is there a duct that can trap the given band?
|
||||
|
||||
Pass the band frequency in GHz. Returns 1.0 if any duct has
|
||||
min_freq_ghz <= band_ghz, 0.0 otherwise, nil if no data.
|
||||
"""
|
||||
@spec duct_usable_for_band(float, float, DateTime.t(), float) :: float | nil
|
||||
def duct_usable_for_band(lat, lon, valid_time, band_ghz) do
|
||||
with %HrrrNativeProfile{best_duct_band_ghz: f} when is_float(f) <-
|
||||
find_nearest_native(lat, lon, valid_time) do
|
||||
if f <= band_ghz, do: 1.0, else: 0.0
|
||||
else
|
||||
_ -> nil
|
||||
end
|
||||
end
|
||||
|
||||
@doc "Duct usable for 10 GHz."
|
||||
def duct_usable_10ghz(lat, lon, vt), do: duct_usable_for_band(lat, lon, vt, 10.0)
|
||||
|
||||
@doc "Duct usable for 24 GHz."
|
||||
def duct_usable_24ghz(lat, lon, vt), do: duct_usable_for_band(lat, lon, vt, 24.0)
|
||||
|
||||
@doc "Duct usable for 47 GHz."
|
||||
def duct_usable_47ghz(lat, lon, vt), do: duct_usable_for_band(lat, lon, vt, 47.0)
|
||||
|
||||
defp find_nearest_native(lat, lon, valid_time) do
|
||||
dlat = 0.07
|
||||
dlon = 0.07
|
||||
|
|
|
|||
153
lib/microwaveprop/propagation/duct.ex
Normal file
153
lib/microwaveprop/propagation/duct.ex
Normal file
|
|
@ -0,0 +1,153 @@
|
|||
defmodule Microwaveprop.Propagation.Duct do
|
||||
@moduledoc """
|
||||
Detects atmospheric ducts from a native-level HRRR profile by
|
||||
computing the modified refractivity M-profile and finding regions
|
||||
where dM/dh < 0.
|
||||
|
||||
This replaces the scalar `min_refractivity_gradient` with physical
|
||||
duct geometry: base height, top height, thickness, M-deficit, and
|
||||
the minimum frequency the duct can trap. The per-band trapped
|
||||
frequency is the key improvement — a 50 m duct can trap 24 GHz but
|
||||
not 3 GHz, and the old scalar approach had no way to express that.
|
||||
|
||||
References:
|
||||
- ITU-R P.453-14: refractivity formula N = 77.6*P/T + 3.73e5*e/T²
|
||||
- ITU-R P.834-9: ducting defined by regions where dM/dh < 0
|
||||
- Bean & Dutton (1966): minimum trapped wavelength λ_min ≈ ...
|
||||
"""
|
||||
|
||||
@doc """
|
||||
Compute the refractivity N at each level of a native profile.
|
||||
|
||||
Returns `[{height_m, N}, ...]` sorted ascending by height.
|
||||
Uses ITU-R P.453 with water vapor pressure derived from specific
|
||||
humidity.
|
||||
"""
|
||||
def refractivity_profile(%{heights_m: heights, temp_k: temps, spfh: spfhs, pressure_pa: pressures}) do
|
||||
[heights, temps, spfhs, pressures]
|
||||
|> Enum.zip()
|
||||
|> Enum.map(fn {h, t, q, p} ->
|
||||
# Water vapor pressure: e = q*P / (0.622 + 0.378*q) [Pa]
|
||||
e = q * p / (0.622 + 0.378 * q)
|
||||
# Convert to hPa for the P.453 formula
|
||||
p_hpa = p / 100.0
|
||||
e_hpa = e / 100.0
|
||||
|
||||
n = 77.6 * p_hpa / t + 3.73e5 * e_hpa / (t * t)
|
||||
{h, n}
|
||||
end)
|
||||
end
|
||||
|
||||
@doc """
|
||||
Convert an N-profile to the modified refractivity M-profile.
|
||||
|
||||
M = N + 157 * h (where h is in km)
|
||||
|
||||
In a standard atmosphere M always increases with height. A duct
|
||||
exists wherever M decreases.
|
||||
"""
|
||||
def m_profile(n_profile) do
|
||||
Enum.map(n_profile, fn {h, n} ->
|
||||
{h, n + 157.0 * h / 1000.0}
|
||||
end)
|
||||
end
|
||||
|
||||
@doc """
|
||||
Detect ducts from an M-profile.
|
||||
|
||||
A duct is a contiguous region where M decreases with height (dM/dh < 0).
|
||||
Returns a list of `%{base_m, top_m, thickness_m, m_deficit}` maps,
|
||||
where `m_deficit` is the total M decrease (positive number) across the
|
||||
duct — the strength of the trapping.
|
||||
"""
|
||||
def detect_ducts(m_profile) when length(m_profile) < 2, do: []
|
||||
|
||||
def detect_ducts(m_profile) do
|
||||
m_profile
|
||||
|> Enum.chunk_every(2, 1, :discard)
|
||||
|> Enum.reduce({[], nil}, fn [{h1, m1}, {h2, m2}], {ducts, current_duct} ->
|
||||
if m2 < m1 do
|
||||
# M is decreasing — we're in a duct
|
||||
case current_duct do
|
||||
nil ->
|
||||
{ducts, %{base_m: h1, base_m_val: m1, top_m: h2, min_m_val: m2}}
|
||||
|
||||
duct ->
|
||||
{ducts, %{duct | top_m: h2, min_m_val: min(duct.min_m_val, m2)}}
|
||||
end
|
||||
else
|
||||
# M is increasing — close any open duct
|
||||
case current_duct do
|
||||
nil ->
|
||||
{ducts, nil}
|
||||
|
||||
duct ->
|
||||
finished = finalize_duct(duct)
|
||||
{[finished | ducts], nil}
|
||||
end
|
||||
end
|
||||
end)
|
||||
|> then(fn {ducts, current} ->
|
||||
case current do
|
||||
nil -> Enum.reverse(ducts)
|
||||
duct -> Enum.reverse([finalize_duct(duct) | ducts])
|
||||
end
|
||||
end)
|
||||
end
|
||||
|
||||
defp finalize_duct(%{base_m: base, top_m: top, base_m_val: m_base, min_m_val: m_min}) do
|
||||
%{
|
||||
base_m: base,
|
||||
top_m: top,
|
||||
thickness_m: top - base,
|
||||
m_deficit: m_base - m_min
|
||||
}
|
||||
end
|
||||
|
||||
@doc """
|
||||
Minimum trapped frequency (GHz) for a duct of given thickness and M-deficit.
|
||||
|
||||
Uses the waveguide approximation from Bean & Dutton:
|
||||
λ_max = 2.5 * d * sqrt(ΔM * 1e-6) (meters)
|
||||
f_min = c / λ_max
|
||||
|
||||
where d is duct thickness in meters, ΔM is the M-deficit.
|
||||
Returns a very high frequency (999 GHz) for degenerate ducts.
|
||||
"""
|
||||
def min_trapped_frequency_ghz(%{thickness_m: d, m_deficit: delta_m})
|
||||
when d > 0 and delta_m > 0 do
|
||||
# Maximum trapped wavelength (meters)
|
||||
lambda_max = 2.5 * d * :math.sqrt(delta_m * 1.0e-6)
|
||||
|
||||
if lambda_max > 1.0e-6 do
|
||||
# f = c / λ, convert to GHz
|
||||
3.0e8 / lambda_max / 1.0e9
|
||||
else
|
||||
999.0
|
||||
end
|
||||
end
|
||||
|
||||
def min_trapped_frequency_ghz(_), do: 999.0
|
||||
|
||||
@doc """
|
||||
Full analysis pipeline: native profile → ducts with trapped frequencies.
|
||||
|
||||
Returns `%{ducts: [duct_map], best_duct_band_ghz: float | nil}`
|
||||
where each duct_map includes `:min_freq_ghz` and
|
||||
`best_duct_band_ghz` is the lowest min_freq across all detected ducts.
|
||||
"""
|
||||
def analyze(profile) do
|
||||
ducts =
|
||||
profile
|
||||
|> refractivity_profile()
|
||||
|> m_profile()
|
||||
|> detect_ducts()
|
||||
|> Enum.map(fn duct ->
|
||||
Map.put(duct, :min_freq_ghz, min_trapped_frequency_ghz(duct))
|
||||
end)
|
||||
|
||||
best = ducts |> Enum.map(& &1.min_freq_ghz) |> Enum.min(fn -> nil end)
|
||||
|
||||
%{ducts: ducts, best_duct_band_ghz: best}
|
||||
end
|
||||
end
|
||||
|
|
@ -14,6 +14,7 @@ defmodule Mix.Tasks.HrrrNativeDeriveFields do
|
|||
|
||||
import Ecto.Query
|
||||
|
||||
alias Microwaveprop.Propagation.Duct
|
||||
alias Microwaveprop.Propagation.Inversion
|
||||
alias Microwaveprop.Repo
|
||||
alias Microwaveprop.Weather.HrrrNativeProfile
|
||||
|
|
@ -65,27 +66,38 @@ defmodule Mix.Tasks.HrrrNativeDeriveFields do
|
|||
level_count: profile.level_count
|
||||
}
|
||||
|
||||
case Inversion.find_inversion_top(p) do
|
||||
{:ok, %{height_m: top_h, level_idx: top_idx, base_idx: base_idx}} ->
|
||||
ri = Inversion.bulk_richardson(p, base_idx, top_idx)
|
||||
shear = Inversion.shear_magnitude(p, base_idx, top_idx)
|
||||
theta_e_jump = ThetaE.theta_e_jump(p, base_idx, top_idx)
|
||||
# Phase 2: turbulence features
|
||||
inversion_fields =
|
||||
case Inversion.find_inversion_top(p) do
|
||||
{:ok, %{height_m: top_h, level_idx: top_idx, base_idx: base_idx}} ->
|
||||
ri = Inversion.bulk_richardson(p, base_idx, top_idx)
|
||||
shear = Inversion.shear_magnitude(p, base_idx, top_idx)
|
||||
theta_e_jump = ThetaE.theta_e_jump(p, base_idx, top_idx)
|
||||
|
||||
[
|
||||
inversion_top_m: top_h,
|
||||
bulk_richardson: ri,
|
||||
shear_at_top_ms: shear,
|
||||
theta_e_jump_k: theta_e_jump
|
||||
]
|
||||
[
|
||||
inversion_top_m: top_h,
|
||||
bulk_richardson: ri,
|
||||
shear_at_top_ms: shear,
|
||||
theta_e_jump_k: theta_e_jump
|
||||
]
|
||||
|
||||
:none ->
|
||||
# No inversion — still fill nulls so we know the row was processed
|
||||
[
|
||||
inversion_top_m: nil,
|
||||
bulk_richardson: nil,
|
||||
shear_at_top_ms: nil,
|
||||
theta_e_jump_k: nil
|
||||
]
|
||||
end
|
||||
:none ->
|
||||
[
|
||||
inversion_top_m: nil,
|
||||
bulk_richardson: nil,
|
||||
shear_at_top_ms: nil,
|
||||
theta_e_jump_k: nil
|
||||
]
|
||||
end
|
||||
|
||||
# Phase 4: duct geometry
|
||||
duct_result = Duct.analyze(p)
|
||||
|
||||
duct_fields = [
|
||||
ducts: duct_result.ducts,
|
||||
best_duct_band_ghz: duct_result.best_duct_band_ghz
|
||||
]
|
||||
|
||||
inversion_fields ++ duct_fields
|
||||
end
|
||||
end
|
||||
|
|
|
|||
167
test/microwaveprop/propagation/duct_test.exs
Normal file
167
test/microwaveprop/propagation/duct_test.exs
Normal file
|
|
@ -0,0 +1,167 @@
|
|||
defmodule Microwaveprop.Propagation.DuctTest do
|
||||
use ExUnit.Case, async: true
|
||||
|
||||
alias Microwaveprop.Propagation.Duct
|
||||
|
||||
defp profile(levels) do
|
||||
%{
|
||||
level_count: length(levels),
|
||||
heights_m: Enum.map(levels, &elem(&1, 0)),
|
||||
temp_k: Enum.map(levels, &elem(&1, 1)),
|
||||
spfh: Enum.map(levels, &elem(&1, 2)),
|
||||
pressure_pa: Enum.map(levels, &elem(&1, 3))
|
||||
}
|
||||
end
|
||||
|
||||
describe "refractivity_profile/1" do
|
||||
test "computes N at each level with plausible values" do
|
||||
p =
|
||||
profile([
|
||||
{10.0, 295.0, 0.012, 101_000.0},
|
||||
{500.0, 290.0, 0.006, 95_000.0},
|
||||
{1000.0, 283.0, 0.003, 90_000.0}
|
||||
])
|
||||
|
||||
n_profile = Duct.refractivity_profile(p)
|
||||
assert length(n_profile) == 3
|
||||
|
||||
# N should decrease with height in a normal atmosphere
|
||||
[{_, n0}, {_, n1}, {_, n2}] = n_profile
|
||||
assert n0 > n1
|
||||
assert n1 > n2
|
||||
# Typical surface N = 300-400
|
||||
assert n0 > 250 and n0 < 500
|
||||
end
|
||||
end
|
||||
|
||||
describe "m_profile/1" do
|
||||
test "M = N + 157*h where h is in km" do
|
||||
n_profile = [{0.0, 350.0}, {1000.0, 280.0}]
|
||||
m_profile = Duct.m_profile(n_profile)
|
||||
|
||||
[{h0, m0}, {1000.0, m1}] = m_profile
|
||||
assert h0 == 0.0
|
||||
assert_in_delta m0, 350.0, 0.01
|
||||
# M at 1km = 280 + 157*1.0 = 437
|
||||
assert_in_delta m1, 437.0, 0.01
|
||||
end
|
||||
end
|
||||
|
||||
describe "detect_ducts/1" do
|
||||
test "returns empty list when M increases monotonically (no duct)" do
|
||||
# Standard atmosphere: M always increases with height
|
||||
m_profile = [
|
||||
{0.0, 350.0},
|
||||
{100.0, 365.0},
|
||||
{300.0, 397.0},
|
||||
{500.0, 428.0},
|
||||
{1000.0, 507.0}
|
||||
]
|
||||
|
||||
assert Duct.detect_ducts(m_profile) == []
|
||||
end
|
||||
|
||||
test "detects a surface-based duct where M decreases" do
|
||||
# Surface-based duct: M decreases from 0 to 100m, then resumes increasing
|
||||
m_profile = [
|
||||
{0.0, 380.0},
|
||||
{50.0, 370.0},
|
||||
{100.0, 360.0},
|
||||
{200.0, 375.0},
|
||||
{500.0, 428.0}
|
||||
]
|
||||
|
||||
ducts = Duct.detect_ducts(m_profile)
|
||||
assert length(ducts) == 1
|
||||
|
||||
[duct] = ducts
|
||||
assert duct.base_m == 0.0
|
||||
assert duct.top_m == 100.0
|
||||
assert duct.thickness_m == 100.0
|
||||
assert duct.m_deficit > 0
|
||||
end
|
||||
|
||||
test "detects an elevated duct" do
|
||||
# Elevated duct between 500-700m
|
||||
m_profile = [
|
||||
{0.0, 350.0},
|
||||
{100.0, 366.0},
|
||||
{300.0, 397.0},
|
||||
{500.0, 430.0},
|
||||
{600.0, 425.0},
|
||||
{700.0, 420.0},
|
||||
{1000.0, 500.0}
|
||||
]
|
||||
|
||||
ducts = Duct.detect_ducts(m_profile)
|
||||
assert length(ducts) == 1
|
||||
|
||||
[duct] = ducts
|
||||
assert duct.base_m == 500.0
|
||||
assert duct.top_m == 700.0
|
||||
assert duct.thickness_m == 200.0
|
||||
end
|
||||
|
||||
test "detects multiple ducts" do
|
||||
m_profile = [
|
||||
{0.0, 380.0},
|
||||
{50.0, 370.0},
|
||||
{100.0, 390.0},
|
||||
{500.0, 430.0},
|
||||
{600.0, 425.0},
|
||||
{700.0, 440.0}
|
||||
]
|
||||
|
||||
ducts = Duct.detect_ducts(m_profile)
|
||||
assert length(ducts) == 2
|
||||
end
|
||||
end
|
||||
|
||||
describe "min_trapped_frequency_ghz/1" do
|
||||
test "a narrow duct with weak deficit traps only high frequencies" do
|
||||
# 10m duct, 1 M-unit deficit → λ_max ≈ 0.025m → f ≈ 12 GHz
|
||||
duct = %{thickness_m: 10.0, m_deficit: 1.0}
|
||||
f = Duct.min_trapped_frequency_ghz(duct)
|
||||
assert_in_delta f, 12.0, 1.0
|
||||
end
|
||||
|
||||
test "a moderate duct traps sub-GHz frequencies" do
|
||||
# 50m duct, 10 M-units → λ_max ≈ 0.395m → f ≈ 0.76 GHz
|
||||
duct = %{thickness_m: 50.0, m_deficit: 10.0}
|
||||
f = Duct.min_trapped_frequency_ghz(duct)
|
||||
assert_in_delta f, 0.76, 0.1
|
||||
end
|
||||
|
||||
test "a thick strong duct traps very low frequencies" do
|
||||
# 200m, 20 M-units → λ_max ≈ 2.24m → f ≈ 0.13 GHz
|
||||
duct = %{thickness_m: 200.0, m_deficit: 20.0}
|
||||
f = Duct.min_trapped_frequency_ghz(duct)
|
||||
assert f < 0.2
|
||||
end
|
||||
|
||||
test "thicker ducts trap lower frequencies than thinner ones" do
|
||||
thin = Duct.min_trapped_frequency_ghz(%{thickness_m: 30.0, m_deficit: 5.0})
|
||||
thick = Duct.min_trapped_frequency_ghz(%{thickness_m: 150.0, m_deficit: 5.0})
|
||||
assert thick < thin
|
||||
end
|
||||
end
|
||||
|
||||
describe "analyze/1" do
|
||||
test "full pipeline from native profile to duct list with frequencies" do
|
||||
# Profile with a surface-based duct (strong inversion)
|
||||
p =
|
||||
profile([
|
||||
{10.0, 300.0, 0.018, 101_325.0},
|
||||
{50.0, 301.0, 0.017, 100_800.0},
|
||||
{100.0, 302.0, 0.012, 100_300.0},
|
||||
{150.0, 298.0, 0.008, 99_800.0},
|
||||
{300.0, 293.0, 0.005, 97_000.0},
|
||||
{500.0, 288.0, 0.003, 95_000.0}
|
||||
])
|
||||
|
||||
result = Duct.analyze(p)
|
||||
assert is_list(result.ducts)
|
||||
assert is_float(result.best_duct_band_ghz) or is_nil(result.best_duct_band_ghz)
|
||||
end
|
||||
end
|
||||
end
|
||||
Loading…
Add table
Reference in a new issue