prop/lib/microwaveprop/terrain/viewshed.ex

defmodule Microwaveprop.Terrain.Viewshed do
  @moduledoc false

  alias Microwaveprop.Terrain.Srtm
  alias Microwaveprop.Terrain.TerrainAnalysis

  @earth_radius_km 6371.0
  @default_angular_step 2
  @default_max_range_km 50

  @doc "Compute destination lat/lon given origin, bearing (degrees), and distance (km)."
  def destination_point(lat, lon, bearing_deg, dist_km) do
    lat_rad = deg_to_rad(lat)
    lon_rad = deg_to_rad(lon)
    brg_rad = deg_to_rad(bearing_deg)
    d = dist_km / @earth_radius_km

    lat2 =
      :math.asin(
        :math.sin(lat_rad) * :math.cos(d) +
          :math.cos(lat_rad) * :math.sin(d) * :math.cos(brg_rad)
      )

    lon2 =
      lon_rad +
        :math.atan2(
          :math.sin(brg_rad) * :math.sin(d) * :math.cos(lat_rad),
          :math.cos(d) - :math.sin(lat_rad) * :math.sin(lat2)
        )

    {rad_to_deg(lat2), rad_to_deg(lon2)}
  end

  @doc """
  Find the max clear distance along a ray from TerrainAnalysis points.
  Skips endpoints (first/last), returns the dist_km of the last clear
  interior point before the first obstruction.
  """
  def find_reach_km(points, max_range_km) do
    interior = Enum.slice(points, 1..-2//1)

    case find_first_obstructed_index(interior) do
      nil ->
        max_range_km

      0 ->
        0.0

      idx ->
        Enum.at(interior, idx - 1).dist_km
    end
  end

  @doc """
  Compute effective reach accounting for both terrain and propagation conditions.
  Higher propagation scores indicate ducting potential, which lets signals
  propagate beyond terrain obstructions via atmospheric waveguide.
  """
  def effective_reach_km(analysis, max_range_km, score) do
    case analysis.verdict do
      "CLEAR" ->
        max_range_km

      "FRESNEL_MINOR" ->
        max_range_km * 0.9

      "FRESNEL_PARTIAL" ->
        max_range_km * 0.7

      "BLOCKED" ->
        terrain_factor = terrain_reach_factor(analysis.diffraction_db)
        ducting_factor = ducting_reach_factor(score)
        max_range_km * max(terrain_factor, ducting_factor)
    end
  end

  # Convert diffraction loss to a range reduction factor (0.0–1.0).
  # Mild diffraction still allows significant propagation;
  # heavy blockage attenuates sharply.
  defp terrain_reach_factor(db) when db <= 3, do: 0.8
  defp terrain_reach_factor(db) when db <= 6, do: 0.5
  defp terrain_reach_factor(db) when db <= 12, do: 0.3
  defp terrain_reach_factor(db) when db <= 20, do: 0.15
  defp terrain_reach_factor(_db), do: 0.05

  # Higher propagation scores mean stronger ducting potential —
  # signals can ride atmospheric layers over terrain obstacles.
  # Returns a floor factor so blocked paths still get range
  # proportional to atmospheric conditions.
  defp ducting_reach_factor(score) when score >= 80, do: 0.7
  defp ducting_reach_factor(score) when score >= 65, do: 0.5
  defp ducting_reach_factor(score) when score >= 50, do: 0.3
  defp ducting_reach_factor(score) when score >= 33, do: 0.15
  defp ducting_reach_factor(_score), do: 0.05

  @doc """
  Compute a terrain viewshed from a point. Returns a map with :origin
  and :boundary (list of %{bearing, reach_km, lat, lon}).

  Options:
    - :freq_ghz — frequency for Fresnel zone calc (default 10.0)
    - :max_range_km — max ray distance (default 50)
    - :ant_height_m — antenna height at both ends (default 2.4)
    - :angular_step — degrees between rays (default 2)
    - :tiles_dir — SRTM tiles directory (default from config)
  """
  def compute(lat, lon, opts \\ []) do
    freq_ghz = Keyword.get(opts, :freq_ghz, 10.0)
    max_range_km = Keyword.get(opts, :max_range_km, @default_max_range_km)
    ant_height_m = Keyword.get(opts, :ant_height_m, 2.4)
    score = Keyword.get(opts, :score, 50)
    angular_step = Keyword.get(opts, :angular_step, @default_angular_step)
    tiles_dir = Keyword.get(opts, :tiles_dir, srtm_tiles_dir())

    bearings = Enum.to_list(0..359//angular_step)

    boundary =
      bearings
      |> Task.async_stream(
        fn bearing ->
          compute_ray(lat, lon, bearing, freq_ghz, max_range_km, ant_height_m, score, tiles_dir)
        end,
        max_concurrency: System.schedulers_online(),
        timeout: 30_000,
        on_timeout: :kill_task
      )
      |> Enum.map(fn
        {:ok, result} -> result
        {:exit, _} -> nil
      end)
      |> Enum.reject(&is_nil/1)

    %{origin: %{lat: lat, lon: lon}, boundary: boundary}
  end

  @doc """
  Analyse a single ray's profile. Public for testing.
  Returns %{reach_km: float, verdict: string}.
  """
  def analyse_ray(profile, dist_km, freq_ghz, ant_ht_a_m, ant_ht_b_m) do
    analysis = TerrainAnalysis.analyse(profile, dist_km, freq_ghz, ant_ht_a: ant_ht_a_m, ant_ht_b: ant_ht_b_m)
    reach_km = find_reach_km(analysis.points, dist_km)
    %{reach_km: reach_km, verdict: analysis.verdict}
  end

  defp compute_ray(origin_lat, origin_lon, bearing, freq_ghz, max_range_km, ant_height_m, score, tiles_dir) do
    # Check terrain within the radio horizon where terrain features matter.
    # Beyond this, propagation is atmospheric (ducting/scatter) and terrain
    # doesn't block — earth curvature is handled by the atmosphere.
    terrain_check_km = min(radio_horizon_km(ant_height_m) * 2.0, max_range_km)
    {end_lat, end_lon} = destination_point(origin_lat, origin_lon, bearing, terrain_check_km)

    case Srtm.fetch_elevation_profile(origin_lat, origin_lon, end_lat, end_lon, tiles_dir) do
      {:ok, profile} ->
        analysis =
          TerrainAnalysis.analyse(profile, terrain_check_km, freq_ghz, ant_ht_a: ant_height_m, ant_ht_b: ant_height_m)

        reach_km = effective_reach_km(analysis, max_range_km, score)
        {reach_lat, reach_lon} = destination_point(origin_lat, origin_lon, bearing, reach_km)

        %{
          bearing: bearing,
          reach_km: reach_km,
          lat: reach_lat,
          lon: reach_lon
        }

      {:error, _} ->
        {end_lat2, end_lon2} = destination_point(origin_lat, origin_lon, bearing, max_range_km)
        %{bearing: bearing, reach_km: max_range_km, lat: end_lat2, lon: end_lon2}
    end
  end

  # Radio horizon distance in km for a given antenna height (K=4/3 atmosphere)
  defp radio_horizon_km(height_m) do
    :math.sqrt(2 * (4 / 3) * @earth_radius_km * 1000 * height_m) / 1000
  end

  defp srtm_tiles_dir do
    Application.get_env(:microwaveprop, :srtm_tiles_dir, Path.expand("~/srtm/tiles"))
  end

  defp find_first_obstructed_index(interior) do
    Enum.find_index(interior, & &1.obstructed)
  end

  defp deg_to_rad(deg), do: deg * :math.pi() / 180
  defp rad_to_deg(rad), do: rad * 180 / :math.pi()
end