--- Time and astronomy calculations.
-- Ported from "Calendrical Calculations" (4th edition)
-- by Nachum Dershowitz and Edward M. Reingold.
-- Original Lisp code (CALENDRICA 4.0) is Apache 2.0 licensed.
-- @module calendrica-astro
-- @release 0.1 2026-07-19

local M = {}

local basic     = require("calendrica-basic")
local gregorian = require("calendrica-gregorian")


-- === Time unit helpers ===

--- `x` hours as a fraction of a day.
-- @tparam number x Hours.
-- @treturn number Fraction of a day.
function M.hr(x)
  return x / 24
end
local hr = M.hr

--- `x` minutes as a fraction of a day.
-- @tparam number x Minutes.
-- @treturn number Fraction of a day.
function M.mn(x)
  return x / 24 / 60
end
local mn = M.mn

--- `x` seconds as a fraction of a day.
-- @tparam number x Seconds.
-- @treturn number Fraction of a day.
function M.sec(x)
  return x / 24 / 60 / 60
end
local sec = M.sec

--- `x` meters (identity; for type clarity).
-- @tparam number x Meters.
-- @treturn number Distance in meters.
function M.mt(x)
  return x
end
local mt = M.mt

--- `x` degrees (identity; for type clarity).
-- @tparam number x Degrees.
-- @treturn number Angle in degrees.
function M.deg(x)
  return x
end

--- `x` arcminutes expressed in degrees.
-- @tparam number x Arcminutes.
-- @treturn number Degrees.
function M.mins(x)
  return x / 60
end
local mins = M.mins

--- `x` arcseconds expressed in degrees.
-- @tparam number x Arcseconds.
-- @treturn number Degrees.
function M.secs(x)
  return x / 3600
end
local secs = M.secs

--- Angle from `d` degrees, `m` arcminutes, `s` arcseconds.
-- @tparam number d Degrees.
-- @tparam number m Arcminutes.
-- @tparam number s Arcseconds.
-- @treturn number Angle in degrees.
function M.angle(d, m, s)
  return d + (m + s / 60) / 60
end
local angle = M.angle


-- === Trigonometry in degrees ===

-- Convert angle theta from radians to degrees.
local function degrees_from_radians(theta)
  return (theta / math.pi / (1/180)) % 360
end

-- Convert angle theta from degrees to radians.
local function radians_from_degrees(theta)
  return (theta % 360) * math.pi * (1/180)
end

--- Sine of `theta` (given in degrees).
-- @tparam number theta Angle in degrees.
-- @treturn number Sine value.
function M.sin_degrees(theta)
  return math.sin(radians_from_degrees(theta))
end
local sin_degrees = M.sin_degrees

--- Cosine of `theta` (given in degrees).
-- @tparam number theta Angle in degrees.
-- @treturn number Cosine value.
function M.cos_degrees(theta)
  return math.cos(radians_from_degrees(theta))
end
local cos_degrees = M.cos_degrees

--- Tangent of `theta` (given in degrees).
-- @tparam number theta Angle in degrees.
-- @treturn number Tangent value.
function M.tan_degrees(theta)
  return math.tan(radians_from_degrees(theta))
end
local tan_degrees = M.tan_degrees

--- Arctangent of `y`/`x` in degrees; returns BOGUS if both are 0.
-- @tparam number y Numerator.
-- @tparam number x Denominator.
-- @treturn number Angle in degrees, or BOGUS.
function M.arctan_degrees(y, x)
  if x == 0 and y == 0 then
    return basic.BOGUS
  end
  if x == 0 then
    return basic.sign(y) * 90
  end
  local alpha = degrees_from_radians(math.atan(y / x))
  if x >= 0 then
    return alpha % 360
  else
    return (alpha + 180) % 360
  end
end
local arctan_degrees = M.arctan_degrees

--- Arcsine of `x` in degrees.
-- @tparam number x Sine value (−1..1).
-- @treturn number Angle in degrees.
function M.arcsin_degrees(x)
  return degrees_from_radians(math.asin(x))
end
local arcsin_degrees = M.arcsin_degrees

--- Arccosine of `x` in degrees.
-- @tparam number x Cosine value (−1..1).
-- @treturn number Angle in degrees.
function M.arccos_degrees(x)
  return degrees_from_radians(math.acos(x))
end
local arccos_degrees = M.arccos_degrees


-- === Location ===

--- Construct a location from latitude, longitude, elevation, and time zone.
-- @tparam number latitude  Degrees north.
-- @tparam number longitude Degrees east.
-- @tparam number elevation Meters above sea level.
-- @tparam number zone      Hours offset from UT.
-- @treturn table {latitude, longitude, elevation, zone}
function M.location(latitude, longitude, elevation, zone)
  return {latitude, longitude, elevation, zone}
end
local location = M.location

--- Latitude (degrees north) of location `loc`.
-- @tparam table loc Location.
-- @treturn number Latitude.
function M.latitude(loc)   return loc[1] end
local latitude = M.latitude

--- Longitude (degrees east) of location `loc`.
-- @tparam table loc Location.
-- @treturn number Longitude.
function M.longitude(loc)  return loc[2] end
local longitude = M.longitude

--- Elevation (meters) of location `loc`.
-- @tparam table loc Location.
-- @treturn number Elevation.
function M.elevation(loc)  return loc[3] end
local elevation = M.elevation

--- Time zone offset (hours from UT) of location `loc`.
-- @tparam table loc Location.
-- @treturn number Zone offset.
function M.zone(loc)       return loc[4] end
local zone = M.zone

local MECCA = location(angle(21, 25, 24), angle(39, 49, 24), mt(298), hr(3))

local PADUA = location(angle(45, 24, 28), angle(11, 53, 9), mt(18), hr(1))

--- Angle (clockwise from North) to face `focus` when standing in `loc`.
-- Subject to errors near focus and its antipode.
-- @tparam table loc   Observer's location.
-- @tparam table focus Target location.
-- @treturn number Bearing in degrees.
function M.direction(loc, focus)
  local phi       = latitude(loc)
  local phi_prime = latitude(focus)
  local psi       = longitude(loc)
  local psi_prime = longitude(focus)
  local y = sin_degrees(psi_prime - psi)
  local x = cos_degrees(phi) * tan_degrees(phi_prime)
          - sin_degrees(phi) * cos_degrees(psi - psi_prime)
  if (x == 0 and y == 0) or phi_prime == 90 then
    return 0
  elseif phi_prime == -90 then
    return 180
  else
    return arctan_degrees(y, x)
  end
end


-- === Time conversions ===

--- Standard time from universal time `tee_u` at `loc`.
-- @tparam number tee_u Universal time moment.
-- @tparam table  loc   Location.
-- @treturn number Standard time moment.
function M.standard_from_universal(tee_u, loc)
  return tee_u + zone(loc)
end
local standard_from_universal = M.standard_from_universal

--- Universal time from standard time `tee_s` at `loc`.
-- @tparam number tee_s Standard time moment.
-- @tparam table  loc   Location.
-- @treturn number Universal time moment.
function M.universal_from_standard(tee_s, loc)
  return tee_s - zone(loc)
end
local universal_from_standard = M.universal_from_standard

-- Difference between UT and local mean time at longitude phi.
local function zone_from_longitude(phi)
  return phi / 360
end

