Determining the Date For a Phase of Earth's Moon

[Eddie W. Shore]

Introduction

The program MOONDATE determines when a desired phase of Earth's moon given month and year. The four phases of the moon available are:

New Moon (0.00)
1st Quarter (0.25)
Full Moon (0.50)
3rd Quarter (0.75)

The following equations are used:

Let y = year, m = month, t = type (see the above)

y' = y + (m -1)/12
k = integer( (y - 2000) * 12.3685) + t/4
t = k/1236.85

J = 2451550.09765 + 29.530588853 * k + (1.337 * 10^-4) * t^2 - (1.5 * 10^-7) *t^3 + (7.3 * 10^-10) * t^4

J is the Julian Date and will need to be converted to the Gregorian Date. With the DATE+ function, the HP Prime makes this easy. Look to Dieter's post on this thread for details, link: http://www.hpmuseum.org/forum/thread-121...ulian+date

Caution: This program is designed to work with all dates after January 1, 2000. If you choose any dates before, you may have to adjust the month

Code:

HP Prime Program MOONDATE

EXPORT MOONDATE()
BEGIN
// EWS 2019-02-19
LOCAL P,Y,K,M,T,l,c,s;
LOCAL J,D,N,R;
l:={"New","First Qtr","Full",
"Third Qtr"};
INPUT(
{Y,
{M,{1,2,3,4,5,6,7,8,9,10,11,12}}
,{c,l}},
"Moon Phase Date",
{"Year: ","Month:","Stage:"});
s:=(c-1)/4;
N:=M;

REPEAT 
R:=Y+(N-1)/12;
K:=IP((R-2000)*12.3685)+s;
T:=K/1236.85;
J:=2451550.09765+29.530588853*K
+1.337ᴇ−4*T^2-1.5ᴇ−7*T^3
+7.3ᴇ−10*T^4;
J:=IP(J)-2451545;
D:=DATEADD(2000.0101,J);
N:=N+1;
UNTIL D≥(Y+M/100);

RETURN D;
END;

Note: You can use an alternate of 1/1/2000 (date code 2000.0101) with corresponding Julian Date 2451545.

Example:

January 2019:
New: 2019.0105 (1/5/2019)
1st Qtr: 2019.0113 (1/13/2019)
Full: 2019.0120 (1/20/2019)
3rd Qtr: 2019.0127 (1/27/2019)

Source:
Meenus, Jean. "Astronomical Algorithms" Willmann-Bell, Inc.: Richmond, VA 1991 ISBN 0-943396-35-2

Blog post: https://edspi31415.blogspot.com/2019/02/...se-of.html