prop/lib/microwaveprop/beacons/range_estimate.ex

defmodule Microwaveprop.Beacons.RangeEstimate do
  @moduledoc """
  Estimates a beacon's reception range at several signal-strength tiers.

  For each tier we solve a link budget of

      Rx_dBm = EIRP_dBm + Rx_gain_dBi - FSPL(d, f) - atm_loss_per_km * d

  for the distance `d` at which the received power equals the tier threshold.
  `EIRP_dBm` comes directly from the beacon's stored `power_mw` — the field
  already represents EIRP (TX power × antenna gain). Free-space path loss uses
  the standard formula and atmospheric absorption comes from the band's
  O₂/H₂O coefficients in `BandConfig`.

  The result is then multiplied by a propagation-score factor derived from the
  latest `propagation_scores` row at the beacon's lat/lon — score 50 → 1.0x,
  score 0 → 0.5x, score 100 → 1.5x — so current HRRR conditions shift the
  rings in or out.
  """

  alias Microwaveprop.Propagation
  alias Microwaveprop.Propagation.BandConfig

  # Signal-strength tiers and their RX sensitivity thresholds (dBm).
  # Colors match the map_live / propagation_map_hook palette.
  @tiers [
    %{label: "Excellent", rx_dbm: -100, color: "#00ffa3"},
    %{label: "Good", rx_dbm: -115, color: "#7dffd4"},
    %{label: "Marginal", rx_dbm: -125, color: "#ffe566"},
    %{label: "Weak CW", rx_dbm: -135, color: "#ff9044"},
    %{label: "Detection", rx_dbm: -145, color: "#ff4f4f"}
  ]

  # Assume the receiving station is an average amateur microwave station
  # (dish/horn + low-noise preamp).
  @rx_gain_dbi 20.0

  @doc """
  Convert a power in milliwatts to dBm. Returns `-999.9` for non-positive input.
  """
  @spec mw_to_dbm(number()) :: float()
  def mw_to_dbm(mw) when is_number(mw) and mw > 0, do: 10.0 * :math.log10(mw)
  def mw_to_dbm(_), do: -999.9

  @doc """
  Returns the closest configured band frequency (in MHz) to the given beacon
  frequency. e.g. `nearest_band_mhz(10368.1) == 10_000`.
  """
  @spec nearest_band_mhz(number()) :: integer()
  def nearest_band_mhz(freq_mhz) when is_number(freq_mhz) do
    BandConfig.all_freqs()
    |> Enum.min_by(fn b -> abs(b - freq_mhz) end)
  end

  @doc """
  Estimate a beacon's reception range.

  Returns a map with band info, current score, and a list of rings sorted
  weakest-distance first (strongest RX tier → shortest radius).
  """
  @spec estimate(Microwaveprop.Beacons.Beacon.t()) :: map()
  def estimate(beacon) do
    band_mhz = nearest_band_mhz(beacon.frequency_mhz)
    band_config = BandConfig.get(band_mhz)

    detail = Propagation.point_detail(band_mhz, beacon.lat, beacon.lon)
    score = (detail && detail.score) || 50
    valid_time = detail && detail.valid_time

    f_mhz = beacon.frequency_mhz * 1.0
    eirp_dbm = mw_to_dbm(beacon.power_mw || 0.0)

    atm_per_km = atm_loss_per_km(band_config)
    score_mult = 0.5 + score / 100.0

    rings =
      @tiers
      |> Enum.map(fn tier ->
        d_phys = solve_range(eirp_dbm, @rx_gain_dbi, tier.rx_dbm, f_mhz, atm_per_km)
        radius = Float.round(d_phys * score_mult, 1)

        %{
          label: tier.label,
          rx_dbm: tier.rx_dbm,
          color: tier.color,
          radius_km: radius
        }
      end)
      |> Enum.filter(fn ring -> ring.radius_km > 0.5 end)
      |> Enum.sort_by(& &1.radius_km)

    %{
      beacon_id: beacon.id,
      band_mhz: band_mhz,
      band_label: band_config && band_config.label,
      score: score,
      score_mult: Float.round(score_mult, 2),
      valid_time: valid_time,
      eirp_dbm: Float.round(eirp_dbm, 1),
      atm_per_km: Float.round(atm_per_km, 3),
      rings: rings
    }
  end

  # Total dB/km atmospheric attenuation from O2 + water vapor. Uses a moderate
  # absolute humidity of 10 g/m³ as a default — the HRRR-derived score already
  # captures humidity variability, so using a fixed value here keeps the
  # physics clean.
  defp atm_loss_per_km(nil), do: 0.0

  defp atm_loss_per_km(band_config) do
    o2 = Map.get(band_config, :o2_db_km, 0.0)
    h2o_coeff = Map.get(band_config, :h2o_coeff, 0.0)
    o2 + h2o_coeff * 10.0
  end

  # Solve `FSPL(d) + atm_per_km * d = budget_db` for d via bisection.
  # `budget_db = EIRP + Rx_gain - threshold`.
  defp solve_range(eirp_dbm, rx_gain, threshold_dbm, f_mhz, atm_per_km) do
    budget = eirp_dbm + rx_gain - threshold_dbm
    log_f_const = 20.0 * :math.log10(f_mhz) + 32.44

    cond do
      budget <= log_f_const ->
        # Even at 1 km the budget is already negative → ring is effectively 0.
        0.0

      true ->
        bisect(0.01, 5000.0, budget, log_f_const, atm_per_km, 50)
    end
  end

  defp bisect(lo, hi, _budget, _log_f, _atm, 0), do: (lo + hi) / 2.0

  defp bisect(lo, hi, budget, log_f, atm, iters) do
    mid = (lo + hi) / 2.0
    val = 20.0 * :math.log10(mid) + log_f + atm * mid

    if val > budget do
      bisect(lo, mid, budget, log_f, atm, iters - 1)
    else
      bisect(mid, hi, budget, log_f, atm, iters - 1)
    end
  end
end