prop/lib/microwaveprop/moon.ex

defmodule Microwaveprop.Moon do
  @moduledoc """
  Low-precision lunar ephemeris.

  Uses the Astronomical Almanac "Low-Precision Formulas for the Moon"
  (Section D.4) — 6 longitude terms, 4 latitude terms, and 4 parallax
  terms. Accurate to ~0.3° on position and ~200 km on distance, which
  is well inside the beamwidth of any amateur EME dish at every
  supported microwave band (1° beam at 10 GHz for a 1 m dish; wider
  below that).

  Exposes two coordinates:

  - `geocentric/1` — the Moon's right ascension, declination, and
    distance from Earth's centre at a UTC instant.
  - `observer_position/3` — topocentric azimuth / elevation /
    slant range for an observer at `(lat, lon)` degrees.

  Both are pure functions of UTC plus observer coordinates; no side
  effects, no DB access. Callers drive them on a timer to refresh
  antenna aim points.
  """

  @earth_radius_km 6378.14

  @doc """
  Julian Day for a Gregorian `DateTime`. Caller is responsible for
  passing UT; TAI/UT1 offset (< 1 s) is ignored — the low-precision
  lunar series doesn't resolve it anyway.
  """
  @spec julian_day(DateTime.t()) :: float()
  def julian_day(%DateTime{} = dt) do
    day_frac = (dt.hour + dt.minute / 60 + dt.second / 3600) / 24

    {y, m} =
      if dt.month <= 2 do
        {dt.year - 1, dt.month + 12}
      else
        {dt.year, dt.month}
      end

    a = div(y, 100)
    b = 2 - a + div(a, 4)

    trunc(365.25 * (y + 4716)) + trunc(30.6001 * (m + 1)) + dt.day + day_frac + b - 1524.5
  end

  @doc """
  Geocentric ecliptic → equatorial position of the Moon at a UTC
  instant. Returns `%{ra, dec, distance_km}` where `ra` and `dec`
  are in degrees.
  """
  @spec geocentric(DateTime.t()) :: %{ra: float(), dec: float(), distance_km: float()}
  def geocentric(%DateTime{} = dt) do
    jd = julian_day(dt)
    t = (jd - 2_451_545.0) / 36_525.0

    # Mean-element arguments (degrees, wrapped before trig conversion).
    l = 218.32 + 481_267.881 * t
    m_prime = 134.9 + 477_198.85 * t
    m = 357.53 + 35_999.05 * t
    d = 297.85 + 445_267.11 * t
    f = 93.27 + 483_202.03 * t

    m_prime_r = deg_to_rad(m_prime)
    m_r = deg_to_rad(m)
    d_r = deg_to_rad(d)
    f_r = deg_to_rad(f)

    # Ecliptic longitude λ (degrees).
    lambda =
      l + 6.29 * :math.sin(m_prime_r) -
        1.27 * :math.sin(m_prime_r - 2 * d_r) +
        0.66 * :math.sin(2 * d_r) +
        0.21 * :math.sin(2 * m_prime_r) -
        0.19 * :math.sin(m_r) -
        0.11 * :math.sin(2 * f_r)

    # Ecliptic latitude β (degrees).
    beta =
      5.13 * :math.sin(f_r) +
        0.28 * :math.sin(m_prime_r + f_r) -
        0.28 * :math.sin(m_prime_r - f_r) -
        0.17 * :math.sin(2 * d_r - f_r)

    # Horizontal parallax π (degrees).
    parallax =
      0.9508 +
        0.0518 * :math.cos(m_prime_r) +
        0.0095 * :math.cos(m_prime_r - 2 * d_r) +
        0.0078 * :math.cos(2 * d_r) +
        0.0028 * :math.cos(2 * m_prime_r)

    distance_km = @earth_radius_km / :math.sin(deg_to_rad(parallax))

    # Mean obliquity of the ecliptic (IAU 1980).
    epsilon = 23.4393 - 0.0130 * t

    lambda_r = deg_to_rad(lambda)
    beta_r = deg_to_rad(beta)
    eps_r = deg_to_rad(epsilon)

    sin_dec =
      :math.sin(beta_r) * :math.cos(eps_r) +
        :math.cos(beta_r) * :math.sin(eps_r) * :math.sin(lambda_r)

    dec = rad_to_deg(:math.asin(sin_dec))

    ra =
      (:math.sin(lambda_r) * :math.cos(eps_r) - :math.tan(beta_r) * :math.sin(eps_r))
      |> :math.atan2(:math.cos(lambda_r))
      |> rad_to_deg()
      |> normalize_360()

    %{ra: ra, dec: dec, distance_km: distance_km}
  end

  @doc """
  Topocentric Moon position for an observer at `(lat, lon)` degrees,
  UTC `dt`. Returns azimuth ∈ [0, 360), elevation ∈ [-90, 90], and
  the geocentric slant range (the topocentric correction is below
  the module's stated accuracy floor).
  """
  @spec observer_position(DateTime.t(), float(), float()) ::
          %{azimuth: float(), elevation: float(), distance_km: float()}
  def observer_position(%DateTime{} = dt, lat_deg, lon_deg) do
    g = geocentric(dt)

    jd = julian_day(dt)
    t = (jd - 2_451_545.0) / 36_525.0

    # Greenwich Mean Sidereal Time — USNO polynomial (degrees).
    gmst =
      normalize_360(
        280.46061837 + 360.98564736629 * (jd - 2_451_545.0) + 0.000387933 * t * t -
          t * t * t / 38_710_000.0
      )

    lst = normalize_360(gmst + lon_deg)
    hour_angle = normalize_360(lst - g.ra)

    lat_r = deg_to_rad(lat_deg)
    dec_r = deg_to_rad(g.dec)
    ha_r = deg_to_rad(hour_angle)

    sin_el =
      :math.sin(lat_r) * :math.sin(dec_r) +
        :math.cos(lat_r) * :math.cos(dec_r) * :math.cos(ha_r)

    elevation = rad_to_deg(:math.asin(sin_el))

    azimuth =
      -:math.sin(ha_r)
      |> :math.atan2(:math.tan(dec_r) * :math.cos(lat_r) - :math.sin(lat_r) * :math.cos(ha_r))
      |> rad_to_deg()
      |> normalize_360()

    %{azimuth: azimuth, elevation: elevation, distance_km: g.distance_km}
  end

  @doc """
  Radial velocity of the Moon relative to an observer at `(lat, lon)`
  in km/s at UTC instant `dt`. Positive = the Moon is moving away
  (increasing slant range), negative = approaching.

  Computed via centred finite difference on the geocentric distance,
  which tracks the dominant rate (the Earth-Moon elliptical orbit
  plus the observer's diurnal baseline contribution is below 100 m/s).
  Accurate enough to size the ±kHz Doppler shift EME operators care
  about at the VHF/microwave bands this project supports.
  """
  @spec radial_velocity_km_s(DateTime.t(), float(), float()) :: float()
  def radial_velocity_km_s(%DateTime{} = dt, lat_deg, lon_deg) do
    before = DateTime.shift(dt, second: -1)
    later = DateTime.shift(dt, second: 1)

    d_before = observer_position(before, lat_deg, lon_deg).distance_km
    d_later = observer_position(later, lat_deg, lon_deg).distance_km

    # Centred difference in km per 2 s → km/s.
    (d_later - d_before) / 2.0
  end

  # --- Helpers --------------------------------------------------------

  defp deg_to_rad(d), do: d * :math.pi() / 180.0
  defp rad_to_deg(r), do: r * 180.0 / :math.pi()

  defp normalize_360(x) do
    r = :math.fmod(x, 360.0)
    if r < 0, do: r + 360.0, else: r
  end
end