// Ported from SunCalc 2.0.1, BSD-2-Clause. 60+60 Meeus periodic terms.

///|
const EarthRadius : Double = 6378.14

///|
priv struct MoonCoords {
  ra : Double
  dec : Double
  dist : Double
}

///|
pub(all) struct MoonPosition {
  azimuth : Double
  altitude : Double
  distance : Double
  parallactic_angle : Double
} derive(Eq, Debug)

///|
fn nutation_obliquity(t : Double) -> (Double, Double) {
  let om = Rad * (125.04452 - 1934.136261 * t)
  let ls = Rad * (280.4665 + 36000.7698 * t)
  let lm = Rad * (218.3165 + 481267.8813 * t)
  let dpsi = (
      -17.20 * @math.sin(om) -
      1.32 * @math.sin(2.0 * ls) -
      0.23 * @math.sin(2.0 * lm) +
      0.21 * @math.sin(2.0 * om)
    ) /
    3600.0
  let deps = (
      9.20 * @math.cos(om) +
      0.57 * @math.cos(2.0 * ls) +
      0.10 * @math.cos(2.0 * lm) -
      0.09 * @math.cos(2.0 * om)
    ) /
    3600.0
  let eps0 = 23.439291 - t * (0.0130042 + t * (0.00000016 - t * 0.000000504))
  (dpsi, Rad * (eps0 + deps))
}

///|
fn moon_coords(d : Double) -> MoonCoords {
  let t = d / 36525.0
  let lp = 218.3164477 +
    t *
    (481267.88123421 + t * (-0.0015786 + t * (1.0 / 538841.0 - t / 65194000.0)))
  let dd = 297.8501921 +
    t *
    (445267.1114034 + t * (-0.0018819 + t * (1.0 / 545868.0 - t / 113065000.0)))
  let m = 357.5291092 + t * (35999.0502909 + t * (-0.0001536 + t / 24490000.0))
  let mp = 134.9633964 +
    t *
    (477198.8675055 + t * (0.0087414 + t * (1.0 / 69699.0 - t / 14712000.0)))
  let f = 93.2720950 +
    t *
    (
      483202.0175233 +
      t * (-0.0036539 + t * (-1.0 / 3526000.0 + t / 863310000.0))
    )
  let a1 = 119.75 + 131.849 * t
  let a2 = 53.09 + 479264.290 * t
  let a3 = 313.45 + 481266.484 * t
  let e = 1.0 - t * (0.002516 + t * 0.0000074)
  let dr = Rad * dd
  let mr = Rad * m
  let mpr = Rad * mp
  let fr = Rad * f
  let mut sl = 0.0
  let mut sr = 0.0
  let mut sb = 0.0
  for i = 0; i < moon_lon.length(); i = i + 6 {
    let k = moon_lon[i + 1]
    let arg = moon_lon[i].to_double() * dr +
      k.to_double() * mr +
      moon_lon[i + 2].to_double() * mpr +
      moon_lon[i + 3].to_double() * fr
    let ef = if k == 1 || k == -1 {
      e
    } else if k == 2 || k == -2 {
      e * e
    } else {
      1.0
    }
    sl += moon_lon[i + 4].to_double() * ef * @math.sin(arg)
    sr += moon_lon[i + 5].to_double() * ef * @math.cos(arg)
  }
  for i = 0; i < moon_lat.length(); i = i + 5 {
    let k = moon_lat[i + 1]
    let arg = moon_lat[i].to_double() * dr +
      k.to_double() * mr +
      moon_lat[i + 2].to_double() * mpr +
      moon_lat[i + 3].to_double() * fr
    let ef = if k == 1 || k == -1 {
      e
    } else if k == 2 || k == -2 {
      e * e
    } else {
      1.0
    }
    sb += moon_lat[i + 4].to_double() * ef * @math.sin(arg)
  }
  let a1r = Rad * a1
  let lpr = Rad * lp
  sl += 3958.0 * @math.sin(a1r) +
    1962.0 * @math.sin(lpr - fr) +
    318.0 * @math.sin(Rad * a2)
  sb += -2235.0 * @math.sin(lpr) +
    382.0 * @math.sin(Rad * a3) +
    175.0 * @math.sin(a1r - fr) +
    175.0 * @math.sin(a1r + fr) +
    127.0 * @math.sin(lpr - mpr) -
    115.0 * @math.sin(lpr + mpr)
  let (dpsi, eps) = nutation_obliquity(t)
  let l = Rad * (lp + sl / 1000000.0 + dpsi)
  let b = Rad * (sb / 1000000.0)
  {
    ra: @math.atan2(
      @math.sin(l) * @math.cos(eps) - @math.tan(b) * @math.sin(eps),
      @math.cos(l),
    ),
    dec: @math.asin(
      @math.sin(b) * @math.cos(eps) +
      @math.cos(b) * @math.sin(eps) * @math.sin(l),
    ),
    dist: 385000.56 + sr / 1000.0,
  }
}

///|
/// Topocentric apparent Moon position; distance is geocentric kilometres.
pub fn moon_position(date : Instant, site : Observer) -> MoonPosition {
  let phi = Rad * site.lat
  let d = to_days(date)
  let c = moon_coords(to_tt(d))
  let h_angle = sidereal_time(d, Rad * -site.lng) - c.ra
  let h_geo = altitude(h_angle, phi, c.dec)
  let h = h_geo - @math.asin(EarthRadius / c.dist * @math.cos(h_geo))
  let pa = @math.atan2(
    @math.sin(h_angle),
    @math.tan(phi) * @math.cos(c.dec) - @math.sin(c.dec) * @math.cos(h_angle),
  )
  {
    azimuth: azimuth(h_angle, phi, c.dec),
    altitude: (h + refraction(h)) / Rad,
    distance: c.dist,
    parallactic_angle: pa / Rad,
  }
}