--- Local time from universal time `tee_u` at `loc`.
-- @tparam number tee_u Universal time moment.
-- @tparam table  loc   Location.
-- @treturn number Local time moment.
function M.local_from_universal(tee_u, loc)
  return tee_u + zone_from_longitude(longitude(loc))
end
local local_from_universal = M.local_from_universal

--- Universal time from local time `tee_l` at `loc`.
-- @tparam number tee_l Local time moment.
-- @tparam table  loc   Location.
-- @treturn number Universal time moment.
function M.universal_from_local(tee_l, loc)
  return tee_l - zone_from_longitude(longitude(loc))
end
local universal_from_local = M.universal_from_local

local equation_of_time  -- forward declaration (used by apparent_from_local/local_from_apparent before definition)

-- Standard time from local time tee_l at loc.
local function standard_from_local(tee_l, loc)
  return standard_from_universal(universal_from_local(tee_l, loc), loc)
end

-- Local time from standard time tee_s at loc.
local function local_from_standard(tee_s, loc)
  return local_from_universal(universal_from_standard(tee_s, loc), loc)
end

-- Sundial time from local time tee_l at loc.
local function apparent_from_local(tee_l, loc)
  return tee_l + equation_of_time(universal_from_local(tee_l, loc))
end

-- Local time from sundial time tee at loc.
local function local_from_apparent(tee, loc)
  return tee - equation_of_time(universal_from_local(tee, loc))
end

--- True (apparent) time at universal time `tee_u` at `loc`.
-- @tparam number tee_u Universal time moment.
-- @tparam table  loc   Location.
-- @treturn number Apparent time moment.
function M.apparent_from_universal(tee_u, loc)
  return apparent_from_local(local_from_universal(tee_u, loc), loc)
end

-- Universal time from sundial time tee at loc.
local function universal_from_apparent(tee, loc)
  return universal_from_local(local_from_apparent(tee, loc), loc)
end

--- Universal time of true (apparent) midnight of fixed `date` at `loc`.
-- @tparam number date Fixed date.
-- @tparam table  loc  Location.
-- @treturn number Universal time moment.
function M.midnight(date, loc)
  return universal_from_apparent(date, loc)
end

--- Universal time of midday on fixed `date` at `loc`.
-- @tparam number date Fixed date.
-- @tparam table  loc  Location.
-- @treturn number Universal time moment.
function M.midday(date, loc)
  return universal_from_apparent(date + hr(12), loc)
end
local midday = M.midday


-- === Dynamical time ===

-- Noon at start of Gregorian year 2000.
local J2000 = hr(12) + gregorian.gregorian_new_year(2000)

local julian_centuries   -- forward declaration (used by obliquity etc. before definition)

--- Ephemeris correction (dynamical − universal time) at moment `tee`.
-- @tparam number tee Moment.
-- @treturn number Correction in days.
function M.ephemeris_correction(tee)
  local year = gregorian.gregorian_year_from_fixed(math.floor(tee))
  local c = gregorian.gregorian_date_difference(
    gregorian.gregorian_date(1900, gregorian.JANUARY, 1),
    gregorian.gregorian_date(year, gregorian.JULY, 1)
  ) / 36525
  local y2000  = year - 2000
  local y1700  = year - 1700
  local y1600  = year - 1600
  local y1000  = (year - 1000) / 100
  local y0     = year / 100
  local y1820  = (year - 1820) / 100
  local c2051  = (1/86400) * (-20 + 32 * ((year - 1820) / 100)^2
                   + 0.5628 * (2150 - year))
  local c2006  = (1/86400) * basic.poly(y2000, {62.92, 0.32217, 0.005589})
  local c1987  = (1/86400) * basic.poly(y2000,
                   {63.86, 0.3345, -0.060374, 0.0017275,
                    0.000651814, 0.00002373599})
  local c1900  = basic.poly(c,
                   {-0.00002, 0.000297, 0.025184, -0.181133,
                    0.553040, -0.861938, 0.677066, -0.212591})
  local c1800  = basic.poly(c,
                   {-0.000009, 0.003844, 0.083563, 0.865736,
                    4.867575, 15.845535, 31.332267, 38.291999,
                    28.316289, 11.636204, 2.043794})
  local c1700  = (1/86400) * basic.poly(y1700,
                   {8.118780842, -0.005092142, 0.003336121, -0.0000266484})
  local c1600  = (1/86400) * basic.poly(y1600,
                   {120, -0.9808, -0.01532, 0.000140272128})
  local c500   = (1/86400) * basic.poly(y1000,
                   {1574.2, -556.01, 71.23472, 0.319781,
                    -0.8503463, -0.005050998, 0.0083572073})
  local c0     = (1/86400) * basic.poly(y0,
                   {10583.6, -1014.41, 33.78311, -5.952053,
                    -0.1798452, 0.022174192, 0.0090316521})
  local other  = (1/86400) * basic.poly(y1820, {-20, 0, 32})
  if     year >= 2051 and year <= 2150 then return c2051
  elseif year >= 2006 and year <= 2050 then return c2006
  elseif year >= 1987 and year <= 2005 then return c1987
  elseif year >= 1900 and year <= 1986 then return c1900
  elseif year >= 1800 and year <= 1899 then return c1800
  elseif year >= 1700 and year <= 1799 then return c1700
  elseif year >= 1600 and year <= 1699 then return c1600
  elseif year >= 500  and year <= 1599 then return c500
  elseif year >  -500 and year < 500   then return c0
  else                                       return other
  end
end
local ephemeris_correction = M.ephemeris_correction

-- Universal moment from Dynamical time tee.
local function universal_from_dynamical(tee)
  return tee - ephemeris_correction(tee)
end

-- Dynamical time at Universal moment tee_u.
local function dynamical_from_universal(tee_u)
  return tee_u + ephemeris_correction(tee_u)
end

-- Julian centuries since 2000 at moment tee.
julian_centuries = function(tee)
  return (dynamical_from_universal(tee) - J2000) / 36525
end


-- === Sidereal time ===

--- Mean sidereal time of day from moment `tee` expressed as hour angle.
-- Adapted from Meeus, "Astronomical Algorithms," 2nd edn., 1998, p. 88.
-- @tparam number tee Universal time moment.
-- @treturn number Hour angle in degrees.
function M.sidereal_from_moment(tee)
  local c = (tee - J2000) / 36525
  return basic.poly(c, {280.46061837, 36525 * 360.98564736629,
                         0.000387933, -1/38710000}) % 360
end
local sidereal_from_moment = M.sidereal_from_moment


-- === Astronomical constants ===

local MEAN_TROPICAL_YEAR = 365.242189

local MEAN_SIDEREAL_YEAR = 365.25636

local MEAN_SYNODIC_MONTH = 29.530588861

local SPRING = 0

local SUMMER = 90

local AUTUMN = 180

local WINTER = 270

local MORNING = true

local EVENING = false

local NEW = 0

local FIRST_QUARTER = 90

local FULL = 180

local LAST_QUARTER = 270


-- === Solar calculations ===

-- Obliquity of ecliptic at moment tee.
local function obliquity(tee)
  local c = julian_centuries(tee)
  return angle(23, 26, 21.448)
    + basic.poly(c, {0,
                     angle(0, 0, -46.8150),
                     angle(0, 0, -0.00059),
                     angle(0, 0,  0.001813)})
end

