prop/lib/microwaveprop/propagation/path_analysis.ex

defmodule Microwaveprop.Propagation.PathAnalysis do
  @moduledoc """
  Path-level analysis for contact display: duct detection, propagation
  mechanism explanation, data source summaries, and elevation profiles.

  Wraps Propagation.Duct and Propagation.MechanismClassifier for the
  algorithm-heavy work while providing display-ready results for the
  ContactLive.Show LiveView.

  Extracted from ContactLive.Show (2026-07) to eliminate ~750 lines of
  duplicated duct-detection and mechanism-classification logic.
  """

  alias Microwaveprop.Repo
  alias Microwaveprop.Terrain.ElevationClient
  alias Microwaveprop.Terrain.TerrainAnalysis

  require Logger

  # ── Gradient thresholds (centralized — shared with Propagation.Duct) ──
  # ITU-R P.453 / P.834 definitions:
  #   dN/dh < -157 N/km → ducting (M-profile decreases with height)
  #   dN/dh < -100 N/km → super-refractive
  #   dN/dh < -79  N/km → enhanced refraction
  #   dN/dh ≥ -79  N/km → standard atmosphere
  @ducting_threshold -157
  @super_refractive_threshold -100
  @enhanced_refraction_threshold -79

  @earth_radius_m 6_371_000.0

  # ── Public API ───────────────────────────────────────────────────────

  @doc """
  Compute the elevation profile and duct layers for a contact path.

  Returns a map with `:points`, `:freq_mhz`, `:dist_km`, `:k_factor`,
  `:fwd_az`, `:rev_az`, `:fwd_el`, `:rev_el`, `:verdict`,
  `:first_obstruction_km`, and `:ducts`, or nil when coordinates are
  unavailable or elevation data cannot be fetched.
  """
  @spec compute_elevation_profile(map(), list(), list()) :: map() | nil
  def compute_elevation_profile(contact, hrrr_path, soundings) do
    with %{"lat" => lat1} <- contact.pos1,
         lon1 when is_number(lon1) <- contact.pos1["lon"],
         %{"lat" => lat2} <- contact.pos2,
         lon2 when is_number(lon2) <- contact.pos2["lon"],
         {:ok, profile} <- safe_fetch_elevation(lat1, lon1, lat2, lon2) do
      freq_ghz = band_to_ghz(contact.band)
      dist_km = Microwaveprop.Geo.haversine_km(lat1, lon1, lat2, lon2)
      # Default to 10 ft (3.048 m) AGL when the operator didn't record an
      # antenna height — better than pretending the antenna is on the dirt.
      default_ht_m = 3.048
      ant_ht_a_m = ft_to_m(contact.height1_ft) || default_ht_m
      ant_ht_b_m = ft_to_m(contact.height2_ft) || default_ht_m
      analysis = TerrainAnalysis.analyse(profile, dist_km, freq_ghz, ant_ht_a: ant_ht_a_m, ant_ht_b: ant_ht_b_m)

      fwd_az = initial_bearing(lat1, lon1, lat2, lon2)
      rev_az = initial_bearing(lat2, lon2, lat1, lon1)

      tx_elev = hd(profile).elev + ant_ht_a_m
      rx_elev = List.last(profile).elev + ant_ht_b_m
      fwd_el = elevation_angle(tx_elev, rx_elev, dist_km * 1000)
      rev_el = elevation_angle(rx_elev, tx_elev, dist_km * 1000)

      first_obs =
        case Enum.find(analysis.points, & &1.obstructed) do
          nil -> nil
          p -> p.dist_km
        end

      ducts =
        hrrr_path
        |> extract_ducts(soundings, tx_elev)
        |> merge_nearby_ducts()
        |> Enum.sort_by(& &1.strength, :desc)
        |> Enum.take(3)
        |> mark_likely_duct(tx_elev, rx_elev)

      %{
        points:
          Enum.map(analysis.points, fn p ->
            %{elev: p.elev, beam: p.beam, r1: p.r1, dist_km: p.dist_km}
          end),
        freq_mhz: freq_ghz * 1000,
        dist_km: dist_km,
        k_factor: analysis.k_factor,
        fwd_az: fwd_az,
        rev_az: rev_az,
        fwd_el: fwd_el,
        rev_el: rev_el,
        verdict: analysis.verdict,
        first_obstruction_km: first_obs,
        ducts: ducts
      }
    else
      _ -> nil
    end
  end

  @doc """
  Build a summary map of the data sources available for this path.
  """
  @spec build_data_sources(map() | nil, list(), map() | nil, map(), map() | nil) :: map()
  def build_data_sources(hrrr, hrrr_path, terrain, weather, elevation_profile) do
    %{
      hrrr: build_hrrr_source(hrrr, hrrr_path),
      elevation: build_elevation_source(elevation_profile),
      terrain: build_terrain_source(terrain),
      obs_count: length(weather.surface_observations),
      sounding_count: length(weather.soundings)
    }
  end

  @doc """
  Build the propagation summary text for the contact.
  """
  @spec build_summary(map() | nil, map() | nil, map() | nil, list(), float() | nil, integer(), map()) ::
          String.t()
  def build_summary(terrain, elevation_profile, hrrr, soundings, dist_km, band_mhz, contact) do
    terrain_status = terrain_summary(terrain, elevation_profile, contact)
    mechanism = propagation_mechanism(terrain, elevation_profile, hrrr, soundings, dist_km)
    band_note = band_summary(band_mhz, hrrr)

    [terrain_status, mechanism, band_note]
    |> Enum.reject(&is_nil/1)
    |> Enum.join(" ")
  end

  @doc """
  Build the detail notes list for the propagation analysis panel.
  """
  @spec build_details(map() | nil, list(), list(), map() | nil, integer(), map()) :: [String.t()]
  def build_details(hrrr, hrrr_path, soundings, elevation_profile, _band_mhz, contact) do
    Enum.reject(
      [
        antenna_heights_detail(contact),
        refractivity_detail(hrrr),
        boundary_layer_detail(hrrr),
        path_duct_count_detail(hrrr_path),
        sounding_ducting_detail(soundings),
        duct_layer_detail(elevation_profile)
      ],
      &is_nil/1
    )
  end

  @doc """
  Human-readable description of a detected duct layer.
  """
  @spec duct_description(map()) :: String.t()
  def duct_description(%{source: source, base_m_msl: base, top_m_msl: top, strength: strength}) do
    base_ft = round(base * 3.28084)
    top_ft = round(top * 3.28084)

    src =
      case source do
        "sounding" -> "sounding-detected"
        "inferred" -> "estimated"
        _ -> "detected"
      end

    "#{src} layer at #{base_ft}-#{top_ft} ft, #{strength} M-units"
  end

  def duct_description(_), do: "detected layer"

  @doc """
  Determine the most likely propagation mechanism text for a contact,
  considering terrain, ducting conditions, and path length.
  """
  @spec propagation_mechanism(map() | nil, map() | nil, map() | nil, list(), float() | nil) ::
          String.t() | nil
  def propagation_mechanism(terrain, elevation_profile, hrrr, soundings, dist_km) do
    blocked? = terrain && terrain.verdict == "BLOCKED"
    ducting? = has_ducting?(hrrr, soundings)
    enhanced_refraction? = enhanced_refraction?(hrrr)
    long_path? = dist_km && dist_km > 100

    ducts = (elevation_profile && elevation_profile.ducts) || []
    likely_duct = Enum.find(ducts, & &1[:likely])

    props = %{
      blocked?: blocked?,
      ducting?: ducting?,
      enhanced_refraction?: enhanced_refraction?,
      long_path?: long_path?,
      likely_duct: likely_duct,
      terrain: terrain,
      hrrr: hrrr
    }

    if blocked? do
      blocked_mechanism(props)
    else
      clear_path_mechanism(props)
    end
  end

  @doc """
  Returns true if HRRR or sounding data indicates ducting conditions.
  """
  @spec has_ducting?(map() | nil, list()) :: boolean()
  def has_ducting?(hrrr, soundings) do
    hrrr_ducting = hrrr && hrrr.ducting_detected
    sounding_ducting = is_list(soundings) && Enum.any?(soundings, & &1.ducting_detected)
    hrrr_ducting || sounding_ducting
  end

  @doc """
  Returns true if the HRRR refractivity gradient indicates enhanced
  (super-refractive) conditions.
  """
  @spec enhanced_refraction?(map() | nil) :: boolean()
  def enhanced_refraction?(nil), do: false

  def enhanced_refraction?(hrrr) do
    is_number(hrrr.min_refractivity_gradient) && hrrr.min_refractivity_gradient < @super_refractive_threshold
  end

  @doc """
  Translates a refractivity gradient (N-units/km) into a human label.

  Uses ITU-R P.834 conventions (same thresholds as Propagation.Duct).
  """
  @spec refractivity_label(integer()) :: String.t()
  def refractivity_label(g) when g < @ducting_threshold, do: "ducting"
  def refractivity_label(g) when g < @super_refractive_threshold, do: "super-refractive"
  def refractivity_label(g) when g < @enhanced_refraction_threshold, do: "enhanced"
  def refractivity_label(_), do: "normal"

  # ── Duct extraction pipeline ────────────────────────────────────────

  # Extract duct layers from best available source across all HRRR path profiles:
  # 1. HRRR explicit ducts (M-profile detected) from any path point
  # 2. Sounding explicit ducts (high vertical resolution, most reliable)
  # 3. Inferred from strongest HRRR refractivity gradient along path
  defp extract_ducts(hrrr_path, soundings, surface_elev) do
    hrrr_ducts = extract_hrrr_ducts(hrrr_path, surface_elev)

    cond do
      hrrr_ducts != [] ->
        hrrr_ducts

      (sounding_ducts = extract_sounding_ducts(soundings, surface_elev)) != [] ->
        sounding_ducts

      true ->
        # Use the profile with the strongest (most negative) gradient
        best =
          hrrr_path
          |> Enum.filter(&is_number(&1.min_refractivity_gradient))
          |> Enum.min_by(& &1.min_refractivity_gradient, fn -> nil end)

        infer_duct_from_gradient(best, surface_elev)
    end
  end

  defp extract_hrrr_ducts(hrrr_path, surface_elev) when is_list(hrrr_path) do
    hrrr_path
    |> Enum.flat_map(&duct_characteristics_to_maps(&1.duct_characteristics, surface_elev, "hrrr"))
    |> Enum.uniq_by(fn d -> {round(d.base_m_msl), round(d.top_m_msl)} end)
  end

  defp duct_characteristics_to_maps(ducts, surface_elev, source) when is_list(ducts) and ducts != [] do
    Enum.map(ducts, fn d ->
      %{
        base_m_msl: surface_elev + (d["base"] || 0),
        top_m_msl: surface_elev + (d["top"] || 0),
        strength: d["strength"],
        source: source
      }
    end)
  end

  defp duct_characteristics_to_maps(_ducts, _surface_elev, _source), do: []

  defp extract_sounding_ducts(soundings, surface_elev) when is_list(soundings) do
    soundings
    |> Enum.filter(& &1.ducting_detected)
    |> Enum.flat_map(&duct_characteristics_to_maps(&1.duct_characteristics, surface_elev, "sounding"))
    |> Enum.uniq_by(fn d -> {round(d.base_m_msl), round(d.top_m_msl)} end)
  end

  defp extract_sounding_ducts(_, _), do: []

  # HRRR pressure levels are too coarse (~250m) to detect thin ducting layers
  # via the M-profile method. When min_refractivity_gradient indicates enhanced
  # propagation, infer a duct layer within the boundary layer.
  defp infer_duct_from_gradient(nil, _surface_elev), do: []

  defp infer_duct_from_gradient(hrrr, surface_elev) do
    grad = hrrr.min_refractivity_gradient
    hpbl = hrrr.hpbl_m

    if is_number(grad) and grad < @super_refractive_threshold and is_number(hpbl) and hpbl > 0 do
      duct_top = min(hpbl * 0.6, 500)
      strength = abs(grad) / 10

      [
        %{
          base_m_msl: surface_elev,
          top_m_msl: surface_elev + duct_top,
          strength: Float.round(strength, 1),
          source: "inferred"
        }
      ]
    else
      []
    end
  end

  # Merge duct layers that overlap or are within 100m of each other
  defp merge_nearby_ducts(ducts) do
    ducts
    |> Enum.sort_by(& &1.base_m_msl)
    |> Enum.reduce([], fn duct, acc ->
      case acc do
        [prev | rest] when duct.base_m_msl <= prev.top_m_msl + 100 ->
          merged = %{
            prev
            | top_m_msl: max(prev.top_m_msl, duct.top_m_msl),
              strength: max(prev.strength, duct.strength)
          }

          [merged | rest]

        _ ->
          [duct | acc]
      end
    end)
    |> Enum.reverse()
  end

  # Mark the duct most likely carrying the signal. The signal enters at
  # surface level from each end, so the duct whose base is closest to
  # the average endpoint elevation is the most probable propagation path.
  defp mark_likely_duct([], _tx_elev, _rx_elev), do: []

  defp mark_likely_duct(ducts, tx_elev, rx_elev) do
    avg_elev = (tx_elev + rx_elev) / 2

    {_, likely_idx} =
      ducts
      |> Enum.with_index()
      |> Enum.min_by(fn {d, _idx} ->
        # Distance from average endpoint elevation to the duct layer
        # Prefer ducts that the endpoints are actually inside of
        if avg_elev >= d.base_m_msl and avg_elev <= d.top_m_msl do
          -d.strength
        else
          min(abs(d.base_m_msl - avg_elev), abs(d.top_m_msl - avg_elev))
        end
      end)

    Enum.with_index(ducts, fn d, i ->
      Map.put(d, :likely, i == likely_idx)
    end)
  end

  # ── Elevation fetch helper ──────────────────────────────────────────

  defp safe_fetch_elevation(lat1, lon1, lat2, lon2) do
    ElevationClient.fetch_elevation_profile(lat1, lon1, lat2, lon2, 256)
  rescue
    e ->
      Logger.warning("Elevation profile failed: #{Exception.message(e)}")
      {:error, :elevation_unavailable}
  end

  # ── Conversion helpers ──────────────────────────────────────────────

  defp band_to_ghz(nil), do: 10.0

  defp band_to_ghz(band) do
    band
    |> Decimal.to_float()
    |> Kernel./(1000)
  end

  defp ft_to_m(nil), do: nil
  defp ft_to_m(ft) when is_integer(ft), do: ft * 0.3048

  # ── Geodesy helpers ─────────────────────────────────────────────────

  defp initial_bearing(lat1, lon1, lat2, lon2) do
    rlat1 = deg_to_rad(lat1)
    rlat2 = deg_to_rad(lat2)
    dlon = deg_to_rad(lon2 - lon1)

    y = :math.sin(dlon) * :math.cos(rlat2)

    x =
      :math.cos(rlat1) * :math.cos(rlat2) -
        :math.cos(rlat1) * :math.sin(rlat2) * :math.cos(dlon)

    bearing = y |> :math.atan2(x) |> rad_to_deg()
    Float.round(rem_float(bearing + 360, 360), 1)
  end

  defp elevation_angle(h_tx, h_rx, dist_m) do
    k_r = 4 / 3 * @earth_radius_m
    angle = :math.atan2(h_rx - h_tx, dist_m) - dist_m / (2 * k_r)
    Float.round(rad_to_deg(angle), 2)
  end

  defp deg_to_rad(deg), do: deg * :math.pi() / 180
  defp rad_to_deg(rad), do: rad * 180 / :math.pi()
  defp rem_float(a, b), do: a - Float.floor(a / b) * b

  # ── Data sources building ───────────────────────────────────────────

  defp build_hrrr_source(nil, _), do: nil

  defp build_hrrr_source(hrrr, hrrr_path) do
    levels = if hrrr.profile, do: length(hrrr.profile), else: 0

    %{
      profile_count: length(hrrr_path),
      levels: levels,
      valid_time: Calendar.strftime(hrrr.valid_time, "%Y-%m-%d %H:%M UTC"),
      lat: :erlang.float_to_binary(hrrr.lat / 1, decimals: 2),
      lon: :erlang.float_to_binary(hrrr.lon / 1, decimals: 2)
    }
  end

  defp build_elevation_source(nil), do: nil

  defp build_elevation_source(ep) do
    %{
      samples: length(ep.points),
      dist: Microwaveprop.Format.distance_km(ep.dist_km),
      source: "SRTM 1-arcsecond",
      duct_count: length(ep.ducts)
    }
  end

  defp build_terrain_source(nil), do: nil

  defp build_terrain_source(terrain) do
    %{
      samples: terrain.sample_count,
      verdict: terrain.verdict,
      max_elev: :erlang.float_to_binary((terrain.max_elevation_m || 0) / 1, decimals: 0),
      diffraction: :erlang.float_to_binary((terrain.diffraction_db || 0) / 1, decimals: 1)
    }
  end

  # ── Terrain summary text ────────────────────────────────────────────

  defp terrain_summary(nil, nil, _contact), do: "No terrain data available."

  defp terrain_summary(terrain, elevation_profile, contact) do
    verdict = (elevation_profile && elevation_profile.verdict) || (terrain && terrain.verdict)
    describe_verdict(verdict, elevation_profile, contact)
  end

  defp describe_verdict("CLEAR", _, _), do: "Line of sight is clear between stations."

  defp describe_verdict("FRESNEL_MINOR", _, _), do: "Line of sight is clear but with minor Fresnel zone encroachment."

  defp describe_verdict("FRESNEL_PARTIAL", _, _), do: "Line of sight has partial Fresnel zone obstruction."

  defp describe_verdict("BLOCKED", elevation_profile, contact) do
    obs_dist = elevation_profile && elevation_profile.first_obstruction_km
    blocked_message(obs_dist, contact)
  end

  defp describe_verdict(_, _, _), do: nil

  defp blocked_message(nil, _contact), do: "Path is terrain-obstructed."

  defp blocked_message(obs_dist, contact) do
    origin = contact && contact.station1
    origin_label = if origin, do: origin, else: "station 1"

    "Path is terrain-obstructed at #{:erlang.float_to_binary(obs_dist * 0.621371, decimals: 1)} mi (#{:erlang.float_to_binary(obs_dist, decimals: 1)} km) from #{origin_label}."
  end

  # ── Mechanism text generation ───────────────────────────────────────

  defp blocked_mechanism(%{ducting?: true, likely_duct: duct}) when not is_nil(duct) do
    "Signal likely propagated via atmospheric ducting (#{duct_description(duct)}), bending over the terrain obstruction."
  end

  defp blocked_mechanism(%{ducting?: true}) do
    "Ducting conditions detected — signal likely propagated through a tropospheric duct, bypassing the terrain obstruction."
  end

  defp blocked_mechanism(%{enhanced_refraction?: true, hrrr: hrrr}) do
    grad = round(hrrr.min_refractivity_gradient)

    "Enhanced refraction (dN/dh = #{grad} N-units/km) likely bent the signal over the obstruction. Conditions are super-refractive but below full ducting threshold."
  end

  defp blocked_mechanism(%{terrain: terrain}) do
    diff_db = terrain && terrain.diffraction_db
    blocked_diffraction_message(diff_db)
  end

  defp blocked_diffraction_message(diff_db) when is_number(diff_db) and diff_db > 0 do
    "Path is obstructed with #{:erlang.float_to_binary(diff_db, decimals: 1)} dB diffraction loss. Contact may have been enabled by knife-edge diffraction or tropospheric scatter."
  end

  defp blocked_diffraction_message(_) do
    "Path is obstructed. Contact likely enabled by tropospheric scatter or transient ducting conditions."
  end

  defp clear_path_mechanism(%{ducting?: true, long_path?: true, likely_duct: duct}) when not is_nil(duct) do
    "Ducting conditions present (#{duct_description(duct)}) — likely extending range beyond normal line of sight."
  end

  defp clear_path_mechanism(%{ducting?: true, long_path?: true}) do
    "Ducting conditions present — likely contributing to extended range."
  end

  defp clear_path_mechanism(%{enhanced_refraction?: true, long_path?: true}) do
    "Enhanced refraction present — may have contributed to extended range."
  end

  defp clear_path_mechanism(_), do: nil

  # ── Detail text helpers ─────────────────────────────────────────────

  defp antenna_heights_detail(%{height1_ft: h1, height2_ft: h2} = contact) when is_integer(h1) or is_integer(h2) do
    s1 = contact.station1 || "Station 1"
    s2 = contact.station2 || "Station 2"

    parts =
      Enum.reject(
        [if(is_integer(h1), do: "#{s1} @ #{h1} ft AGL"), if(is_integer(h2), do: "#{s2} @ #{h2} ft AGL")],
        &is_nil/1
      )

    "Antenna heights: #{Enum.join(parts, ", ")}."
  end

  defp antenna_heights_detail(_), do: nil

  defp path_duct_count_detail(hrrr_path) when is_list(hrrr_path) and hrrr_path != [] do
    total = length(hrrr_path)
    ducting = Enum.count(hrrr_path, &(&1 && &1.ducting_detected))

    if ducting > 0 do
      "Tropo duct detected at #{ducting}/#{total} HRRR samples along the path."
    end
  end

  defp path_duct_count_detail(_), do: nil

  defp refractivity_detail(%{min_refractivity_gradient: grad}) when is_number(grad) do
    g = round(grad)
    label = refractivity_label(g)
    "Refractivity gradient: #{g} N-units/km (#{label})."
  end

  defp refractivity_detail(_), do: nil

  defp boundary_layer_detail(%{hpbl_m: hpbl}) when is_number(hpbl) do
    h = round(hpbl)
    note = if h < 300, do: " — shallow boundary layer favors ducting", else: ""
    "Boundary layer depth: #{h} m#{note}."
  end

  defp boundary_layer_detail(_), do: nil

  defp sounding_ducting_detail(soundings) when is_list(soundings) do
    soundings = Repo.preload(soundings, :station)
    ducting = Enum.filter(soundings, & &1.ducting_detected)

    if ducting == [] do
      nil
    else
      names = Enum.map_join(ducting, ", ", fn s -> s.station.name || s.station.station_code end)
      "Ducting detected by soundings at: #{names}."
    end
  end

  defp sounding_ducting_detail(_), do: nil

  defp duct_layer_detail(%{ducts: ducts}) when is_list(ducts) do
    likely = Enum.find(ducts, & &1[:likely])
    if likely, do: "Most likely propagation duct: #{duct_description(likely)}."
  end

  defp duct_layer_detail(_), do: nil

  # ── Band-specific summary text ──────────────────────────────────────

  defp band_summary(band_mhz, hrrr) when band_mhz <= 12_000 do
    if hrrr && is_number(hrrr.pwat_mm) && hrrr.pwat_mm > 20 do
      "At 10 GHz, the elevated moisture (PWAT #{:erlang.float_to_binary(hrrr.pwat_mm, decimals: 1)} mm) enhances refractivity and ducting potential."
    end
  end

  defp band_summary(band_mhz, hrrr) when band_mhz >= 24_000 do
    if hrrr && is_number(hrrr.pwat_mm) && hrrr.pwat_mm > 25 do
      "At #{round(band_mhz / 1000)} GHz, high moisture (PWAT #{:erlang.float_to_binary(hrrr.pwat_mm, decimals: 1)} mm) causes significant water vapor absorption."
    end
  end

  defp band_summary(_, _), do: nil
end