defmodule Microwaveprop.Propagation.MoonEphemeris do
  @moduledoc """
  Approximate moon-position calculation for EME classification.

  Accuracy target: ±1° — enough to decide "is the moon above the horizon
  for both stations right now?" which is the EME feasibility test. We
  use Jean Meeus's low-precision algorithm (Astronomical Algorithms,
  chapter 47), good to ~10' over historical amateur contact timestamps.

  No external ephemeris table — pure math so it works for any date.
  """

  @deg_to_rad :math.pi() / 180.0
  @rad_to_deg 180.0 / :math.pi()

  @doc """
  Returns the moon's altitude above the horizon in degrees for a given
  observer lat/lon and UTC timestamp. Negative means the moon is below
  the horizon (not reachable).
  """
  @spec altitude_deg(number(), number(), DateTime.t()) :: float()
  def altitude_deg(lat, lon, %DateTime{} = ts) do
    {ra, dec} = moon_ra_dec(ts)
    gst = greenwich_sidereal_hours(ts)
    # Local hour angle in degrees
    lst_deg = :math.fmod(gst * 15.0 + lon, 360.0)
    hour_angle = lst_deg - ra * 15.0

    lat_r = lat * @deg_to_rad
    dec_r = dec * @deg_to_rad
    ha_r = hour_angle * @deg_to_rad

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

    sin_alt |> max(-1.0) |> min(1.0) |> :math.asin() |> Kernel.*(@rad_to_deg)
  end

  @doc """
  `true` if the moon is above the horizon at the given observer + time.
  """
  @spec above_horizon?(number(), number(), DateTime.t()) :: boolean()
  def above_horizon?(lat, lon, %DateTime{} = ts), do: altitude_deg(lat, lon, ts) > 0.0

  # — private math helpers below —

  # Returns {right ascension (hours), declination (deg)} of the moon.
  defp moon_ra_dec(%DateTime{} = ts) do
    jd = julian_day(ts)
    t = (jd - 2_451_545.0) / 36_525.0

    # Mean orbital elements (Meeus §47, degrees).
    l_prime = 218.3164477 + 481_267.88123421 * t
    d = 297.8501921 + 445_267.1114034 * t
    m = 357.5291092 + 35_999.0502909 * t
    m_prime = 134.9633964 + 477_198.8675055 * t
    f = 93.2720950 + 483_202.0175233 * t

    # Longitude correction (leading term only — 6.289°).
    lambda =
      l_prime +
        6.289 * sin_deg(m_prime) -
        1.274 * sin_deg(2 * d - m_prime) +
        0.658 * sin_deg(2 * d) -
        0.186 * sin_deg(m)

    # Latitude (β) leading term.
    beta = 5.128 * sin_deg(f)

    # Obliquity of the ecliptic.
    eps = 23.439 - 0.0000004 * (jd - 2_451_545.0)

    # Convert ecliptic (λ, β) → equatorial (α, δ).
    lam_r = lambda * @deg_to_rad
    beta_r = beta * @deg_to_rad
    eps_r = eps * @deg_to_rad

    sin_alpha = :math.sin(lam_r) * :math.cos(eps_r) - :math.tan(beta_r) * :math.sin(eps_r)
    cos_alpha = :math.cos(lam_r)
    alpha_r = :math.atan2(sin_alpha, cos_alpha)

    sin_delta =
      :math.sin(beta_r) * :math.cos(eps_r) +
        :math.cos(beta_r) * :math.sin(eps_r) * :math.sin(lam_r)

    delta_r = :math.asin(sin_delta)

    ra_hours = :math.fmod(alpha_r * @rad_to_deg / 15.0 + 48.0, 24.0)
    dec_deg = delta_r * @rad_to_deg
    {ra_hours, dec_deg}
  end

  defp julian_day(%DateTime{year: y, month: m, day: d, hour: h, minute: mi, second: s}) do
    {y2, m2} =
      if m <= 2 do
        {y - 1, m + 12}
      else
        {y, m}
      end

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

    jd_day =
      Float.floor(365.25 * (y2 + 4716)) +
        Float.floor(30.6001 * (m2 + 1)) +
        d + b - 1524.5

    frac = (h + mi / 60 + s / 3600) / 24.0
    jd_day + frac
  end

  defp greenwich_sidereal_hours(%DateTime{} = ts) do
    jd = julian_day(ts)
    t = (jd - 2_451_545.0) / 36_525.0

    # Mean sidereal time at Greenwich, in hours.
    gmst_sec =
      67_310.54841 +
        (876_600 * 3600 + 8_640_184.812866) * t +
        0.093104 * t * t -
        6.2e-6 * t * t * t

    hours = gmst_sec / 3600.0
    :math.fmod(:math.fmod(hours, 24.0) + 24.0, 24.0)
  end

  defp sin_deg(deg), do: :math.sin(deg * @deg_to_rad)
end