--- Declination of object at latitude `beta`, longitude `lambda` at moment `tee`.
-- @tparam number tee    Universal time moment.
-- @tparam number beta   Latitude of object (degrees).
-- @tparam number lambda Longitude of object (degrees).
-- @treturn number Declination in degrees.
function M.declination(tee, beta, lambda)
  local varepsilon = obliquity(tee)
  return arcsin_degrees(
    sin_degrees(beta) * cos_degrees(varepsilon)
    + cos_degrees(beta) * sin_degrees(varepsilon) * sin_degrees(lambda)
  )
end
local declination = M.declination

--- Right ascension of object at latitude `beta`, longitude `lambda` at moment `tee`.
-- @tparam number tee    Universal time moment.
-- @tparam number beta   Latitude of object (degrees).
-- @tparam number lambda Longitude of object (degrees).
-- @treturn number Right ascension in degrees.
function M.right_ascension(tee, beta, lambda)
  local varepsilon = obliquity(tee)
  return arctan_degrees(
    sin_degrees(lambda) * cos_degrees(varepsilon)
    - tan_degrees(beta) * sin_degrees(varepsilon),
    cos_degrees(lambda)
  )
end
local right_ascension = M.right_ascension

-- Longitudinal nutation at moment tee.
local function nutation(tee)
  local c     = julian_centuries(tee)
  local cap_A = basic.poly(c, {124.90, -1934.134, 0.002063})
  local cap_B = basic.poly(c, {201.11, 72001.5377, 0.00057})
  return -0.004778 * sin_degrees(cap_A)
       + -0.0003667 * sin_degrees(cap_B)
end

-- Aberration at moment tee.
local function aberration(tee)
  local c = julian_centuries(tee)
  return 0.0000974 * cos_degrees(177.63 + 35999.01848 * c)
       - 0.005575
end

--- Longitude of sun (in degrees) at moment `tee`.
-- Adapted from Bretagnon and Simon, "Planetary Programs and Tables," 1986.
-- @tparam number tee Universal time moment.
-- @treturn number Solar longitude in degrees.
function M.solar_longitude(tee)
  local c = julian_centuries(tee)
  local coefficients = {
    403406, 195207, 119433, 112392, 3891, 2819, 1721,
    660, 350, 334, 314, 268, 242, 234, 158, 132, 129, 114,
    99, 93, 86, 78, 72, 68, 64, 46, 38, 37, 32, 29, 28, 27, 27,
    25, 24, 21, 21, 20, 18, 17, 14, 13, 13, 13, 12, 10, 10, 10, 10,
  }
  local multipliers = {
    0.9287892, 35999.1376958, 35999.4089666,
    35998.7287385, 71998.20261, 71998.4403,
    36000.35726, 71997.4812, 32964.4678,
    -19.4410, 445267.1117, 45036.8840, 3.1008,
    22518.4434, -19.9739, 65928.9345,
    9038.0293, 3034.7684, 33718.148, 3034.448,
    -2280.773, 29929.992, 31556.493, 149.588,
    9037.750, 107997.405, -4444.176, 151.771,
    67555.316, 31556.080, -4561.540,
    107996.706, 1221.655, 62894.167,
    31437.369, 14578.298, -31931.757,
    34777.243, 1221.999, 62894.511,
    -4442.039, 107997.909, 119.066, 16859.071,
    -4.578, 26895.292, -39.127, 12297.536, 90073.778,
  }
  local addends = {
    270.54861, 340.19128, 63.91854, 331.26220,
    317.843, 86.631, 240.052, 310.26, 247.23,
    260.87, 297.82, 343.14, 166.79, 81.53,
    3.50, 132.75, 182.95, 162.03, 29.8,
    266.4, 249.2, 157.6, 257.8, 185.1, 69.9,
    8.0, 197.1, 250.4, 65.3, 162.7, 341.5,
    291.6, 98.5, 146.7, 110.0, 5.2, 342.6,
    230.9, 256.1, 45.3, 242.9, 115.2, 151.8,
    285.3, 53.3, 126.6, 205.7, 85.9, 146.1,
  }
  local lambda = 282.7771834
    + 36000.76953744 * c
    + 0.000005729577951308232
      * basic.sigma(
          { coefficients, multipliers, addends },
          function(x, z, y)
            return x * sin_degrees(y + z * c)
          end
        )
  return (lambda + aberration(tee) + nutation(tee)) % 360
end
local solar_longitude = M.solar_longitude

--- Moment UT of first time at or after `tee` when solar longitude is `lambda` degrees.
-- @tparam number lambda Target solar longitude (degrees).
-- @tparam number tee    Start moment (UT).
-- @treturn number Universal time moment.
function M.solar_longitude_after(lambda, tee)
  local rate = MEAN_TROPICAL_YEAR / 360
  local tau  = tee + rate * ((lambda - solar_longitude(tee)) % 360)
  local a    = math.max(tee, tau - 5)
  local b    = tau + 5
  return basic.invert_angular(
    solar_longitude, lambda, basic.interval_closed(a, b)
  )
end
local solar_longitude_after = M.solar_longitude_after

--- Moment UT of `season` in Gregorian year `g_year`.
-- @tparam number season Solar longitude of season (e.g. SPRING=0).
-- @tparam number g_year Gregorian year.
-- @treturn number Universal time moment.
function M.season_in_gregorian(season, g_year)
  return solar_longitude_after(season, gregorian.gregorian_new_year(g_year))
end

--- Precession at moment `tee` using 0,0 as J2000 coordinates.
-- Adapted from Meeus, "Astronomical Algorithms," 2nd edn., 1998, pp. 136-137.
-- @tparam number tee Universal time moment.
-- @treturn number Precession angle in degrees.
function M.precession(tee)
  local c     = julian_centuries(tee)
  local eta   = basic.poly(c, {0, secs(47.0029), secs(-0.03302),
                                   secs(0.000060)}) % 360
  local cap_P = basic.poly(c, {174.876384, secs(-869.8089),
                                   secs(0.03536)}) % 360
  local p     = basic.poly(c, {0, secs(5029.0966), secs(1.11113),
                                   secs(0.000006)}) % 360
  local cap_A = cos_degrees(eta) * sin_degrees(cap_P)
  local cap_B = cos_degrees(cap_P)
  local arg   = arctan_degrees(cap_A, cap_B)
  return (p + cap_P - arg) % 360
end
local precession = M.precession

--- Sidereal solar longitude at moment `tee`.
-- `sidereal_start` is the precession value at the Hindu reference epoch.
-- @tparam number tee            Universal time moment.
-- @tparam number sidereal_start Precession offset (degrees).
-- @treturn number Sidereal solar longitude in degrees.
function M.sidereal_solar_longitude(tee, sidereal_start)
  return (solar_longitude(tee) - precession(tee) + sidereal_start) % 360
end

--- Approximate moment at or before `tee` when solar longitude last exceeded `lambda`.
-- @tparam number lambda Solar longitude (degrees).
-- @tparam number tee    Upper bound moment (UT).
-- @treturn number Universal time moment.
function M.estimate_prior_solar_longitude(lambda, tee)
  local rate    = MEAN_TROPICAL_YEAR / 360
  local tau     = tee - rate * ((solar_longitude(tee) - lambda) % 360)
  local cap_delta = basic.mod3(solar_longitude(tau) - lambda, -180, 180)
  return math.min(tee, tau - rate * cap_delta)
end


-- === Equation of time ===

