1
$\begingroup$

I have a software tool which calculate the synodic period with degree inputs with two planets. Synodic Period is the temporal interval that it takes for an object to reappear at the same point in relation to two other objects (linear nodes).

But I'm not sure how the calculation done behind the scene.

Here's the settings:

Planets: Sun and Moon (Geocentric Ecliptic)

Synodic Degree: 180

Begin: 12:00AM, 4 January 2022 GMT-5

End result: 1:23PM, 19 January 2022 GMT-5

How to calculate end result datetime with synodic degree between Sun and Moon?

$\endgroup$
6
  • $\begingroup$ It's not easy to accurately calculate the positions of Solar system bodies. We can't guess which method your software tool uses. Maybe it uses VSOP, or the SPICE Toolkit from JPL / NASA. $\endgroup$
    – PM 2Ring
    Jul 2 at 15:32
  • $\begingroup$ @PM2Ring Can it be done using astro.com/swisseph/swephprg.htm Swisspeh ? Accuracy is not my concern. If you know overall theory how to calculate, that's would be great. I would like to know more. $\endgroup$
    – Arcana
    Jul 3 at 16:15
  • $\begingroup$ Yes I had the stuffs I need practically. But I don't know how to compute the equation theoretically. $\endgroup$
    – Arcana
    Jul 4 at 2:11
  • $\begingroup$ I show how to compute the ecliptic longitude (and right ascension) of the Sun to reasonable accuracy in this answer. Computing the Moon's position is much harder. As I said there, lunar theory is complicated! To do what you want, you can make a rough estimate using the mean synodic month length (~29.53059 days), and then use software to do a more accurate calculation. $\endgroup$
    – PM 2Ring
    Jul 4 at 3:16
  • $\begingroup$ Let's assume the accuracy for Sun and Moon are correct, what's the formula to calculate the end time, or theoretically find the end time with brute force on ephemeris data? $\endgroup$
    – Arcana
    Jul 4 at 8:17

1 Answer 1

1
$\begingroup$

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.

Moon elongation, Jan 2022

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.

$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .