prop/lib/microwaveprop/weather/weather_layers.ex

defmodule Microwaveprop.Weather.WeatherLayers do
  @moduledoc "Derives scalar weather map layers from HRRR profile data."

  alias Microwaveprop.Propagation.Duct
  alias Microwaveprop.Weather.SoundingParams

  @doc """
  Derives a map of scalar weather fields from a row containing HRRR profile data.

  Expects a map with keys: `:surface_temp_c`, `:surface_dewpoint_c`,
  `:surface_pressure_mb`, `:surface_refractivity`, `:profile`, `:duct_characteristics`.

  Returns a map with derived fields: `:surface_rh`, `:surface_refractivity`,
  `:temp_850mb`, `:dewpoint_850mb`, `:temp_700mb`, `:dewpoint_700mb`,
  `:lapse_rate`, `:mid_lapse_rate`, `:inversion_strength`,
  `:inversion_base_m`, `:duct_base_m`, `:duct_strength`.
  """
  @spec derive(map()) :: map()
  def derive(row) do
    sorted = sort_profile(row.profile)

    %{
      surface_rh: compute_rh(row.surface_temp_c, row.surface_dewpoint_c),
      surface_refractivity: row.surface_refractivity,
      temp_850mb: level_value(sorted, 850, "tmpc"),
      dewpoint_850mb: level_value(sorted, 850, "dwpc"),
      temp_700mb: level_value(sorted, 700, "tmpc"),
      dewpoint_700mb: level_value(sorted, 700, "dwpc"),
      lapse_rate: compute_lapse_rate(sorted),
      mid_lapse_rate: compute_layer_lapse_rate(sorted, 850, 700),
      inversion_strength: inversion_strength(sorted),
      inversion_base_m: inversion_base_m(sorted),
      duct_base_m: duct_field(row.duct_characteristics, "base"),
      duct_strength: duct_field(row.duct_characteristics, "strength"),
      duct_cutoff_ghz: duct_cutoff_ghz(row.duct_characteristics)
    }
  end

  # Lowest microwave frequency that any detected duct can trap — the
  # "best operationally useful band" for the cell. Prefers the
  # pre-computed `:min_freq_ghz` on native-shape ducts (string or atom
  # key); falls back to `Duct.min_trapped_frequency_ghz/1` when only
  # `thickness_m` + `m_deficit` are present; skips sounding-shape ducts
  # that carry neither (these don't appear on the HRRR weather grid in
  # practice but we tolerate them for safety).
  defp duct_cutoff_ghz(nil), do: nil
  defp duct_cutoff_ghz([]), do: nil

  defp duct_cutoff_ghz(ducts) when is_list(ducts) do
    ducts
    |> Enum.map(&duct_min_freq/1)
    |> Enum.reject(&is_nil/1)
    |> Enum.min(fn -> nil end)
  end

  defp duct_min_freq(duct) when is_map(duct) do
    d = stringify_keys(duct)

    case d["min_freq_ghz"] do
      f when is_number(f) and f > 0 -> f * 1.0
      _ -> duct_min_freq_from_geometry(d)
    end
  end

  defp duct_min_freq(_), do: nil

  defp duct_min_freq_from_geometry(duct) do
    thickness = duct["thickness_m"]
    deficit = duct["m_deficit"]

    if valid_geometry?(thickness, deficit) do
      Duct.min_trapped_frequency_ghz(%{thickness_m: thickness * 1.0, m_deficit: deficit * 1.0})
    end
  end

  defp valid_geometry?(t, d) when is_number(t) and is_number(d) and t > 0 and d > 0, do: true
  defp valid_geometry?(_, _), do: false

  defp sort_profile(nil), do: []
  defp sort_profile([]), do: []

  defp sort_profile(profile) do
    # Accept either legacy string keys (`"pres"`, `"hght"`) or the
    # Rust ProfilesFile shape (`"pres_mb"` / atom `:pres_mb` post-atomization
    # via ProfilesFile.read). `SoundingParams.normalize_profile_entry/1`
    # is the single source of truth for that coercion.
    profile
    |> Enum.map(&SoundingParams.normalize_profile_entry/1)
    |> Enum.filter(fn p -> p["pres"] != nil and p["tmpc"] != nil and p["hght"] != nil end)
    |> Enum.sort_by(fn p -> p["pres"] end, :desc)
  end

  defp compute_rh(t_c, td_c) do
    100.0 * SoundingParams.sat_vap_pres(td_c) / SoundingParams.sat_vap_pres(t_c)
  end

  defp level_value([], _target_pres, _key), do: nil

  defp level_value(sorted, target_pres, key) do
    nearest = Enum.min_by(sorted, fn p -> abs(p["pres"] - target_pres) end)

    if abs(nearest["pres"] - target_pres) <= 25.0 do
      nearest[key]
    end
  end

  defp compute_lapse_rate([]), do: nil
  defp compute_lapse_rate([_]), do: nil

  defp compute_lapse_rate(sorted) do
    surface = hd(sorted)
    top = List.last(sorted)
    dh_km = (top["hght"] - surface["hght"]) / 1000.0

    if dh_km > 0 do
      (surface["tmpc"] - top["tmpc"]) / dh_km
    end
  end

  # Lapse rate over a specific pressure layer (e.g. 850→700 mb), used as
  # a cap-mechanism diagnostic. A steep lapse rate above the boundary
  # layer indicates an elevated mixed layer (EML), which is the classic
  # cap structure: warm boundary layer → very stable inversion → steep
  # dry-adiabatic layer aloft. Returns nil if either bracketing level
  # is missing from the profile.
  defp compute_layer_lapse_rate(sorted, lower_pres, upper_pres) when lower_pres > upper_pres do
    lower = level_entry(sorted, lower_pres)
    upper = level_entry(sorted, upper_pres)

    if lower == nil or upper == nil do
      nil
    else
      dh_km = (upper["hght"] - lower["hght"]) / 1000.0
      if dh_km > 0, do: (lower["tmpc"] - upper["tmpc"]) / dh_km
    end
  end

  defp level_entry([], _target_pres), do: nil

  defp level_entry(sorted, target_pres) do
    nearest = Enum.min_by(sorted, fn p -> abs(p["pres"] - target_pres) end)
    if abs(nearest["pres"] - target_pres) <= 25.0, do: nearest
  end

  defp inversion_strength([]), do: 0.0

  defp inversion_strength(sorted) do
    sorted
    |> Enum.chunk_every(2, 1, :discard)
    |> Enum.reduce(0.0, fn [lower, upper], max_str ->
      dt = upper["tmpc"] - lower["tmpc"]

      if dt > max_str do
        dt
      else
        max_str
      end
    end)
  end

  defp inversion_base_m([]), do: nil

  defp inversion_base_m(sorted) do
    sfc_hght = hd(sorted)["hght"]

    result =
      sorted
      |> Enum.chunk_every(2, 1, :discard)
      |> Enum.reduce({0.0, nil}, fn [lower, upper], {max_str, _base} = acc ->
        dt = upper["tmpc"] - lower["tmpc"]

        if dt > 0 and dt > max_str do
          {dt, lower["hght"] - sfc_hght}
        else
          acc
        end
      end)

    case result do
      {str, base} when str > 0 -> base
      _ -> nil
    end
  end

  defp duct_field(nil, _key), do: nil
  defp duct_field([], _key), do: nil

  # Ducts arrive in two shapes:
  #
  # * `SoundingParams.detect_ducts` emits string keys `"base"` / `"strength"`.
  # * HrrrNativeClient and the Rust f00 pipeline emit atom keys
  #   `:base_m` / `:top_m` / `:thickness_m` / `:min_freq_ghz`, which get
  #   atomized by `ProfilesFile.read/1` via the `@mp_atom_keys` whitelist.
  #
  # Read both shapes so `/weather` renders Duct Base / Duct Strength
  # layers regardless of producer. `:thickness_m` stands in for
  # "strength" when the native shape is in use — thicker ducts trap a
  # wider frequency range, so it's the most meaningful scalar.
  defp duct_field(ducts, "base") do
    ducts
    |> Enum.map(&duct_value(&1, ["base", "base_m"]))
    |> Enum.reject(&is_nil/1)
    |> Enum.min(fn -> nil end)
  end

  defp duct_field(ducts, "strength") do
    ducts
    |> Enum.map(&duct_value(&1, ["strength", "thickness_m"]))
    |> Enum.reject(&is_nil/1)
    |> Enum.max(fn -> nil end)
  end

  # First present key wins. Keys are normalized to strings first because
  # a duct read back through `ProfilesFile.read/1` carries atom keys
  # (`:base_m`) while a `Msgpax`-decoded or sounding-derived one carries
  # strings — indexing an atom-keyed map with a string silently returns
  # nil, which blanks the layer.
  defp duct_value(duct, keys) when is_map(duct) do
    d = stringify_keys(duct)
    Enum.find_value(keys, fn key -> d[key] end)
  end

  defp duct_value(_duct, _keys), do: nil

  defp stringify_keys(map), do: Map.new(map, fn {k, v} -> {to_string(k), v} end)
end