--- Equation of time (solar noon offset) at moment `tee`.
-- @tparam number tee Moment.
-- @treturn number Fraction of a day.
function M.equation_of_time(tee)
  local c           = julian_centuries(tee)
  local lambda      = basic.poly(c, {280.46645, 36000.76983, 0.0003032})
  local anomaly     = basic.poly(c, {357.52910, 35999.05030,
                                     -0.0001559, -0.00000048})
  local eccentricity = basic.poly(c, {0.016708617, -0.000042037,
                                      -0.0000001236})
  local varepsilon  = obliquity(tee)
  local y           = tan_degrees(varepsilon / 2)^2
  local equation    = (1 / (2 * math.pi))
    * (  y * sin_degrees(2 * lambda)
       - 2 * eccentricity * sin_degrees(anomaly)
       + 4 * eccentricity * y * sin_degrees(anomaly)
         * cos_degrees(2 * lambda)
       - 0.5 * y * y * sin_degrees(4 * lambda)
       - 1.25 * eccentricity * eccentricity * sin_degrees(2 * anomaly))
  return basic.sign(equation) * math.min(math.abs(equation), hr(12))
end
equation_of_time = M.equation_of_time


-- === Sunrise / sunset / dawn / dusk ===

-- Sine of angle between sun position at tee and depression alpha at loc.
local function sine_offset(tee, loc, alpha)
  local phi      = latitude(loc)
  local tee_prime = universal_from_local(tee, loc)
  local delta    = declination(tee_prime, 0, solar_longitude(tee_prime))
  return tan_degrees(phi) * tan_degrees(delta)
       + sin_degrees(alpha) / (cos_degrees(delta) * cos_degrees(phi))
end

-- Local time near tee when sun's depression is alpha at loc; bogus if none.
local function approx_moment_of_depression(tee, loc, alpha, early)
  local try  = sine_offset(tee, loc, alpha)
  local date = basic.fixed_from_moment(tee)
  local alt
  if alpha >= 0 then
    alt = early and date or date + 1
  else
    alt = date + hr(12)
  end
  local value = math.abs(try) > 1 and sine_offset(alt, loc, alpha) or try
  if math.abs(value) <= 1 then
    local offset = basic.mod3(
      arcsin_degrees(value) / 360,
      hr(-12), hr(12)
    )
    return local_from_apparent(
      date + (early and (hr(6) - offset) or (hr(18) + offset)),
      loc
    )
  else
    return basic.BOGUS
  end
end

-- Local time near approx when sun's depression is alpha at loc; bogus if none.
local function moment_of_depression(approx, loc, alpha, early)
  local tee = approx_moment_of_depression(approx, loc, alpha, early)
  if tee == basic.BOGUS then
    return basic.BOGUS
  elseif math.abs(approx - tee) < sec(30) then
    return tee
  else
    return moment_of_depression(tee, loc, alpha, early)
  end
end

--- Standard time of dawn on fixed `date` at `loc` with depression angle `alpha`.
-- Returns BOGUS if there is no dawn.
-- @tparam number date  Fixed date.
-- @tparam table  loc   Location.
-- @tparam number alpha Depression angle (degrees).
-- @treturn number Standard time, or BOGUS.
function M.dawn(date, loc, alpha)
  local result = moment_of_depression(date + hr(6), loc, alpha, MORNING)
  if result == basic.BOGUS then
    return basic.BOGUS
  else
    return standard_from_local(result, loc)
  end
end
local dawn = M.dawn

--- Standard time of dusk on fixed `date` at `loc` with depression angle `alpha`.
-- Returns BOGUS if there is no dusk.
-- @tparam number date  Fixed date.
-- @tparam table  loc   Location.
-- @tparam number alpha Depression angle (degrees).
-- @treturn number Standard time, or BOGUS.
function M.dusk(date, loc, alpha)
  local result = moment_of_depression(date + hr(18), loc, alpha, EVENING)
  if result == basic.BOGUS then
    return basic.BOGUS
  else
    return standard_from_local(result, loc)
  end
end
local dusk = M.dusk

-- Refraction angle at moment tee at loc.
local function refraction(tee, loc) -- luacheck: ignore tee
  local h     = math.max(mt(0), elevation(loc))
  local cap_R = mt(6372000)  -- Radius of Earth.
  local dip   = arccos_degrees(cap_R / (cap_R + h))
  return mins(34) + dip + secs(19) * math.sqrt(h)
end

--- Standard time of sunrise on fixed `date` at `loc`.
-- @tparam number date Fixed date.
-- @tparam table  loc  Location.
-- @treturn number Standard time of sunrise, or BOGUS.
function M.sunrise(date, loc)
  local alpha = refraction(date + hr(6), loc) + mins(16)
  return dawn(date, loc, alpha)
end
local sunrise = M.sunrise

--- Standard time of sunset on fixed `date` at `loc`.
-- @tparam number date Fixed date.
-- @tparam table  loc  Location.
-- @treturn number Standard time of sunset, or BOGUS.
function M.sunset(date, loc)
  local alpha = refraction(date + hr(18), loc) + mins(16)
  return dusk(date, loc, alpha)
end
local sunset = M.sunset

-- Standard time of Jewish dusk (Vilna Gaon) on fixed date at loc.
local function jewish_dusk(date, loc)
  return dusk(date, loc, angle(4, 40, 0))
end

-- Standard time of end of Jewish sabbath (Berthold Cohn) on fixed date at loc.
local function jewish_sabbath_ends(date, loc)
  return dusk(date, loc, angle(7, 5, 0))
end

-- Length of daytime temporal hour on fixed date at loc; bogus if no sunrise/sunset.
local function daytime_temporal_hour(date, loc)
  if sunrise(date, loc) == basic.BOGUS
    or sunset(date, loc) == basic.BOGUS
  then
    return basic.BOGUS
  else
    return (sunset(date, loc) - sunrise(date, loc)) / 12
  end
end

-- Length of nighttime temporal hour on fixed date at loc; bogus if no sunrise/sunset.
local function nighttime_temporal_hour(date, loc)
  if sunrise(date + 1, loc) == basic.BOGUS
    or sunset(date, loc) == basic.BOGUS
  then
    return basic.BOGUS
  else
    return (sunrise(date + 1, loc) - sunset(date, loc)) / 12
  end
end

-- Standard time of temporal moment tee at loc; bogus if temporal hour undefined.
local function standard_from_sundial(tee, loc)
  local date = basic.fixed_from_moment(tee)
  local hour = 24 * basic.time_from_moment(tee)
  local h
  if hour >= 6 and hour <= 18 then
    h = daytime_temporal_hour(date, loc)
  elseif hour < 6 then
    h = nighttime_temporal_hour(date - 1, loc)
  else
    h = nighttime_temporal_hour(date, loc)
  end
  if h == basic.BOGUS then
    return basic.BOGUS
  elseif hour >= 6 and hour <= 18 then
    return sunrise(date, loc) + (hour - 6) * h
  elseif hour < 6 then
    return sunset(date - 1, loc) + (hour + 6) * h
  else
    return sunset(date, loc) + (hour - 18) * h
  end
end

-- Standard time of end of Jewish morning on fixed date at loc.
local function jewish_morning_end(date, loc)
  return standard_from_sundial(date + hr(10), loc)
end

