prop/lib/microwaveprop/propagation/duct.ex

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.
  """
  @spec refractivity_profile(map()) :: [{float(), float()}]
  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.
  """
  @spec m_profile([{float(), float()}]) :: [{float(), float()}]
  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.
  """
  @spec detect_ducts([{float(), float()}]) ::
          [%{base_m: float(), top_m: float(), thickness_m: float(), m_deficit: float()}]
  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}], acc ->
      process_m_pair(acc, h1, m1, h2, m2)
    end)
    |> close_trailing_duct()
  end

  defp process_m_pair({ducts, current_duct}, h1, m1, h2, m2) when m2 < m1 do
    # M is decreasing — we're in a duct
    duct = extend_or_start_duct(current_duct, h1, m1, h2, m2)
    {ducts, duct}
  end

  defp process_m_pair({ducts, nil}, _h1, _m1, _h2, _m2), do: {ducts, nil}

  defp process_m_pair({ducts, duct}, _h1, _m1, _h2, _m2) do
    {[finalize_duct(duct) | ducts], nil}
  end

  defp extend_or_start_duct(nil, h1, m1, h2, m2) do
    %{base_m: h1, base_m_val: m1, top_m: h2, min_m_val: m2}
  end

  defp extend_or_start_duct(duct, _h1, _m1, h2, m2) do
    %{duct | top_m: h2, min_m_val: min(duct.min_m_val, m2)}
  end

  defp close_trailing_duct({ducts, nil}), do: Enum.reverse(ducts)
  defp close_trailing_duct({ducts, duct}), do: Enum.reverse([finalize_duct(duct) | ducts])

  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.
  """
  @spec min_trapped_frequency_ghz(%{thickness_m: float(), m_deficit: float()}) :: float()
  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.
  """
  @spec analyze(map()) ::
          %{
            ducts: [
              %{
                base_m: float(),
                top_m: float(),
                thickness_m: float(),
                m_deficit: float(),
                min_freq_ghz: float()
              }
            ],
            best_duct_band_ghz: float() | nil
          }
  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