I'm pretty sure that you want the
difference in ecliptic longitude of the Moon and Sun. That's what's used in the definition of the phases of the Moon. That is, New Moon (or Dark Moon) is at $0°$, First Quarter is at $90°$, Full Moon is at $180°$, and Last Quarter is at $270°$. This Moon phase angle is also known as the elongation of the Moon.
In Chat you said:
I write my own programming to run every day/hour/minute to get every degree angles from moon or sun. Then compute the formula. I'm not sure this is the way.
That's extreme brute force! Fortunately, there is an easier way.
Here's a plot of the Moon's elongation for the first synodic month of 2022, produced using JPL Horizons. It uses a time step of six hours. The numbers on the $X$ axis are Julian days.
As you can see, the graph is almost a straight line. If you have the $(x, y)$ coordinates of two points on that curve you can estimate the $x$ value for any $y$ value using linear interpolation.
We just need to rearrange the standard equation of a line in two-point form:
$$x = x_0 + \left(\frac{x_1-x_0}{y_1-y_0}\right)(y-y_0)$$
If $(x_0, y_0)$ and $(x_1, y_1)$ are two points on the curve, then $(x, y)$ lies on the line through those points, so it will be very close to the curve. So we can look up the $y$ value in the ephemeris corresponding to that new $x$ value.
We can then repeat the process, replacing whichever of $(x_0, y_0)$ and $(x_1, y_1)$ is furthest from the curve by our new $(x, y)$. We stop when $y$ is close enough to the $y$ value that we want.
There's a slight complication with using this method because elongation is an angle, so it wraps around at $360°=0°$. A simple way to deal with this is to use angles in the range $-180°$ to $180°$, and shift the elongation angle we're searching for to $0°$. Then we just need to use estimates $(x_0, y_0)$ and $(x_1, y_1)$ that are within two weeks of the correct time.
Here's some Python code which performs this process. It gets ephemeris data from Horizons, so it's quite accurate. However, it uses ecliptic longitudes for an observer at the centre of the Earth, and ignores atmospheric refraction. (Horizons can actually calculate longitudes for an observer on the surface of the Earth, if you give it the observer's longitude, latitude and altitude. It can also include an estimate of the effects of refraction).
""" Fetch Sun & Moon ecliptic longitude data from Horizons
to compute the Moon's elongation (phase) for a given UTC date & time.
Also, find the date & time when the elongation has changed by a given delta.
Written by PM 2Ring 2022.07.07
Run as Python!
"""
import re, requests
from functools import lru_cache
from datetime import timedelta
url = "https://ssd.jpl.nasa.gov/api/horizons_file.api"
api_version = "1.0"
base_cmd = """
MAKE_EPHEM=YES
CENTER=500@399
QUANTITIES=31
CAL_FORMAT=BOTH
ANG_FORMAT=DEG
CSV_FORMAT=YES
OBJ_DATA=NO
"""
@lru_cache(maxsize=20)
def fetch_data(target, day, verbose=False):
ttype = "JD" if isinstance(day, float) else "CAL"
cmd = f"""
!$$SOF
{base_cmd}
COMMAND={target}
TLIST_TYPE={ttype}
TLIST='{day}'
!$$EOF
"""
if verbose:
print(cmd)
req = requests.post(url, data={'format': 'text'}, files={'input': ('cmd', cmd)})
m = re.search(r"API VERSION:\s*(\S*)", req.text)
if m is None:
print("Malformed command file, aborting")
print(cmd)
print(req.text)
return None
version = m.group(1)
if version != api_version:
print(f"Warning: API version is {version}, but this script is for {api_version}")
m = re.search(r"(?s)\\\$\\\$SOE(.*)\\\$\\\$EOE", req.text)
if m is None:
print("NO EPHEMERIS")
print(req.text)
return None
if verbose:
print(req.text)
print(' %' * 30)
data = m.group(1)[1:]
row = data.split(', ')
return row[0], row[1], float(row[4])
Horizons body IDs
See https://ssd.jpl.nasa.gov/api/horizons.api?format=text&OBJ_DATA=YES&MAKE_EPHEM=NO&COMMAND=MB
Sun: 10
Moon: 301
def horizons_elong(day):
""" Fetch Sun & Moon ecliptic longitude from Horizons
and thence compute the Moon's elongation
"""
cal, jd, sun_pos = fetch_data(10, day)
moon_pos = fetch_data(301, day)[2]
elong = (moon_pos - sun_pos) % 360
print(cal, jd)
print("Sun ", sun_pos)
print("Moon ", moon_pos)
print("elong", elong, "\n")
return elong, jd
Mean synodic month length
synodic_mean = 29.5305889
Find x given y such that (x, y) is on the line
through (x0, y0) & (x1, y1)
def interpolate(y, x0, y0, x1, y1):
dy = y1 - y0
if dy == 0:
return x0
return x0 + (y - y0) * (x1 - x0) / dy
@interact
def main(start='2000-Jan-06 18:13:38.114', delta=0, auto_update=False):
# If start is purely numeric, treat it as a JD number
try:
start = float(start)
except ValueError:
start = start.strip()
start_elong, start_jd = horizons_elong(start)
if delta == 0:
return
start_jd = float(start_jd)
delta = float(delta)
print("delta", delta)
goal = (start_elong + delta) % 360
print("goal ", goal, "\n")
# Shift coordinates to put goal at 0
def func(jd):
return (180 + horizons_elong(jd)[0] - goal) % 360 - 180
# Initial estimate
x0 = start_jd + (delta / 360) * synodic_mean
# Step forward an hour
x1 = x0 + 1 / 24
y0, y1 = func(x0), func(x1)
# Refine the estimate
for i in range(9):
print("\nRound", i)
x = interpolate(0, x0, y0, x1, y1)
y = func(x)
print(f"{x=}\n{y=}")
if abs(y) < 5e-7:
break
# Replace the point which is furthest from the X axis
if abs(y1) > abs(y0):
x1, y1 = x, y
else:
x0, y0 = x, y
d = x - start_jd
print(f"\n{d} days")
print(timedelta(days=d))
Dates may be given in most of the time formats that Horizons accepts, except for Modified Julian Day number. However, they must be given in UTC, you cannot specify a different timezone or time scale. Note that Horizons uses the Julian calendar for dates prior to 1582-Oct-15.
If you give a standard Julian day number, do not give the JD prefix, just enter the plain number.
In theory, this program should work for any dates in the range 9999 BC to 9999 AD, but there seems to be a minor Horizons bug affecting dates earlier than JD $-999999$ (BC 7451-Feb-25).
delta can be any numeric value. It can also be a numeric expression, eg 5*360 + 90. If delta is zero, it just prints the Sun & Moon longitudes and Moon elongation for that date.
This program uses the current value of the mean synodic month length, $29.5305889$. The mean length is actually increasing, from ~$29.53056$ in 10,000 BC to ~$29.53061$ in 10,000 AD.
Here's the output of the program for the data given in the question.
2022-Jan-04 05:00:00.000 2459583.708333333
Sun 283.8004537
Moon 303.8864265
elong 20.085972800000036
delta 180.0
goal 200.08597280000004
2022-Jan-18 23:22:01.440 2459598.473627778
Sun 298.8430904
Moon 129.8813202
elong 191.0382298
2022-Jan-19 00:22:01.440 2459598.515294444
Sun 298.8854946
Moon 130.3958492
elong 191.51035459999997
Round 0
2022-Jan-19 18:31:51.408 2459599.272122778
Sun 299.6556721
Moon 139.7942933
elong 200.1386212
x=2459599.272122779
y=0.052648399999952744
Round 1
2022-Jan-19 18:25:12.408 2459599.267504722
Sun 299.6509729
Moon 139.7366397
elong 200.0856668
x=2459599.267504725
y=-0.00030600000002323213
Round 2
2022-Jan-19 18:25:14.714 2459599.267531412
Sun 299.6510001
Moon 139.7369729
elong 200.08597280000004
x=2459599.267531411
y=0.0
15.559198077768087 days
15 days, 13:25:14.713919
And here's a link to a live version running on the SageMathCell server. The program itself is plain Python. It uses some Sage features for getting the input data, but they can easily be replaced with standard Python functions if you want to run the program on the command line. It uses the popular requests library to fetch data from Horizons.
