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 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)
    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, 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_m, 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, tiles_dir) do
    {end_lat, end_lon} = destination_point(origin_lat, origin_lon, bearing, max_range_km)

    case Srtm.fetch_elevation_profile(origin_lat, origin_lon, end_lat, end_lon, tiles_dir) do
      {:ok, profile} ->
        result = analyse_ray(profile, max_range_km, freq_ghz, ant_height_m, ant_height_m)
        {reach_lat, reach_lon} = destination_point(origin_lat, origin_lon, bearing, result.reach_km)

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

      {:error, _} ->
        %{bearing: bearing, reach_km: max_range_km, lat: end_lat, lon: end_lon}
    end
  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