-- Standard time of asr (Hanafi rule) on fixed date at loc; bogus if no asr.
local function asr(date, loc)
  local noon    = midday(date, loc)
  local phi     = latitude(loc)
  local delta   = declination(noon, 0, solar_longitude(noon))
  local altitude = arcsin_degrees(
    cos_degrees(delta) * cos_degrees(phi)
    + sin_degrees(delta) * sin_degrees(phi)
  )
  local h = basic.mod3(
    arctan_degrees(tan_degrees(altitude),
                     1 + 2 * tan_degrees(altitude)),
    -90, 90
  )
  if altitude <= 0 then
    return basic.BOGUS
  else
    return dusk(date, loc, -h)
  end
end

-- Standard time of asr (Shafi'i rule) on fixed date at loc; bogus if no asr.
local function alt_asr(date, loc) -- luacheck: ignore
  local noon    = midday(date, loc)
  local phi     = latitude(loc)
  local delta   = declination(noon, 0, solar_longitude(noon))
  local altitude = arcsin_degrees(
    cos_degrees(delta) * cos_degrees(phi)
    + sin_degrees(delta) * sin_degrees(phi)
  )
  local h = basic.mod3(
    arctan_degrees(tan_degrees(altitude),
                     1 + tan_degrees(altitude)),
    -90, 90
  )
  if altitude <= 0 then
    return basic.BOGUS
  else
    return dusk(date, loc, -h)
  end
end


-- === Italian hours ===

-- Local time of dusk in Padua, Italy on date of moment tee.
local function local_zero_hour(tee)
  local date = basic.fixed_from_moment(tee)
  return local_from_standard(
    dusk(date, PADUA, angle(0, 16, 0)) + mn(30),
    PADUA
  )
end

-- Italian time corresponding to local time tee_l.
local function italian_from_local(tee_l) -- luacheck: ignore
  local date = basic.fixed_from_moment(tee_l)
  local z0   = local_zero_hour(tee_l - 1)
  local z    = local_zero_hour(tee_l)
  if tee_l > z then
    return tee_l + (date + 1 - z)
  else
    return tee_l + (date - z0)
  end
end

-- Local time corresponding to Italian time tee.
local function local_from_italian(tee) -- luacheck: ignore
  local date = basic.fixed_from_moment(tee)
  local z    = local_zero_hour(tee - 1)
  return tee - (date - z)
end


-- === Lunar calculations ===

-- Mean longitude of moon in degrees at moment given in Julian centuries c.
local function mean_lunar_longitude(c)
  return basic.poly(c, {218.3164477, 481267.88123421,
                         -0.0015786, 1/538841, -1/65194000}) % 360
end

-- Elongation of moon in degrees at moment given in Julian centuries c.
local function lunar_elongation(c)
  return basic.poly(c, {297.8501921, 445267.1114034,
                         -0.0018819, 1/545868, -1/113065000}) % 360
end

-- Mean anomaly of sun in degrees at moment given in Julian centuries c.
local function solar_anomaly(c)
  return basic.poly(c, {357.5291092, 35999.0502909,
                         -0.0001536, 1/24490000}) % 360
end

-- Mean anomaly of moon in degrees at moment given in Julian centuries c.
local function lunar_anomaly(c)
  return basic.poly(c, {134.9633964, 477198.8675055,
                         0.0087414, 1/69699, -1/14712000}) % 360
end

-- Moon's argument of latitude in degrees at moment given in Julian centuries c.
local function moon_node(c)
  return basic.poly(c, {93.2720950, 483202.0175233,
                         -0.0036539, -1/3526000, 1/863310000}) % 360
end

-- Angular distance of the lunar node from the equinoctial point on fixed date.
local function lunar_node(date) -- luacheck: ignore
  return basic.mod3(moon_node(julian_centuries(date)), -90, 90)
end

--- Longitude of moon in degrees at moment `tee`.
-- Adapted from "Astronomical Algorithms" by Jean Meeus, 2nd edn., 1998, pp. 338-342.
-- @tparam number tee Moment (Julian day or fixed).
-- @treturn number Lunar longitude in degrees.
function M.lunar_longitude(tee)
  local c          = julian_centuries(tee)
  local cap_L_prime = mean_lunar_longitude(c)
  local cap_D      = lunar_elongation(c)
  local cap_M      = solar_anomaly(c)
  local cap_M_prime = lunar_anomaly(c)
  local cap_F      = moon_node(c)
  local cap_E      = basic.poly(c, {1, -0.002516, -0.0000074})
  local args_lunar_elongation = {
    0,2,2,0,0,0,2,2,2,2,0,1,0,2,0,0,4,0,4,2,2,1,
    1,2,2,4,2,0,2,2,1,2,0,0,2,2,2,4,0,3,2,4,0,2,
    2,2,4,0,4,1,2,0,1,3,4,2,0,1,2,
  }
  local args_solar_anomaly = {
    0,0,0,0,1,0,0,-1,0,-1,1,0,1,0,0,0,0,0,0,1,1,
    0,1,-1,0,0,0,1,0,-1,0,-2,1,2,-2,0,0,-1,0,0,1,
    -1,2,2,1,-1,0,0,-1,0,1,0,1,0,0,-1,2,1,0,
  }
  local args_lunar_anomaly = {
    1,-1,0,2,0,0,-2,-1,1,0,-1,0,1,0,1,1,-1,3,-2,
    -1,0,-1,0,1,2,0,-3,-2,-1,-2,1,0,2,0,-1,1,0,
    -1,2,-1,1,-2,-1,-1,-2,0,1,4,0,-2,0,2,1,-2,-3,
    2,1,-1,3,
  }
  local args_moon_node = {
    0,0,0,0,0,2,0,0,0,0,0,0,0,-2,2,-2,0,0,0,0,0,
    0,0,0,0,0,0,0,2,0,0,0,0,0,0,-2,2,0,2,0,0,0,0,
    0,0,-2,0,0,0,0,-2,-2,0,0,0,0,0,0,0,
  }
  local sine_coeff = {
    6288774,1274027,658314,213618,-185116,-114332,
    58793,57066,53322,45758,-40923,-34720,-30383,
    15327,-12528,10980,10675,10034,8548,-7888,
    -6766,-5163,4987,4036,3994,3861,3665,-2689,
    -2602,2390,-2348,2236,-2120,-2069,2048,-1773,
    -1595,1215,-1110,-892,-810,759,-713,-700,691,
    596,549,537,520,-487,-399,-381,351,-340,330,
    327,-323,299,294,
  }
  local correction = (1/1000000) * basic.sigma(
    {sine_coeff, args_lunar_elongation, args_solar_anomaly,
       args_lunar_anomaly, args_moon_node},
    function(v, w, x, y, z)
      return v * cap_E^math.abs(x)
           * sin_degrees(w*cap_D + x*cap_M + y*cap_M_prime + z*cap_F)
    end
  )
  local venus = (3958/1000000)
    * sin_degrees(119.75 + c * 131.849)
  local jupiter = (318/1000000)
    * sin_degrees(53.09 + c * 479264.29)
  local flat_earth = (1962/1000000)
    * sin_degrees(cap_L_prime - cap_F)
  return (cap_L_prime + correction + venus + jupiter
          + flat_earth + nutation(tee)) % 360
end
local lunar_longitude = M.lunar_longitude

-- Sidereal lunar longitude at moment tee; sidereal_start defined in Modern Hindu calendar.
local function sidereal_lunar_longitude(tee, sidereal_start) -- luacheck: ignore
  return (lunar_longitude(tee) - precession(tee) + sidereal_start) % 360
end

-- Moment of n-th new moon after/before the new moon of Jan 11, 1 CE.
local function nth_new_moon(n)
  local n0     = 24724  -- Months from RD 0 until j2000.
  local k      = n - n0
  local c      = k / 1236.85
  local approx = J2000 + basic.poly(c,
    {5.09766, MEAN_SYNODIC_MONTH * 1236.85,
     0.00015437, -0.000000150, 0.00000000073})
  local cap_E  = basic.poly(c, {1, -0.002516, -0.0000074})
  local solar_anom = basic.poly(c,
    {2.5534, 1236.85 * 29.10535670, -0.0000014, -0.00000011})
  local lunar_anom = basic.poly(c,
    {201.5643, 385.81693528 * 1236.85,
     0.0107582, 0.00001238, -0.000000058})
  local moon_arg = basic.poly(c,
    {160.7108, 390.67050284 * 1236.85,
     -0.0016118, -0.00000227, 0.000000011})
  local cap_omega = basic.poly(c,
    {124.7746, -1.56375588 * 1236.85, 0.0020672, 0.00000215})
  local e_factor    = {0,1,0,0,1,1,2,0,0,1,0,1,1,1,0,0,0,0,0,0,0,0,0,0}
  local solar_coeff = {0,1,0,0,-1,1,2,0,0,1,0,1,1,-1,2,0,3,1,0,1,-1,-1,1,0}
  local lunar_coeff = {1,0,2,0,1,1,0,1,1,2,3,0,0,2,1,2,0,1,2,1,1,1,3,4}
  local moon_coeff  = {0,0,0,2,0,0,0,-2,2,0,0,2,-2,0,0,-2,0,-2,2,2,2,-2,0,0}
  local sine_coeff  = {
    -0.40720, 0.17241, 0.01608, 0.01039,
     0.00739,-0.00514, 0.00208,-0.00111,
    -0.00057, 0.00056,-0.00042, 0.00042,
     0.00038,-0.00024,-0.00007, 0.00004,
     0.00004, 0.00003, 0.00003,-0.00003,
     0.00003,-0.00002,-0.00002, 0.00002,
  }
  local correction = -0.00017 * sin_degrees(cap_omega)
    + basic.sigma(
        { sine_coeff, e_factor, solar_coeff, lunar_coeff, moon_coeff },
        function(v, w, x, y, z)
          return v * cap_E^w
               * sin_degrees(x*solar_anom + y*lunar_anom + z*moon_arg)
        end
      )
  local add_const  = {
    251.88, 251.83, 349.42,  84.66, 141.74, 207.14, 154.84,
     34.52, 207.19, 291.34, 161.72, 239.56, 331.55,
  }
  local add_coeff  = {
    0.016321, 26.651886, 36.412478, 18.206239, 53.303771,
    2.453732,  7.306860, 27.261239,  0.121824,  1.844379,
    24.198154, 25.513099, 3.592518,
  }
  local add_factor = {
    0.000165, 0.000164, 0.000126, 0.000110, 0.000062,
    0.000060, 0.000056, 0.000047, 0.000042, 0.000040,
    0.000037, 0.000035, 0.000023,
  }
  local extra = 0.000325 * sin_degrees(
    basic.poly(c, {299.77, 132.8475848, -0.009173})
  )
  local additional = basic.sigma(
    { add_const, add_coeff, add_factor },
    function(i, j, l) return l * sin_degrees(i + j * k) end
  )
  return universal_from_dynamical(approx + correction + extra + additional)
end

--- Lunar phase as an angle in degrees at moment `tee`.
-- 0=new moon, 90=first quarter, 180=full moon, 270=last quarter.
-- @tparam number tee Moment.
-- @treturn number Lunar phase in degrees [0, 360).
function M.lunar_phase(tee)
  local phi = (lunar_longitude(tee) - solar_longitude(tee)) % 360
  local t0  = nth_new_moon(0)
  local n   = math.floor(0.5 + (tee - t0) / MEAN_SYNODIC_MONTH)
  local phi_prime = 360 * (((tee - nth_new_moon(n)) / MEAN_SYNODIC_MONTH) % 1)
  if math.abs(phi - phi_prime) > 180 then
    return phi_prime
  else
    return phi
  end
end
local lunar_phase = M.lunar_phase

--- Moment UT of last new moon before `tee`.
-- @tparam number tee Moment.
-- @treturn number Moment UT of preceding new moon.
function M.new_moon_before(tee)
  local t0  = nth_new_moon(0)
  local phi = lunar_phase(tee)
  local n   = math.floor(0.5 + (tee - t0) / MEAN_SYNODIC_MONTH
                              - phi / 360)
  return nth_new_moon(
    basic.final(n - 1, function(k) return nth_new_moon(k) < tee end)
  )
end

--- Moment UT of first new moon at or after `tee`.
-- @tparam number tee Moment.
-- @treturn number Moment UT of next new moon.
function M.new_moon_at_or_after(tee)
  local t0  = nth_new_moon(0)
  local phi = lunar_phase(tee)
  local n   = math.floor(0.5 + (tee - t0) / MEAN_SYNODIC_MONTH
                              - phi / 360)
  return nth_new_moon(
    basic.next(n, function(k) return nth_new_moon(k) >= tee end)
  )
end

--- Moment UT of the last time at or before `tee` when the lunar phase was `phi` degrees.
-- @tparam number phi Lunar phase in degrees.
-- @tparam number tee Moment.
-- @treturn number Moment UT.
function M.lunar_phase_at_or_before(phi, tee)
  local tau = tee - MEAN_SYNODIC_MONTH * (1/360)
                  * ((lunar_phase(tee) - phi) % 360)
  local a   = tau - 2
  local b   = math.min(tee, tau + 2)
  return basic.invert_angular(lunar_phase, phi, basic.interval_closed(a, b))
end

--- Moment UT of the next time at or after `tee` when the lunar phase is `phi` degrees.
-- @tparam number phi Lunar phase in degrees.
-- @tparam number tee Moment.
-- @treturn number Moment UT.
function M.lunar_phase_at_or_after(phi, tee)
  local tau = tee + MEAN_SYNODIC_MONTH * (1/360)
                  * ((phi - lunar_phase(tee)) % 360)
  local a   = math.max(tee, tau - 2)
  local b   = tau + 2
  return basic.invert_angular(lunar_phase, phi, basic.interval_closed(a, b))
end

--- Latitude of moon in degrees at moment `tee`.
-- Adapted from "Astronomical Algorithms" by Jean Meeus, 2nd edn., 1998, pp. 338-342.
-- @tparam number tee Moment.
-- @treturn number Lunar latitude in degrees.
function M.lunar_latitude(tee)
  local c           = julian_centuries(tee)
  local cap_L_prime = mean_lunar_longitude(c)
  local cap_D       = lunar_elongation(c)
  local cap_M       = solar_anomaly(c)
  local cap_M_prime = lunar_anomaly(c)
  local cap_F       = moon_node(c)
  local cap_E       = basic.poly(c, {1, -0.002516, -0.0000074})
  local args_lunar_elongation = {
    0,0,0,2,2,2,2,0,2,0,2,2,2,2,2,2,2,0,4,0,0,0,
    1,0,0,0,1,0,4,4,0,4,2,2,2,2,0,2,2,2,2,4,2,2,
    0,2,1,1,0,2,1,2,0,4,4,1,4,1,4,2,
  }
  local args_solar_anomaly = {
    0,0,0,0,0,0,0,0,0,0,-1,0,0,1,-1,-1,-1,1,0,1,
    0,1,0,1,1,1,0,0,0,0,0,0,0,0,-1,0,0,0,0,1,1,
    0,-1,-2,0,1,1,1,1,1,0,-1,1,0,-1,0,0,0,-1,-2,
  }
  local args_lunar_anomaly = {
    0,1,1,0,-1,-1,0,2,1,2,0,-2,1,0,-1,0,-1,-1,-1,
    0,0,-1,0,1,1,0,0,3,0,-1,1,-2,0,2,1,-2,3,2,-3,
    -1,0,0,1,0,1,1,0,0,-2,-1,1,-2,2,-2,-1,1,1,-1,
    0,0,
  }
  local args_moon_node = {
    1,1,-1,-1,1,-1,1,1,-1,-1,-1,-1,1,-1,1,1,-1,-1,
    -1,1,3,1,1,1,-1,-1,-1,1,-1,1,-3,1,-3,-1,-1,1,
    -1,1,-1,1,1,1,1,-1,3,-1,-1,1,-1,-1,1,-1,1,-1,
    -1,-1,-1,-1,-1,1,
  }
  local sine_coeff = {
    5128122,280602,277693,173237,55413,46271,32573,
    17198,9266,8822,8216,4324,4200,-3359,2463,2211,
    2065,-1870,1828,-1794,-1749,-1565,-1491,-1475,
    -1410,-1344,-1335,1107,1021,833,777,671,607,
    596,491,-451,439,422,421,-366,-351,331,315,
    302,-283,-229,223,223,-220,-220,-185,181,
    -177,176,166,-164,132,-119,115,107,
  }
  local beta = (1/1000000) * basic.sigma(
    { sine_coeff, args_lunar_elongation, args_solar_anomaly,
       args_lunar_anomaly, args_moon_node },
    function(v, w, x, y, z)
      return v * cap_E^math.abs(x)
           * sin_degrees(w*cap_D + x*cap_M + y*cap_M_prime + z*cap_F)
    end
  )
  local venus = (175/1000000)
    * (sin_degrees(119.75 + c * 131.849 + cap_F)
       + sin_degrees(119.75 + c * 131.849 - cap_F))
  local flat_earth = (-2235/1000000) * sin_degrees(cap_L_prime)
    + (127/1000000) * sin_degrees(cap_L_prime - cap_M_prime)
    + (-115/1000000) * sin_degrees(cap_L_prime + cap_M_prime)
  local extra = (382/1000000)
    * sin_degrees(313.45 + c * 481266.484)
  return beta + venus + flat_earth + extra
end
local lunar_latitude = M.lunar_latitude

--- Geocentric altitude of moon in degrees at `tee` at `loc`, ignoring parallax and refraction.
-- @tparam number tee Moment.
-- @tparam table  loc Location.
-- @treturn number Altitude in degrees.
function M.lunar_altitude(tee, loc)
  local phi    = latitude(loc)
  local psi    = longitude(loc)
  local lambda = lunar_longitude(tee)
  local beta   = lunar_latitude(tee)
  local alpha  = right_ascension(tee, beta, lambda)
  local delta  = declination(tee, beta, lambda)
  local theta0 = sidereal_from_moment(tee)
  local cap_H  = (theta0 + psi - alpha) % 360  -- local hour angle (note: -(-psi) = +psi... wait,
  --original uses (- psi) meaning negate psi... see below)
  -- Original: (mod (- theta0 (- psi) alpha) 360)
  -- (- psi) negates psi, so: theta0 - (-psi) - alpha = theta0 + psi - alpha
  -- But standard formula uses west-positive longitudes... keeping as original.
  local altitude = arcsin_degrees(
    sin_degrees(phi) * sin_degrees(delta)
    + cos_degrees(phi) * cos_degrees(delta) * cos_degrees(cap_H)
  )
  return basic.mod3(altitude, -180, 180)
end
local lunar_altitude = M.lunar_altitude

-- Distance to moon in meters at moment tee.
local function lunar_distance(tee)
  local c           = julian_centuries(tee)
  local cap_D       = lunar_elongation(c)
  local cap_M       = solar_anomaly(c)
  local cap_M_prime = lunar_anomaly(c)
  local cap_F       = moon_node(c)
  local cap_E       = basic.poly(c, {1, -0.002516, -0.0000074})
  local args_lunar_elongation = {
    0,2,2,0,0,0,2,2,2,2,0,1,0,2,0,0,4,0,4,2,2,1,
    1,2,2,4,2,0,2,2,1,2,0,0,2,2,2,4,0,3,2,4,0,2,
    2,2,4,0,4,1,2,0,1,3,4,2,0,1,2,2,
  }
  local args_solar_anomaly = {
    0,0,0,0,1,0,0,-1,0,-1,1,0,1,0,0,0,0,0,0,1,1,
    0,1,-1,0,0,0,1,0,-1,0,-2,1,2,-2,0,0,-1,0,0,1,
    -1,2,2,1,-1,0,0,-1,0,1,0,1,0,0,-1,2,1,0,0,
  }
  local args_lunar_anomaly = {
    1,-1,0,2,0,0,-2,-1,1,0,-1,0,1,0,1,1,-1,3,-2,
    -1,0,-1,0,1,2,0,-3,-2,-1,-2,1,0,2,0,-1,1,0,
    -1,2,-1,1,-2,-1,-1,-2,0,1,4,0,-2,0,2,1,-2,-3,
    2,1,-1,3,-1,
  }
  local args_moon_node = {
    0,0,0,0,0,2,0,0,0,0,0,0,0,-2,2,-2,0,0,0,0,0,
    0,0,0,0,0,0,0,2,0,0,0,0,0,0,-2,2,0,2,0,0,0,0,
    0,0,-2,0,0,0,0,-2,-2,0,0,0,0,0,0,0,-2,
  }
  local cosine_coeff = {
    -20905355,-3699111,-2955968,-569925,48888,-3149,
    246158,-152138,-170733,-204586,-129620,108743,
    104755,10321,0,79661,-34782,-23210,-21636,24208,
    30824,-8379,-16675,-12831,-10445,-11650,14403,
    -7003,0,10056,6322,-9884,5751,0,-4950,4130,0,
    -3958,0,3258,2616,-1897,-2117,2354,0,0,-1423,
    -1117,-1571,-1739,0,-4421,0,0,0,0,1165,0,0,
    8752,
  }
  local correction = basic.sigma(
    { cosine_coeff, args_lunar_elongation, args_solar_anomaly,
       args_lunar_anomaly, args_moon_node },
    function(v, w, x, y, z)
      return v * cap_E^math.abs(x)
           * cos_degrees(w*cap_D + x*cap_M + y*cap_M_prime + z*cap_F)
    end
  )
  return mt(385000560) + correction
end

--- Parallax of moon in degrees at `tee` at `loc`.
-- @tparam number tee Moment.
-- @tparam table  loc Location.
-- @treturn number Parallax in degrees.
function M.lunar_parallax(tee, loc)
  local geo     = lunar_altitude(tee, loc)
  local cap_delta = lunar_distance(tee)
  local alt     = mt(6378140) / cap_delta
  local arg     = alt * cos_degrees(geo)
  return arcsin_degrees(arg)
end
local lunar_parallax = M.lunar_parallax

-- Topocentric altitude of moon in degrees at tee at loc, ignoring refraction.
local function topocentric_lunar_altitude(tee, loc)
  return lunar_altitude(tee, loc) - lunar_parallax(tee, loc)
end

-- Geocentric apparent lunar diameter in degrees at moment tee.
local function lunar_diameter(tee)
  return (1792367000/9) / lunar_distance(tee)
end

-- Observed altitude of upper limb of moon at tee at loc, including refraction and elevation.
local function observed_lunar_altitude(tee, loc)
  return topocentric_lunar_altitude(tee, loc)
    + refraction(tee, loc)
    + mins(16)
end

--- Standard time of moonset on fixed `date` at `loc`.
-- @tparam number date Fixed date.
-- @tparam table  loc  Location.
-- @treturn number Standard time of moonset, or BOGUS.
function M.moonset(date, loc)
  local tee    = universal_from_standard(date, loc)
  local waxing = lunar_phase(tee) < 180
  local alt    = observed_lunar_altitude(tee, loc)
  local lat    = latitude(loc)
  local offset = alt / (4 * (90 - math.abs(lat)))
  local approx
  if waxing then
    approx = offset > 0 and tee + offset or tee + 1 + offset
  else
    approx = tee - offset + 0.5
  end
  local set = basic.binary_search(
    approx - hr(6), approx + hr(6),
    function(x) return observed_lunar_altitude(x, loc) < 0 end,
    function(l, u) return (u - l) < mn(1) end
  )
  if set < tee + 1 then
    return math.max(standard_from_universal(set, loc), date)
  else
    return basic.BOGUS
  end
end

--- Standard time of moonrise on fixed `date` at `loc`.
-- @tparam number date Fixed date.
-- @tparam table  loc  Location.
-- @treturn number Standard time of moonrise, or BOGUS.
function M.moonrise(date, loc)
  local tee    = universal_from_standard(date, loc)
  local waning = lunar_phase(tee) > 180
  local alt    = observed_lunar_altitude(tee, loc)
  local lat    = latitude(loc)
  local offset = alt / (4 * (90 - math.abs(lat)))
  local approx
  if waning then
    approx = offset > 0 and tee + 1 - offset or tee - offset
  else
    approx = tee + 0.5 + offset
  end
  local rise = basic.binary_search(
    approx - hr(6), approx + hr(6),
    function(x) return observed_lunar_altitude(x, loc) > 0 end,
    function(l, u) return (u - l) < mn(1) end
  )
  if rise < tee + 1 then
    return math.max(standard_from_universal(rise, loc), date)
  else
    return basic.BOGUS
  end
end



-- === Exports ===

--- Solar longitude at autumnal equinox (degrees).
M.AUTUMN                         = AUTUMN
--- Signifies evening (false).
M.EVENING                        = EVENING
--- Lunar phase at first quarter (degrees).
M.FIRST_QUARTER                  = FIRST_QUARTER
--- Lunar phase at full moon (degrees).
M.FULL                           = FULL
--- Lunar phase at last quarter (degrees).
M.LAST_QUARTER                   = LAST_QUARTER
--- Mean length of a sidereal year in days.
M.MEAN_SIDEREAL_YEAR             = MEAN_SIDEREAL_YEAR
--- Mean length of a synodic month in days.
M.MEAN_SYNODIC_MONTH             = MEAN_SYNODIC_MONTH
--- Mean length of a tropical year in days.
M.MEAN_TROPICAL_YEAR             = MEAN_TROPICAL_YEAR
--- Location of Mecca.
M.MECCA                          = MECCA
--- Signifies morning (true).
M.MORNING                        = MORNING
--- Lunar phase at new moon (degrees).
M.NEW                            = NEW
--- Location of Padua, Italy.
M.PADUA                          = PADUA
--- Solar longitude at vernal equinox (degrees).
M.SPRING                         = SPRING
--- Solar longitude at summer solstice (degrees).
M.SUMMER                         = SUMMER
--- Solar longitude at winter solstice (degrees).
M.WINTER                         = WINTER

--- Time zone offset (fraction of day) for longitude `phi` (degrees).
-- @function zone_from_longitude
-- @tparam number phi Longitude in degrees.
-- @treturn number Zone offset (fraction of day).
M.zone_from_longitude = zone_from_longitude

--- Universal time from dynamical time `tee` (subtracts ephemeris correction).
-- @function universal_from_dynamical
-- @tparam number tee Dynamical time moment.
-- @treturn number Universal time moment.
M.universal_from_dynamical = universal_from_dynamical

--- Dynamical time from universal time `tee_u` (adds ephemeris correction).
-- @function dynamical_from_universal
-- @tparam number tee_u Universal time moment.
-- @treturn number Dynamical time moment.
M.dynamical_from_universal = dynamical_from_universal

--- Length of daytime temporal hour on fixed `date` at `loc`.
-- Returns BOGUS if there is no sunrise or sunset.
-- @function daytime_temporal_hour
-- @tparam number date Fixed date.
-- @tparam table  loc  Location.
-- @treturn number Duration in days, or BOGUS.
M.daytime_temporal_hour = daytime_temporal_hour

--- Length of nighttime temporal hour on fixed `date` at `loc`.
-- Returns BOGUS if there is no sunrise or sunset.
-- @function nighttime_temporal_hour
-- @tparam number date Fixed date.
-- @tparam table  loc  Location.
-- @treturn number Duration in days, or BOGUS.
M.nighttime_temporal_hour = nighttime_temporal_hour

--- Standard time of temporal moment `tee` at `loc` (sundial time).
-- Returns BOGUS if temporal hour is undefined.
-- @function standard_from_sundial
-- @tparam number tee Moment.
-- @tparam table  loc Location.
-- @treturn number Standard time moment, or BOGUS.
M.standard_from_sundial = standard_from_sundial

--- Standard time of Jewish dusk on fixed `date` at `loc` (sun 4°40' below horizon).
-- @function jewish_dusk
-- @tparam number date Fixed date.
-- @tparam table  loc  Location.
-- @treturn number Standard time moment.
M.jewish_dusk = jewish_dusk

--- Standard time of end of Jewish Sabbath on fixed `date` at `loc` (Berthold Cohn; sun 7°5' below horizon).
-- @function jewish_sabbath_ends
-- @tparam number date Fixed date.
-- @tparam table  loc  Location.
-- @treturn number Standard time moment.
M.jewish_sabbath_ends = jewish_sabbath_ends

--- Standard time of end of Jewish morning on fixed `date` at `loc` (10th temporal hour).
-- @function jewish_morning_end
-- @tparam number date Fixed date.
-- @tparam table  loc  Location.
-- @treturn number Standard time moment.
M.jewish_morning_end = jewish_morning_end

--- Standard time of Asr prayer (Hanafi rule) on fixed `date` at `loc`.
-- Returns BOGUS if no Asr occurs.
-- @function asr
-- @tparam number date Fixed date.
-- @tparam table  loc  Location.
-- @treturn number Standard time moment, or BOGUS.
M.asr = asr

--- Distance from Earth to Moon in meters at moment `tee`.
-- @function lunar_distance
-- @tparam number tee Moment.
-- @treturn number Distance in meters.
M.lunar_distance = lunar_distance

--- Geocentric apparent lunar diameter in degrees at moment `tee`.
-- @function lunar_diameter
-- @tparam number tee Moment.
-- @treturn number Diameter in degrees.
M.lunar_diameter = lunar_diameter

return M