// Ported from SunCalc 2.0.1, BSD-2-Clause; see LICENSE and THIRD_PARTY.md.

///|
priv struct Equatorial {
  ra : Double
  dec : Double
}

// Espenak/Meeus piecewise TT-UT approximation, years 1900-2150.

///|
fn delta_t(d : Double) -> Double {
  let y = 2000.0 + d / 365.2425
  if y < 1920.0 {
    let t = y - 1900.0
    return -2.79 +
      t * (1.494119 + t * (-0.0598939 + t * (0.0061966 - t * 0.000197)))
  }
  if y < 1941.0 {
    let t = y - 1920.0
    return 21.20 + t * (0.84493 + t * (-0.076100 + t * 0.0020936))
  }
  if y < 1961.0 {
    let t = y - 1950.0
    return 29.07 + t * (0.407 + t * (-1.0 / 233.0 + t / 2547.0))
  }
  if y < 1986.0 {
    let t = y - 1975.0
    return 45.45 + t * (1.067 + t * (-1.0 / 260.0 - t / 718.0))
  }
  if y < 2005.0 {
    let t = y - 2000.0
    return 63.86 +
      t *
      (
        0.3345 +
        t *
        (-0.060374 + t * (0.0017275 + t * (0.000651814 + t * 0.00002373599)))
      )
  }
  if y < 2050.0 {
    let t = y - 2000.0
    return 62.92 + t * (0.32217 + t * 0.005589)
  }
  let t = (y - 1820.0) / 100.0
  -20.0 + 32.0 * t * t - 0.5628 * (2150.0 - y)
}

///|
fn to_tt(d : Double) -> Double {
  d + delta_t(d) / 86400.0
}

///|
fn azimuth(h : Double, phi : Double, dec : Double) -> Double {
  (
    @math.atan2(
      @math.sin(h),
      @math.cos(h) * @math.sin(phi) - @math.tan(dec) * @math.cos(phi),
    ) /
    Rad +
    540.0
  ) %
  360.0
}

///|
fn altitude(h : Double, phi : Double, dec : Double) -> Double {
  @math.asin(
    clamp_unit(
      @math.sin(phi) * @math.sin(dec) +
      @math.cos(phi) * @math.cos(dec) * @math.cos(h),
    ),
  )
}

///|
fn sidereal_time(d : Double, lw : Double) -> Double {
  Rad * (280.46061837 + 360.98564736629 * d) - lw
}

///|
fn refraction(alt : Double) -> Double {
  let h = if alt < 0.0 { 0.0 } else { alt }
  0.0002967 / @math.tan(h + 0.00312536 / (h + 0.08901179))
}

///|
fn sun_coords(d : Double) -> Equatorial {
  let t = d / 36525.0
  let l0 = Rad * (280.46646 + t * (36000.76983 + t * 0.0003032))
  let m = Rad * (357.52911 + t * (35999.05029 - t * 0.0001537))
  let sm = @math.sin(m)
  let cm = @math.cos(m)
  let c = Rad *
    (
      (1.914602 - t * (0.004817 + t * 0.000014)) * sm +
      (0.019993 - 0.000101 * t) * 2.0 * sm * cm +
      0.000289 * sm * (3.0 - 4.0 * sm * sm)
    )
  let om = Rad * (125.04 - 1934.136 * t)
  let l = l0 + c - Rad * (0.00569 + 0.00478 * @math.sin(om))
  let e = Rad *
    (23.439291 - t * (0.0130042 + t * (0.00000016 - t * 0.000000504))) +
    Rad * 0.00256 * @math.cos(om)
  {
    ra: @math.atan2(@math.cos(e) * @math.sin(l), @math.cos(l)),
    dec: @math.asin(@math.sin(e) * @math.sin(l)),
  }
}

///|
/// SunCalc v2 getPosition. Apparent refracted altitude; degrees, north clockwise.
pub fn sun_position(date : Instant, site : Observer) -> Position {
  let d = to_days(date)
  let c = sun_coords(to_tt(d))
  let h = sidereal_time(d, Rad * -site.lng) - c.ra
  let alt = altitude(h, Rad * site.lat, c.dec)
  {
    azimuth: azimuth(h, Rad * site.lat, c.dec),
    altitude: (alt + refraction(alt)) / Rad,
  }
}