Calculating Moon Phase/Age of New Moon - Accuracy not feasible

[ironfelix717]

I am not an astronomer. Nor is this discussion point a matter of Access VBA, as much as it is an astronomy/science discussion.

I set out today to implement a simple way of calculating the moon phase (percent of illumination, 100% = full moon).

Originally, I found a PHP script here:
calculating-moon-phase

Re-wrote it in VBA. Cool. Until i checked the accuracy compared to the original authors calculator (matches) and then to various web calculators out there. That script was off. By quite a bit. Trashed that. And decided to just learn the math behind calculating the moon phases.

I wasn't looking for "NASA Level Accuracy", simply, the ability to get the date when a full moon is present within accuracy of +/- 12 hours.

I then found this resource to learn how to calculate the new moon age. Turns out its actually quite easy.
https://www.subsystems.us moonphase.pdf

So, I followed their example. Built a script. And wah-lah. Doesn't match. Granted, I used a different "Base Line" date (i'm not going back to B.C.). The date I used for a new moon baseline is 1-6-2000 (see script below). The result for moon age (period age) should be the same. Not way off from their example. Perhaps I have an error somewhere. Doesn't look like it. Could be a date conversion or something.

Regardless, the PDF example's result does not match any of the calculators im basing my results on. I trust the web calculators more than this PDF by far.
https://www.almanac.com/astronomy/moon/calendar/MN/Minneapolis/2000-01
http://www.astronomyknowhow.com/month-percentage.php

In short, I find that you cannot actually calculate this with a simple math formula like the PDF and other's who have attempted this claim. There are either other astronomical or time factors that influence this. Its not feasible to get the accuracy that I request. If you believe i'm wrong, I welcome you to prove that

Best regards

Function MoonPhase(Target As Date) As String
Dim BaseLine            As Date
Dim SinceBase           As Long
Dim PeriodDays          As Long
Dim MoonAge             As Double
Dim MoonCycles          As Double
Dim remainder           As Double


BaseLine = "1-6-2000"  'first new moon of year 200, arbitrary baseline date of a new moon.

SinceBase = Abs(DateDiff("d", BaseLine, Target)) + 2  'days since baseline date

'days in a period of full moon
PeriodDays = 29.53058

'# of mooncycles that occurred between the target and baseline
MoonCycles = SinceBase / PeriodDays

'Calculate the remainder of the fraction
    'mod function is trash and doesn't actually give a remainder today?
    'MoonCycles = SinceBase Mod PeriodSeconds
remainder = "." & Split(CStr(MoonCycles), ".")(1)


'calculate age (days) of the remainder:
MoonAge = remainder * 29.53058
Debug.Print MoonAge

End Function

 

[plog]

I think you are failing in your remainder calculation. Your only using the decimal part , but you also need an integer. You need the remaining value of SinceBase after you subtract out the most whole perioddays that you can without SinceBase falling below 0.

For 50 days after 1-6-2000 you need the remainder to be 20.46942 not .46942 which you are doing now:

WholePeriods As Int

WholePeriods = Int(SinceBase/PeriodDays)

remainder = SinceBase - (WholePeriods*PeriodDays)

 

[The_Doc_Man]

So far as I recall, calculating a lunar phase requires a couple of factors. First, the moon's phase is based on an ellipse, so you have to broaden two phases and narrow two phases. Second, the major axis of the ellipse undergoes a 19-year cycle called nutation. Third, the phase depends on the position of the Earth (by a very small amount) because the Earth is not stationary. Therefore you have a period when the Moon is speeding ahead of the Earth and a period where it appears to be slower because you have a circle riding a circle - an epicycle. When you do all of that, you end up with a serious set of a function of a function of a function. Given that the computation depends on elliptical orbits, you can't really use a "standard" trig function like Sine or Cosine for ANY of the components. The best you can do is an approximation.

The most promising methods seem to take a date and compare it to a reference date where the New Moon is known. Compute the number of cycles and then compute the remainder. From that, look up phases in some kind of "remainder range" table that has the ellipses taken into account. So it is a "divide and conquer" method - compute cycles, then compute remainder to find phase. But of course, there might be some interpolation required for some of the phases.

Therefore, I agree with your conclusion that accuracy will be limited at best. However, you should be able to get close. I looked up some computations and I found a LOT more sites that include various factors in the computation.

 

[ironfelix717]

I think you are failing in your remainder calculation. Your only using the decimal part , but you also need an integer. You need the remaining value of SinceBase after you subtract out the most whole perioddays that you can without SinceBase falling below 0.

For 50 days after 1-6-2000 you need the remainder to be 20.46942 not .46942 which you are doing now:

WholePeriods As Int

WholePeriods = Int(SinceBase/PeriodDays)

remainder = SinceBase - (WholePeriods*PeriodDays)

I'm not following your logic.

50 days after 1-6-2000 would be: 50/29.53 = 1.693. .693 is the remainder.

 

[plog]

The remainder in a base 10 system. You want the remainder in a base 29.53 system.

50 days after 1-6-2000 is 29.53 days + 20.47 days Or 1 cycle plus 20.47 days. You then use the 20.47/29.53 to get your percentage.

 

[ironfelix717]

@plog
I'll have to revisit this manual script another time, but yeah I see your logic now.

@The_Doc_Man
I appreciate you looking into this, as well as the lesson on astrology .


I decided to go a different route. I found a python resource HERE that is apparently backed by some credible sources/authors.
Should suffice my needs for incredible accuracy for the purpose of simply knowing what day the moon is full

I modified the example provided in the post to suit my specific needs.
I will call from VBA to grab what I need for CLI.

Here is the .py for anyone who would be interested in using this.

Use Like:
To find the percentage of the moon phase (where 100% = full moon, 50% = half moon, not the same as percent illuminated - see py file)
Must provide 3 args in the order of : month, day, year separated by spaces only

Find phase percentage on Jan 24, 2022:
CMD: MoonFinder.py 1 24 2022

To find a full moon date/time for a given month and year
Must provide 2 args in the order of: month, day separated by spaces only

Find full moon(s) for july 2022 (may return 2 moons depending on the time!)
CMD: MoonFinder.py 7 2022


Thanks to contributors and the writers of the PyEphem library.

 

[The_Doc_Man]

I appreciate you looking into this, as well as the lesson on astrology

Ah, but the truth is, we Aquarians are scientists and thus don't believe in astrology.