defmodule Microwaveprop.Geo do
  @moduledoc "Geodetic helper functions: distances, bearings, coordinate math."

  alias Microwaveprop.Radio.Maidenhead

  @earth_radius_km 6371.0

  @doc "Great-circle distance between two points in km."
  @spec haversine_km(number(), number(), number(), number()) :: float()
  def haversine_km(lat1, lon1, lat2, lon2) do
    dlat = deg_to_rad(lat2 - lat1)
    dlon = deg_to_rad(lon2 - lon1)
    rlat1 = deg_to_rad(lat1)
    rlat2 = deg_to_rad(lat2)

    a =
      :math.sin(dlat / 2) ** 2 +
        :math.cos(rlat1) * :math.cos(rlat2) * :math.sin(dlon / 2) ** 2

    # `atan2(√a, √(1−a))` is the numerically stable form: when `a`
    # rounds to slightly above 1 the asin form returns NaN, while
    # atan2 stays finite.
    2 * @earth_radius_km * :math.atan2(:math.sqrt(a), :math.sqrt(1 - a))
  end

  @doc "Initial bearing from point 1 to point 2 in degrees (0-360)."
  @spec bearing_deg(number(), number(), number(), number()) :: float()
  def bearing_deg(lat1, lon1, lat2, lon2) do
    lat1r = deg_to_rad(lat1)
    lat2r = deg_to_rad(lat2)
    dlonr = deg_to_rad(lon2 - lon1)
    y = :math.sin(dlonr) * :math.cos(lat2r)
    x = :math.cos(lat1r) * :math.sin(lat2r) - :math.sin(lat1r) * :math.cos(lat2r) * :math.cos(dlonr)
    b = :math.atan2(y, x) * 180 / :math.pi()
    Float.round(:math.fmod(b + 360, 360), 1)
  end

  @doc "Center of a 4-character Maidenhead grid square."
  @spec maidenhead_center(String.t()) :: {float(), float()} | nil
  def maidenhead_center(grid) when is_binary(grid) and byte_size(grid) >= 4 do
    case Maidenhead.to_latlon(grid) do
      {:ok, {lat, lon}} -> {lat, lon}
      :error -> nil
    end
  end

  @doc "4-character Maidenhead grid for a lat/lon."
  @spec latlon_to_grid4(number(), number()) :: String.t()
  def latlon_to_grid4(lat, lon) do
    lon_field = trunc((lon + 180) / 20)
    lat_field = trunc((lat + 90) / 10)
    lon_sq = trunc(rem(trunc(lon + 180), 20) / 2)
    lat_sq = trunc(rem(trunc(lat + 90), 10) / 1)

    <<?A + lon_field, ?A + lat_field, ?0 + lon_sq, ?0 + lat_sq>>
  end

  defp deg_to_rad(d), do: d * :math.pi() / 180
end