calculate easter based on moon positions

easter.html

<h1>Next Easter is <span id="nextEaster"></span></h1>
<h3>Easter is...</h3>
<p>
...the first Sunday after the <i>ecclesiastical</i> full moon that occurs on or soonest after 21 March.
</p>
<p>
In applying the ecclesiastical rules, Christian churches use 21 March as the starting point in determining the date of Easter, from which they find the next full moon, etc.
</p>
<div>above from <a href="https://en.wikipedia.org/wiki/Easter">easter wikipedia page</a></div>
<h3>A full moon is...</h3>
<p>The full moon is the lunar phase when the Moon appears <i>fully illuminated</i> from Earth's perspective.</p>
<p>By definition, its maximum illumination occurs at the moment waxing stops.</p>
<div>above from <a href="https://en.wikipedia.org/wiki/Full_moon">full moon wikipedia page</a></div>

<h2>My code hits <span id="percent"></span>% of <span id="total"></span> test Easters</h2>

Notes:
<ul>
  <li>test data from <a href="http://tlarsen2.tripod.com/thomaslarsen/easterdates.html">here</a></li>
  <li>moon math from <a href="https://github.com/mourner/suncalc/blob/master/suncalc.js">suncalc</a></li>
  <li><a href="https://www.reddit.com/r/ProgrammerHumor/comments/gkl4fv/the_meaning_of_everything/fqsevm6/">comment that encouraged me</a></li>
  <li>I've found the small error comes from some minor full moon detection errors.</li>
</ul>

<script>
// most of this code is from suncalc (except easter)
const PI = Math.PI;
const sin = Math.sin;
const cos = Math.cos;
const tan = Math.tan;
const asin = Math.asin;
const atan = Math.atan2;
const acos = Math.acos;
const rad = PI / 180;
const dayMs = 1000 * 60 * 60 * 24;
const J1970 = 2440588;
const J2000 = 2451545;
const e = rad * 23.43664; // obliquity of the Earth

const toJulian = date => date.valueOf() / dayMs - 0.5 + J1970;
const toDays = date => toJulian(date) - J2000;
const rightAscension = (l, b) => atan(sin(l) * cos(e) - tan(b) * sin(e), cos(l));
const declination = (l, b) => asin(sin(b) * cos(e) + cos(b) * sin(e) * sin(l));
const solarMeanAnomaly = (d) => rad * (357.5291 + 0.98560028 * d);
const dayMonth = d => d.toString().split(' ').slice(1, 3).join(' ');
const tomorrow = date => new Date(date.valueOf() + dayMs);

function eclipticLongitude(M) {
  const C = rad * (1.9148 * sin(M) + 0.02 * sin(2 * M) + 0.0003 * sin(3 * M));
  const P = rad * 102.9372;
  return M + C + P + PI;
}

function sunCoords(d) {
  const M = solarMeanAnomaly(d);
  const L = eclipticLongitude(M);

  return {
    dec: declination(L, 0),
    ra: rightAscension(L, 0)
  };
}

function moonCoords(d) {
  const L = rad * (218.316 + 13.176396 * d); // ecliptic longitude
  const M = rad * (134.963 + 13.064993 * d); // mean anomaly
  const F = rad * (93.272 + 13.229350 * d);  // mean distance
  const l = L + rad * 6.289 * sin(M);        // longitude
  const b = rad * 5.128 * sin(F);            // latitude
  const dt = 385001 - 20905 * cos(M);        // distance to the moon in km

  return {
    ra: rightAscension(l, b),
    dec: declination(l, b),
    dist: dt
  };
}

// suncalc's moon phase code mildly revised
function illumination(date) {
  const d = toDays(date || new Date());
  const s = sunCoords(d);
  const m = moonCoords(d);

  const SUN_DIST = 149597870; // distance from Earth to Sun in km

  const phi = acos(sin(s.dec) * sin(m.dec) + cos(s.dec) * cos(m.dec) * cos(s.ra - m.ra));
  const inc = atan(SUN_DIST * sin(phi), m.dist - SUN_DIST * cos(phi));

  return (1 + cos(inc)) / 2;
};

