Calculate moon illumination given moon age

Question

trying to compute moon illumination percent from moon day.

To get moon day in pseudo code I do:

LUNAR_MONTH = 29.530588853
KNOWN_NEW_MOON = JulianDay("2000-01-06 06:14")
age = (julian_day_now - KNOWN_NEW_MOON).abs % LUNAR_MONTH

So that shows me right now an age of 10.9 days. How do I tell percent illumination?

• PS I've got the above way to calculate moon age from wikipedia.
• PPS Actual Ruby code to do so is in https://stackoverflow.com/q/74987538/520567

Answer 1

The algorithm below produces the fraction of moon illumination based on the Julian Day (rather than the moon age). It is from Meeus' Astronomical Algorithms Second Edition chapter 48.

$$ \begin{align} D &= 297.8501921 + 445267.1114034T - 0.0018819T^2 + \frac{1}{545868.0}T^3 - \frac{1}{113065000.0}T^4 \\ M &= 357.5291092 + 35999.0502909T - 0.0001536T^2 + \frac{1}{24490000.0}T^3 \\ M' &= 134.9633964 + 477198.8675055T + 0.0087414T^2 + \frac{1}{69699}T^3 - \frac{1}{14712000.0}T^4 \\ i &= 180 - D - 6.289 \sin M' + 2.1 \sin M -1.274 \sin(2D - M') -0.658 \sin 2D -0.214 \sin 2M' -0.11 \sin D \\ k &= \frac{1 + \cos i}{2} \end{align} $$

Where $T$ is the number of Julian centuries since J2000: $T=(jd-2451545)/36525.0$. And $k$ is the fraction (from 0.0 to 1.0) of the moon (approximated as a perfect sphere) which is illuminated as viewed from the Geocenter.

$D$, $M$, $M'$ are the Delaunay arguments. Respectively the Moon's mean elongation, the Sun's mean anomaly, and the Moon's mean anomaly.

Here is an example implementation in JavaScript.

As far as an explanation of what all of the contants, etc mean... Sadly there is no deep dark secret explanation guarded by mythical creatures in funny hats who require you to solve puzzles to enter their chamber. Rather, they're just a least squares fit to observations to find the coefficients $a_i$ of a function of the form: $$ f(x) = a_1f_1(x) + a_2f_2(x) + a_2f_2(x) + ... a_if_i(x) $$

For the Delaunay arguments, $f_i(x)$ are just the powers of $T$.

For $i$ the $f_i(x)$ are trig functions of integer multiples of the Delaunay arguments. This method has been used in many ephemeris algorithms, most notably VSOP87. Again, it might seem like there's a deep scientific method involved in figuring out which ones to use. But as described in "Synthetic representation of the Galilean satellites' orbital motions from L1 ephemerides", it's largely a matter of performing a Fourier transform to find good candidates, then trying almost every integer multiple of each argument and just picking the one which fits the observations best. Which solution to pick is largely a balance between computational complexity and accuracy over a desired time frame.