towerops/lib/towerops/workers/coverage_worker.ex

defmodule Towerops.Workers.CoverageWorker do
  @moduledoc """
  Computes an RF coverage prediction for a single coverage record.

  Pipeline:

    1. Resolve antenna and centre coordinates.
    2. Compute a WGS84 bbox of the coverage area from `radius_m`.
    3. Pull a terrain elevation grid for that bbox via `Towerops.Lidar`.
    4. For each output pixel, sample a path profile, run propagation,
       apply the antenna pattern, and produce an RSSI in dBm.
    5. Write GeoTIFF + colored PNG to `priv/static/coverage/<org>/<id>/`.

  Compute is parallelised across BEAM schedulers via `Task.async_stream/3`
  on row bands. Progress is broadcast on
  `Towerops.Coverages.topic/1` and `org_topic/1` so LiveViews can render
  live status.

  Org scoping is enforced: job args carry `organization_id` and the
  worker refuses to run if the loaded coverage's org doesn't match.
  """

  use Oban.Worker,
    queue: :coverage,
    max_attempts: 3,
    unique: [
      fields: [:args],
      keys: [:coverage_id],
      states: [:available, :scheduled, :executing, :retryable]
    ]

  alias Towerops.Coverages
  alias Towerops.Coverages.Antenna
  alias Towerops.Coverages.Buildings
  alias Towerops.Coverages.Coverage
  alias Towerops.Coverages.Notifier
  alias Towerops.Coverages.Profile
  alias Towerops.Coverages.Propagation
  alias Towerops.Coverages.Raster
  alias Towerops.Coverages.TreeCanopy
  alias Towerops.Repo

  require Logger

  # Number of elevation samples per profile. Higher = more accurate
  # diffraction detection, lower = faster compute. 64 is a good
  # trade-off for WISP-scale paths.
  @profile_samples 64

  # SM (Subscriber Module) install heights to compute. Mirrors the
  # cnHeat vertical slider granularity — a WISP-realistic ladder
  # spanning ground-mount through tower-mount.
  @sm_height_tiers_m [1.83, 3.05, 4.57, 6.10, 9.14, 12.19]
  @default_tier_m 3.05

  @impl Oban.Worker
  def perform(%Oban.Job{args: %{"coverage_id" => coverage_id, "organization_id" => organization_id}}) do
    case Repo.get(Coverage, coverage_id) do
      nil ->
        Logger.warning("CoverageWorker: coverage #{coverage_id} not found, skipping")
        :ok

      %Coverage{organization_id: ^organization_id} = coverage ->
        coverage = Repo.preload(coverage, :site)
        run(coverage)

      %Coverage{} ->
        Logger.error(
          "CoverageWorker: organization mismatch for #{coverage_id} " <>
            "(expected #{organization_id}); refusing to run"
        )

        {:cancel, :organization_mismatch}
    end
  end

  defp run(coverage) do
    started = System.monotonic_time(:millisecond)
    log(coverage, "starting compute (cell=#{coverage.cell_size_m} m, radius=#{coverage.radius_m} m)")
    _ = update_progress(coverage, "computing", 1)

    try do
      with {:ok, antenna} <- resolve_antenna(coverage),
           _ = log(coverage, "antenna resolved: #{antenna.slug}"),
           {:ok, {lat, lon}} <- resolve_location(coverage),
           {:ok, bbox} <- compute_bbox(coverage, lat, lon),
           _ = update_progress(coverage, "computing", 5),
           _ = log(coverage, "fetching terrain grid for bbox #{inspect(bbox)}"),
           {:ok, grid} <- fetch_terrain(bbox, coverage.cell_size_m, lat),
           _ =
             log(
               coverage,
               "terrain grid: #{grid.ncols}x#{grid.nrows} cells (#{grid.ncols * grid.nrows} samples)"
             ),
           clutter = Buildings.for_bbox(bbox),
           _ = log(coverage, "loaded #{length(clutter)} building footprints for clutter"),
           canopy = fetch_canopy(bbox, coverage.cell_size_m, lat),
           _ = update_progress(coverage, "computing", 15),
           _ = log(coverage, "starting per-tier compute (#{length(@sm_height_tiers_m)} tiers x LOS+NLOS)"),
           {:ok, tiers} <-
             compute_all_tiers(coverage, antenna, lat, lon, bbox, grid, clutter, canopy) do
        elapsed = System.monotonic_time(:millisecond) - started
        log(coverage, "compute finished in #{elapsed} ms; finalizing")
        finalize(coverage, tiers)
      else
        {:error, reason} -> fail(coverage, reason)
      end
    rescue
      e ->
        Logger.error(
          "CoverageWorker[#{coverage.id}]: unexpected exception: #{Exception.format(:error, e, __STACKTRACE__)}"
        )

        fail(coverage, {:terrain_unavailable, {:exception, Exception.message(e)}})
    end
  end

  defp log(coverage, message) do
    Logger.info("CoverageWorker[#{coverage.id}]: #{message}")
    :ok
  end

  # Best-effort canopy fetch — non-fatal if disabled / out-of-area.
  # Only NLOS sampling consumes it, so a nil here just means
  # "compute as if no trees are present" for that mode.
  defp fetch_canopy(bbox, cell_size_m, lat) do
    cell_size_deg = cell_size_m / (111_000.0 * :math.cos(lat * :math.pi() / 180.0))

    case TreeCanopy.fetch_grid(bbox, cell_size_deg) do
      {:ok, grid} ->
        grid

      {:error, reason} ->
        Logger.info("CoverageWorker: tree canopy unavailable (#{inspect(reason)}); skipping")
        nil
    end
  end

  # Runs compute_pixels twice per SM-height tier (LOS / NLOS) and
  # writes per-tier raster + PNG sets. LOS treats clutter as ground
  # baseline (cnHeat "Height above clutter"); NLOS samples the
  # prefetched building list into the path profile so rooftop
  # diffraction shadows show up ("Height above ground").
  defp compute_all_tiers(coverage, antenna, lat, lon, bbox, grid, clutter, canopy) do
    indexed = Enum.with_index(@sm_height_tiers_m)

    initial = {%{los: %{}, nlos: %{}}, nil}
    env = %{antenna: antenna, lat: lat, lon: lon, bbox: bbox, grid: grid, clutter: clutter, canopy: canopy}

    {results, error} =
      Enum.reduce_while(indexed, initial, fn {height_m, idx}, {acc, _last} ->
        case compute_tier(coverage, env, height_m) do
          {:ok, %{los: los_paths, nlos: nlos_paths}} ->
            _ =
              update_progress(
                coverage,
                "computing",
                15 + round((idx + 1) / length(@sm_height_tiers_m) * 75)
              )

            new_acc = %{
              los: Map.put(acc.los, height_m, los_paths),
              nlos: Map.put(acc.nlos, height_m, nlos_paths)
            }

            {:cont, {new_acc, nil}}

          {:error, reason} ->
            {:halt, {acc, reason}}
        end
      end)

    case {error, Map.get(results.los, @default_tier_m)} do
      {nil, %{} = default} ->
        {:ok, %{tiers: results.los, nlos_tiers: results.nlos, default: default}}

      {nil, nil} ->
        {:error, :default_tier_missing}

      {reason, _} ->
        {:error, reason}
    end
  end

  defp compute_tier(coverage, env, height_m) do
    cov = %{coverage | receiver_height_m: height_m}

    los_ctx = Map.merge(env, %{clutter: [], canopy: nil})
    nlos_ctx = Map.merge(env, %{clutter: env.clutter, canopy: env.canopy})

    log(coverage, "tier h=#{Float.round(height_m, 2)} m: computing LOS pixels")
    t_los = System.monotonic_time(:millisecond)

    with {:ok, los_out} <- compute_pixels(cov, los_ctx, env.lat, env.lon),
         _ =
           log(
             coverage,
             "tier h=#{Float.round(height_m, 2)} m: LOS pixels done in #{System.monotonic_time(:millisecond) - t_los} ms; writing raster"
           ),
         {:ok, los_paths} <-
           Raster.write(coverage, los_out.pixels, los_out.dims, tier_suffix(height_m, :los)),
         _ = log(coverage, "tier h=#{Float.round(height_m, 2)} m: computing NLOS pixels"),
         t_nlos = System.monotonic_time(:millisecond),
         {:ok, nlos_out} <- compute_pixels(cov, nlos_ctx, env.lat, env.lon),
         _ =
           log(
             coverage,
             "tier h=#{Float.round(height_m, 2)} m: NLOS pixels done in #{System.monotonic_time(:millisecond) - t_nlos} ms; writing raster"
           ),
         {:ok, nlos_paths} <-
           Raster.write(coverage, nlos_out.pixels, nlos_out.dims, tier_suffix(height_m, :nlos)) do
      {:ok, %{los: los_paths, nlos: nlos_paths}}
    end
  end

  defp tier_suffix(height_m, mode) do
    cm = round(height_m * 100)

    case mode do
      :los -> "_h#{cm}cm"
      :nlos -> "_h#{cm}cm_nlos"
    end
  end

  defp resolve_antenna(%Coverage{antenna_slug: slug}) do
    case Antenna.get(slug) do
      nil -> {:error, {:unknown_antenna, slug}}
      antenna -> {:ok, antenna}
    end
  end

  defp resolve_location(coverage) do
    case Coverage.location(coverage) do
      nil -> {:error, :missing_location}
      coords -> {:ok, coords}
    end
  end

  defp compute_bbox(%Coverage{radius_m: radius}, lat, lon) do
    # Approximate a WGS84 bbox containing a circle of `radius` metres
    # around (lat, lon). 1° latitude ≈ 111 km; longitude scales with cos(lat).
    deg_lat = radius / 111_000.0
    deg_lon = radius / (111_000.0 * :math.cos(lat * :math.pi() / 180.0))
    {:ok, {lon - deg_lon, lat - deg_lat, lon + deg_lon, lat + deg_lat}}
  end

  defp fetch_terrain({_w, _s, _e, _n} = bbox, cell_size_m, lat) do
    # Convert metres to degrees at the given latitude (longitude axis is
    # narrowest, so use that for the worst-case spacing).
    cell_size_deg = cell_size_m / (111_000.0 * :math.cos(lat * :math.pi() / 180.0))

    case terrain_source().get_elevation_grid(bbox, cell_size_deg) do
      {:ok, grid} -> {:ok, grid}
      {:error, reason} -> {:error, {:terrain_unavailable, reason}}
    end
  end

  defp terrain_source do
    Application.get_env(:towerops, :coverage_terrain_module, Towerops.Lidar)
  end

  defp compute_pixels(coverage, ctx, lat, lon) do
    %{antenna: antenna, grid: grid, bbox: bbox} = ctx
    %{ncols: ncols, nrows: nrows} = grid
    {min_lon, min_lat, max_lon, max_lat} = bbox

    cell_lat = (max_lat - min_lat) / nrows
    cell_lon = (max_lon - min_lon) / ncols

    base_elev = elevation_at_or_zero(grid, lat, lon)
    tx_height = (coverage.height_agl_m || 0.0) + (coverage.height_above_rooftop_m || 0.0)
    tx_z = base_elev + tx_height
    eirp = Coverage.eirp_dbm(coverage, antenna.gain_dbi)

    pixel_ctx = Map.put(ctx, :eirp, eirp)
    tx = {{lat, lon}, tx_z, tx_height}

    log_every = max(1, div(nrows, 10))
    started = System.monotonic_time(:millisecond)

    rows =
      0..(nrows - 1)
      |> Task.async_stream(
        fn r ->
          row_lat = max_lat - (r + 0.5) * cell_lat

          for c <- 0..(ncols - 1) do
            pixel_lon = min_lon + (c + 0.5) * cell_lon
            compute_pixel(coverage, pixel_ctx, tx, {row_lat, pixel_lon})
          end
        end,
        ordered: true,
        max_concurrency: System.schedulers_online(),
        # 10-minute ceiling per row band — surfaces hangs instead of pinning a
        # worker forever. A real row should complete in ms.
        timeout: to_timeout(minute: 10),
        on_timeout: :kill_task
      )
      |> Stream.with_index()
      |> Stream.map(fn
        {{:ok, row}, idx} ->
          if rem(idx + 1, log_every) == 0 do
            elapsed = System.monotonic_time(:millisecond) - started

            Logger.info("CoverageWorker[#{coverage.id}]: pixel rows #{idx + 1}/#{nrows} (#{elapsed} ms elapsed)")
          end

          row

        {{:exit, reason}, idx} ->
          Logger.error("CoverageWorker[#{coverage.id}]: row #{idx} failed: #{inspect(reason)}; filling with NaN")

          List.duplicate(:nan, ncols)
      end)
      |> Enum.to_list()

    pixels = Enum.flat_map(rows, & &1)

    {:ok,
     %{
       pixels: pixels,
       dims: %{
         ncols: ncols,
         nrows: nrows,
         bbox: {min_lat, min_lon, max_lat, max_lon}
       }
     }}
  end

  defp compute_pixel(coverage, ctx, {tx_latlon, _tx_z, _tx_height} = tx, rx_latlon) do
    distance_m = Profile.great_circle_distance_m(tx_latlon, rx_latlon)

    cond do
      distance_m < 1.0 ->
        ctx.eirp

      distance_m > coverage.radius_m * 1.05 ->
        :nan

      true ->
        do_compute_pixel(coverage, ctx, tx, rx_latlon, distance_m)
    end
  end

  defp do_compute_pixel(
         coverage,
         ctx,
         {{tx_lat, tx_lon} = tx_latlon, tx_z, tx_height},
         {rx_lat, rx_lon} = rx_latlon,
         distance_m
       ) do
    %{antenna: antenna, grid: grid, clutter: clutter, eirp: eirp} = ctx
    canopy = Map.get(ctx, :canopy)

    # 1. Path profile (terrain + building/tree clutter from TX → RX).
    samples = max(3, min(@profile_samples, ceil(distance_m / coverage.cell_size_m) + 2))
    profile = Profile.sample(grid, tx_latlon, rx_latlon, samples, clutter: clutter, canopy: canopy)

    # Replace TX endpoint with the ground elevation under the antenna (so the
    # path loss model sees the correct antenna height above terrain).
    [_old_first | rest] = profile
    profile = [tx_z - tx_height | rest]

    # 2. Antenna pattern attenuation at this bearing/elevation.
    bearing_offset = bearing_offset_deg(tx_lat, tx_lon, rx_lat, rx_lon, coverage.azimuth_deg)
    elev_offset = elevation_offset_deg(tx_z, profile, coverage.downtilt_deg, distance_m)
    ant_atten = Antenna.attenuation_db(antenna, bearing_offset, elev_offset)

    # 3. Path loss over the profile.
    loss =
      Propagation.path_loss(
        distance_m,
        coverage.frequency_mhz * 1.0,
        profile,
        tx_height,
        coverage.receiver_height_m || 3.0
      )

    # 4. Foliage attenuation: linear scaling with the user's slider (0..100
    #    maps to 0..15 dB extra loss). A real model would use Weissberger.
    foliage_extra = (coverage.foliage_tuning || 0) * 0.15

    rssi = eirp - ant_atten - loss - foliage_extra

    if rssi < (coverage.rx_threshold_dbm || -90.0) do
      :nan
    else
      max(rssi, -130.0)
    end
  end

  defp elevation_at_or_zero(grid, lat, lon) do
    case Profile.elevation_at(grid, lat, lon) do
      nil -> 0.0
      elev -> elev
    end
  end

  # Bearing from TX to RX in degrees (0 = north, clockwise), then return
  # the offset relative to the antenna's boresight azimuth.
  defp bearing_offset_deg(tx_lat, tx_lon, rx_lat, rx_lon, azimuth_deg) do
    phi1 = tx_lat * :math.pi() / 180.0
    phi2 = rx_lat * :math.pi() / 180.0
    dlam = (rx_lon - tx_lon) * :math.pi() / 180.0

    y = :math.sin(dlam) * :math.cos(phi2)
    x = :math.cos(phi1) * :math.sin(phi2) - :math.sin(phi1) * :math.cos(phi2) * :math.cos(dlam)
    bearing = :math.atan2(y, x) * 180.0 / :math.pi()
    bearing = if bearing < 0, do: bearing + 360.0, else: bearing

    offset = bearing - azimuth_deg

    cond do
      offset > 180.0 -> offset - 360.0
      offset < -180.0 -> offset + 360.0
      true -> offset
    end
  end

  # Elevation angle from antenna to receiver, then offset relative to
  # the configured downtilt. Positive offset = above downtilt direction.
  defp elevation_offset_deg(_tx_z, profile, downtilt_deg, distance_m) do
    rx_z_terrain = List.last(profile) || 0.0
    [first | _] = profile
    delta_z = rx_z_terrain - first
    angle_deg = :math.atan2(delta_z, distance_m) * 180.0 / :math.pi()
    -downtilt_deg - angle_deg
  end

  defp update_progress(coverage, status, pct) do
    {:ok, updated} = Coverages.mark_status(coverage, status, %{progress_pct: pct})
    Coverages.broadcast(updated, {:coverage_status, status_atom(status), pct})
    updated
  end

  defp status_atom("draft"), do: :draft
  defp status_atom("queued"), do: :queued
  defp status_atom("computing"), do: :computing
  defp status_atom("ready"), do: :ready
  defp status_atom("failed"), do: :failed

  defp finalize(coverage, %{tiers: tiers, nlos_tiers: nlos_tiers, default: default}) do
    %{tif: tif, png: png, bbox: {min_lat, min_lon, max_lat, max_lon}} = default

    {:ok, updated} =
      Coverages.mark_status(coverage, "ready", %{
        progress_pct: 100,
        error_message: nil,
        computed_at: DateTime.truncate(DateTime.utc_now(), :second),
        raster_path: tif,
        png_path: png,
        height_tiers: tiers_map(tiers),
        nlos_tiers: tiers_map(nlos_tiers),
        bbox_min_lat: min_lat,
        bbox_min_lon: min_lon,
        bbox_max_lat: max_lat,
        bbox_max_lon: max_lon
      })

    Coverages.broadcast(updated, {:coverage_status, :ready, 100})
    log(coverage, "marked ready (raster_path=#{tif})")
    _ = Notifier.deliver_compute_complete(updated, :ready)
    :ok
  end

  defp tiers_map(tiers) do
    Map.new(tiers, fn {height_m, %{png: tier_png}} ->
      {Float.to_string(height_m * 1.0), tier_png}
    end)
  end

  defp fail(coverage, reason) do
    message = describe_error(reason)
    Logger.warning("CoverageWorker: failed for #{coverage.id}: #{message}")

    {:ok, updated} =
      Coverages.mark_status(coverage, "failed", %{
        error_message: message,
        progress_pct: 0
      })

    Coverages.broadcast(updated, {:coverage_status, :failed, message})
    _ = Notifier.deliver_compute_complete(updated, :failed)

    # Don't let Oban retry — the failure is on the coverage record.
    :ok
  end

  defp describe_error({:terrain_unavailable, :no_tile}),
    do:
      "No terrain data available for this site. The high-resolution LIDAR " <>
        "catalog covers Texas; the nationwide USGS NED (~10 m) fallback " <>
        "covers CONUS / AK / HI / PR. Sites outside those extents are not " <> "supported."

  defp describe_error({:terrain_unavailable, :grid_too_large}),
    do: "Coverage area is too large for the chosen cell size. Use a coarser cell or smaller radius."

  defp describe_error({:terrain_unavailable, :nodata}), do: "LIDAR tiles for this area returned no data."

  defp describe_error({:terrain_unavailable, {:gdal, _code, _output} = reason}) do
    Logger.warning("CoverageWorker: GDAL failure: #{inspect(reason)}")
    "Terrain data fetch failed. The TowerOps team has been notified — please retry shortly."
  end

  defp describe_error({:terrain_unavailable, reason}) do
    Logger.warning("CoverageWorker: terrain fetch failed: #{inspect(reason)}")
    "Terrain data is temporarily unavailable for this site. Please retry shortly."
  end

  defp describe_error({:unknown_antenna, slug}), do: "Antenna #{inspect(slug)} is not in the registry."

  defp describe_error(:missing_location), do: "Site has no coordinates set and no per-coverage override was supplied."

  defp describe_error({:gdal_translate_failed, code, output}),
    do: "GeoTIFF conversion failed (#{code}): #{String.slice(output, 0, 500)}"

  defp describe_error({:gdaldem_failed, code, output}),
    do: "PNG colorisation failed (#{code}): #{String.slice(output, 0, 500)}"

  defp describe_error(other), do: "Compute failed: #{inspect(other)}"
end