From 82bf248ab77818a8f7622598dab9af9b59467475 Mon Sep 17 00:00:00 2001
From: Graham McIntire <graham@mcintire.me>
Date: Fri, 10 Apr 2026 08:31:15 -0500
Subject: [PATCH] Phase 4: Ray-traced duct geometry
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

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).
---
 lib/microwaveprop/backtest/features.ex       |  56 +++++++
 lib/microwaveprop/propagation/duct.ex        | 153 +++++++++++++++++
 lib/mix/tasks/hrrr_native_derive.ex          |  52 +++---
 test/microwaveprop/propagation/duct_test.exs | 167 +++++++++++++++++++
 4 files changed, 408 insertions(+), 20 deletions(-)
 create mode 100644 lib/microwaveprop/propagation/duct.ex
 create mode 100644 test/microwaveprop/propagation/duct_test.exs

diff --git a/lib/microwaveprop/backtest/features.ex b/lib/microwaveprop/backtest/features.ex
index 0903d88c..9b4f30a6 100644
--- a/lib/microwaveprop/backtest/features.ex
+++ b/lib/microwaveprop/backtest/features.ex
@@ -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
diff --git a/lib/microwaveprop/propagation/duct.ex b/lib/microwaveprop/propagation/duct.ex
new file mode 100644
index 00000000..922bff18
--- /dev/null
+++ b/lib/microwaveprop/propagation/duct.ex
@@ -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
diff --git a/lib/mix/tasks/hrrr_native_derive.ex b/lib/mix/tasks/hrrr_native_derive.ex
index 518a797f..f5646146 100644
--- a/lib/mix/tasks/hrrr_native_derive.ex
+++ b/lib/mix/tasks/hrrr_native_derive.ex
@@ -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
diff --git a/test/microwaveprop/propagation/duct_test.exs b/test/microwaveprop/propagation/duct_test.exs
new file mode 100644
index 00000000..5c680422
--- /dev/null
+++ b/test/microwaveprop/propagation/duct_test.exs
@@ -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