function easter(year) {
  let date = new Date(year, 2, 21);
  let prev, next;

  // determine the first day with decreasing illumination where the previous day was max illumination
  do {
    prev = illumination(date);
    next = illumination(date = tomorrow(date));
  } while (next > prev || prev < 0.995);

  // find the next sunday
  while (date.getDay() !== 0) date = tomorrow(date);

  return date;
}

// test data
const easters = `11 April 1700     1 April 1725    29 March 1750    16 April 1775    27 March 1701    21 April 1726    11 April 1751     7 April 1776    16 April 1702    13 April 1727     2 April 1752    30 March 1777     8 April 1703    28 March 1728    22 April 1753    19 April 1778    23 March 1704    17 April 1729    14 April 1754     4 April 1779    12 April 1705     9 April 1730    30 March 1755    26 March 1780     4 April 1706    25 March 1731    18 April 1756    15 April 1781    24 April 1707    13 April 1732    10 April 1757    31 March 1782     8 April 1708     5 April 1733    26 March 1758    20 April 1783    31 March 1709    25 April 1734    15 April 1759    11 April 1784    20 April 1710    10 April 1735     6 April 1760    27 March 1785     5 April 1711     1 April 1736    22 March 1761    16 April 1786    27 March 1712    21 April 1737    11 April 1762     8 April 1787    16 April 1713     6 April 1738     3 April 1763    23 March 1788     1 April 1714    29 March 1739    22 April 1764    12 April 1789    21 April 1715    17 April 1740     7 April 1765     4 April 1790    12 April 1716     2 April 1741    30 March 1766    24 April 1791    28 March 1717    25 March 1742    19 April 1767     8 April 1792    17 April 1718    14 April 1743     3 April 1768    31 March 1793     9 April 1719     5 April 1744    26 March 1769    20 April 1794    31 March 1720    18 April 1745    15 April 1770     5 April 1795    13 April 1721    10 April 1746    31 March 1771    27 March 1796     5 April 1722     2 April 1747    19 April 1772    16 April 1797    28 March 1723    14 April 1748    11 April 1773     8 April 1798    16 April 1724     6 April 1749     3 April 1774    24 March 1799    13 April 1800     3 April 1825    31 March 1850    28 March 1875     5 April 1801    26 March 1826    20 April 1851    16 April 1876    18 April 1802    15 April 1827    11 April 1852     1 April 1877    10 April 1803     6 April 1828    27 March 1853    21 April 1878     1 April 1804    19 April 1829    16 April 1854    13 April 1879    14 April 1805    11 April 1830     8 April 1855    28 March 1880     6 April 1806     3 April 1831    23 March 1856    17 April 1881    29 March 1807    22 April 1832    12 April 1857     9 April 1882    17 April 1808     7 April 1833     4 April 1858    25 March 1883     2 April 1809    30 March 1834    24 April 1859    13 April 1884    22 April 1810    19 April 1835     8 April 1860     5 April 1885    14 April 1811     3 April 1836    31 March 1861    25 April 1886    29 March 1812    26 March 1837    20 April 1862    10 April 1887    18 April 1813    15 April 1838     5 April 1863     1 April 1888    10 April 1814    31 March 1839    27 March 1864    21 April 1889    26 March 1815    19 April 1840    16 April 1865     6 April 1890    14 April 1816    11 April 1841     1 April 1866    29 March 1891     6 April 1817    27 March 1842    21 April 1867    17 April 1892    22 March 1818    16 April 1843    12 April 1868     2 April 1893    11 April 1819     7 April 1844    28 March 1869    25 March 1894     2 April 1820    23 March 1845    17 April 1870    14 April 1895    22 April 1821    12 April 1846     9 April 1871     5 April 1896     7 April 1822     4 April 1847    31 March 1872    18 April 1897    30 March 1823    23 April 1848    13 April 1873    10 April 1898    18 April 1824     8 April 1849     5 April 1874     2 April 1899    15 April 1900    12 April 1925     9 April 1950    30 March 1975     7 April 1901     4 April 1926    25 March 1951    18 April 1976    30 March 1902    17 April 1927    13 April 1952    10 April 1977    12 April 1903     8 April 1928     5 April 1953    26 March 1978     3 April 1904    31 March 1929    18 April 1954    15 April 1979    23 April 1905    20 April 1930    10 April 1955     6 April 1980    15 April 1906     5 April 1931     1 April 1956    19 April 1981    31 March 1907    27 March 1932    21 April 1957    11 April 1982    19 April 1908    16 April 1933     6 April 1958     3 April 1983    11 April 1909     1 April 1934    29 March 1959    22 April 1984    27 March 1910    21 April 1935    17 April 1960     7 April 1985    16 April 1911    12 April 1936     2 April 1961    30 March 1986     7 April 1912    28 March 1937    22 April 1962    19 April 1987    23 March 1913    17 April 1938    14 April 1963     3 April 1988    12 April 1914     9 April 1939    29 March 1964    26 March 1989     4 April 1915    24 March 1940    18 April 1965    15 April 1990    23 April 1916    13 April 1941    10 April 1966    31 March 1991     8 April 1917     5 April 1942    26 March 1967    19 April 1992    31 March 1918    25 April 1943    14 April 1968    11 April 1993    20 April 1919     9 April 1944     6 April 1969     3 April 1994     4 April 1920     1 April 1945    29 March 1970    16 April 1995    27 March 1921    21 April 1946    11 April 1971     7 April 1996    16 April 1922     6 April 1947     2 April 1972    30 March 1997     1 April 1923    28 March 1948    22 April 1973    12 April 1998    20 April 1924    17 April 1949    14 April 1974     4 April 1999    23 April 2000    20 April 2025    10 April 2050     7 April 2075    15 April 2001     5 April 2026     2 April 2051    19 April 2076    31 March 2002    28 March 2027    21 April 2052    11 April 2077    20 April 2003    16 April 2028     6 April 2053     3 April 2078    11 April 2004     1 April 2029    29 March 2054    23 April 2079    27 March 2005    21 April 2030    18 April 2055     7 April 2080    16 April 2006    13 April 2031     2 April 2056    30 March 2081     8 April 2007    28 March 2032    22 April 2057    19 April 2082    23 March 2008    17 April 2033    14 April 2058     4 April 2083    12 April 2009     9 April 2034    30 March 2059    26 March 2084     4 April 2010    25 March 2035    18 April 2060    15 April 2085    24 April 2011    13 April 2036    10 April 2061    31 March 2086     8 April 2012     5 April 2037    26 March 2062    20 April 2087    31 March 2013    25 April 2038    15 April 2063    11 April 2088    20 April 2014    10 April 2039     6 April 2064     3 April 2089     5 April 2015     1 April 2040    29 March 2065    16 April 2090    27 March 2016    21 April 2041    11 April 2066     8 April 2091    16 April 2017     6 April 2042     3 April 2067    30 March 2092     1 April 2018    29 March 2043    22 April 2068    12 April 2093    21 April 2019    17 April 2044    14 April 2069     4 April 2094    12 April 2020     9 April 2045    30 March 2070    24 April 2095     4 April 2021    25 March 2046    19 April 2071    15 April 2096    17 April 2022    14 April 2047    10 April 2072    31 March 2097     9 April 2023     5 April 2048    26 March 2073    20 April 2098    31 March 2024    18 April 2049    15 April 2074    12 April 2099    28 March 2100    22 April 2125    12 April 2150     9 April 2175    17 April 2101    14 April 2126     4 April 2151    31 March 2176     9 April 2102    30 March 2127    23 April 2152    20 April 2177    25 March 2103    18 April 2128    15 April 2153     5 April 2178    13 April 2104    10 April 2129    31 March 2154    28 March 2179     5 April 2105    26 March 2130    20 April 2155    16 April 2180    18 April 2106    15 April 2131    11 April 2156     1 April 2181    10 April 2107     6 April 2132    27 March 2157    21 April 2182     1 April 2108    19 April 2133    16 April 2158    13 April 2183    21 April 2109    11 April 2134     8 April 2159    28 March 2184     6 April 2110     3 April 2135    23 March 2160    17 April 2185    29 March 2111    22 April 2136    12 April 2161     9 April 2186    17 April 2112     7 April 2137     4 April 2162    25 March 2187     2 April 2113    30 March 2138    24 April 2163    13 April 2188    22 April 2114    19 April 2139     8 April 2164     5 April 2189    14 April 2115     3 April 2140    31 March 2165    25 April 2190    29 March 2116    26 March 2141    20 April 2166    10 April 2191    18 April 2117    15 April 2142     5 April 2167     1 April 2192    10 April 2118    31 March 2143    27 March 2168    21 April 2193    26 March 2119    19 April 2144    16 April 2169     6 April 2194    14 April 2120    11 April 2145     1 April 2170    29 March 2195     6 April 2121     3 April 2146    21 April 2171    17 April 2196    29 March 2122    16 April 2147    12 April 2172     9 April 2197    11 April 2123     7 April 2148     4 April 2173    25 March 2198     2 April 2124    30 March 2149    17 April 2174    14 April 2199     6 April 2200    27 March 2225    21 April 2250    18 April 2275    19 April 2201    16 April 2226    13 April 2251     2 April 2276    11 April 2202     8 April 2227    28 March 2252    22 April 2277     3 April 2203    23 March 2228    17 April 2253    14 April 2278    22 April 2204    12 April 2229     9 April 2254    30 March 2279     7 April 2205     4 April 2230    25 March 2255    18 April 2280    30 March 2206    24 April 2231    13 April 2256    10 April 2281    19 April 2207     8 April 2232     5 April 2257    26 March 2282     3 April 2208    31 March 2233    25 April 2258    15 April 2283    26 March 2209    20 April 2234    10 April 2259     6 April 2284    15 April 2210     5 April 2235     1 April 2260    22 March 2285    31 March 2211    27 March 2236    21 April 2261    11 April 2286    19 April 2212    16 April 2237     6 April 2262     3 April 2287    11 April 2213     1 April 2238    29 March 2263    22 April 2288    27 March 2214    21 April 2239    17 April 2264     7 April 2289    16 April 2215    12 April 2240     2 April 2265    30 March 2290     7 April 2216     4 April 2241    25 March 2266    19 April 2291    30 March 2217    17 April 2242    14 April 2267    10 April 2292    12 April 2218     9 April 2243     5 April 2268    26 March 2293     4 April 2219    31 March 2244    18 April 2269    15 April 2294    23 April 2220    13 April 2245    10 April 2270     7 April 2295    15 April 2221     5 April 2246     2 April 2271    19 April 2296    31 March 2222    28 March 2247    21 April 2272    11 April 2297    20 April 2223    16 April 2248     6 April 2273     3 April 2298    11 April 2224     1 April 2249    29 March 2274    16 April 2299`.split(/\s\s+/).map(e => new Date(e));

const missedEasters = easters.filter(d => dayMonth(d) !== dayMonth(easter(d.getFullYear())));

// show the missed easters in console
console.log(missedEasters.map(d => {
    const year = d.getFullYear();
    return `[${year}] ${dayMonth(d)} <=> ${dayMonth(easter(year))}`
  })
  .join('\n')
);
document.getElementById('percent').innerText = Math.round(100 - missedEasters.length/easters.length*100);
document.getElementById('total').innerText = easters.length;

// find this year's easter
let year = new Date().getFullYear();
let nextEaster = easter(year);

// if it has passed, find next year's
if (Date.now() > nextEaster.getTime())
  nextEaster = easter(++year);

document.getElementById('nextEaster').innerText = dayMonth(nextEaster) + ', ' + year;

</script>