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?
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.
const T=(jd-2451545)/36525.0;- what is this? Lets assume I don't care so much to understand all other constants but I'd like to calculate this according to the method I use to get julian day. Otherwise I'm not sure whether my julian day is same as yours. – akostadinov Jan 03 '23 at 10:3500:00, not noon. Now I'm within 4 minutes accuracy compared with mooninfo.org for the next full moon. While this accounts for the inaccuracy, there is still something unclear about yourjdcalculation for me.DateTime.new(2000, 1, 1, 12, 0, 0).ajd.to_f => 2451545.0as in your code. But if I get a unix timestamp of it, then use your formula I get a very different result:DateTime.new(2000, 1, 1, 12, 0, 0).to_time.to_i / 86400000 + 2440587.5 => 2440597.5. Do you have an explanation to this? – akostadinov Jan 04 '23 at 10:55angle = num % 360thenangle += 360 if angle.negative?. Which definition of modulo should be used? Ruby has aflooredmodulo and it seems to be what we need here. But wanted to double check as I'm not familiar with the math you're using. – akostadinov Mar 28 '23 at 13:15DateTime.new(2000, 1, 1, 12, 0, 0).to_time.to_i => 946728000, then in shell$ TZ=UTC date -d "@946728000" => 1.01.2000 (сб) 12:00:00 UTC. – akostadinov Mar 28 '23 at 16:4486400.0, then I get2451545.0. Maybe in your code the timestamp is milliseconds and not seconds like in my case. Good that things are cleared up. That was interesting. Sad only in comments. – akostadinov Mar 28 '23 at 16:49