<p>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 <span class="math-container">$0°$</span>, First Quarter is at <span class="math-container">$90°$</span>, Full Moon is at <span class="math-container">$180°$</span>, and Last Quarter is at <span class="math-container">$270°$</span>. This Moon phase angle is also known as the <a href="https://www.britannica.com/science/elongation-astronomy" rel="nofollow noreferrer">elongation</a> of the Moon.</p>
<p>In Chat you said:</p>
<blockquote>
<p>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.</p>
</blockquote>
<p>That's <em>extreme</em> brute force! Fortunately, there <em>is</em> an easier way.</p>
<p>Here's a plot of the Moon's elongation for the first synodic month of 2022, produced using JPL <a href="https://ssd.jpl.nasa.gov/horizons/" rel="nofollow noreferrer">Horizons</a>. It uses a time step of six hours. The numbers on the <span class="math-container">$X$</span> axis are <a href="https://en.m.wikipedia.org/wiki/Julian_day" rel="nofollow noreferrer">Julian days</a>.</p>
<p><a href="https://i.sstatic.net/MxTsp.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/MxTsp.png" alt="Moon elongation, Jan 2022" /></a></p>
<p>As you can see, the graph is <em>almost</em> a straight line. If you have the <span class="math-container">$(x, y)$</span> coordinates of two points on that curve you can estimate the <span class="math-container">$x$</span> value for any <span class="math-container">$y$</span> value using <a href="https://en.wikipedia.org/wiki/Linear_interpolation" rel="nofollow noreferrer">linear interpolation</a>.</p>
<p>We just need to rearrange the standard equation of a line in two-point form:</p>
<p><span class="math-container">$$x = x_0 + \left(\frac{x_1-x_0}{y_1-y_0}\right)(y-y_0)$$</span></p>
<p>If <span class="math-container">$(x_0, y_0)$</span> and <span class="math-container">$(x_1, y_1)$</span> are two points on the curve, then <span class="math-container">$(x, y)$</span> lies on the line through those points, so it will be very close to the curve. So we can look up the <span class="math-container">$y$</span> value in the ephemeris corresponding to that new <span class="math-container">$x$</span> value.</p>
<p>We can then repeat the process, replacing whichever of <span class="math-container">$(x_0, y_0)$</span> and <span class="math-container">$(x_1, y_1)$</span> is furthest from the curve by our new <span class="math-container">$(x, y)$</span>. We stop when <span class="math-container">$y$</span> is close enough to the <span class="math-container">$y$</span> value that we want.</p>
<p>There's a slight complication with using this method because elongation is an angle, so it wraps around at <span class="math-container">$360°=0°$</span>. A simple way to deal with this is to use angles in the range <span class="math-container">$-180°$</span> to <span class="math-container">$180°$</span>, and shift the elongation angle we're searching for to <span class="math-container">$0°$</span>. Then we just need to use estimates <span class="math-container">$(x_0, y_0)$</span> and <span class="math-container">$(x_1, y_1)$</span> that are within two weeks of the correct time.</p>
<hr />
<p>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).</p>
<pre><code>&quot;&quot;&quot; Fetch Sun &amp; Moon ecliptic longitude data from Horizons
    to compute the Moon's elongation (phase) for a given UTC date &amp; time.
    Also, find the date &amp; time when the elongation has changed by a given delta.

    Written by PM 2Ring 2022.07.07
    Run as Python!
&quot;&quot;&quot;

import re, requests
from functools import lru_cache
from datetime import timedelta

url = &quot;https://ssd.jpl.nasa.gov/api/horizons_file.api&quot;
api_version = &quot;1.0&quot;

base_cmd = &quot;&quot;&quot;
MAKE_EPHEM=YES
CENTER=500@399
QUANTITIES=31
CAL_FORMAT=BOTH
ANG_FORMAT=DEG
CSV_FORMAT=YES
OBJ_DATA=NO
&quot;&quot;&quot;

@lru_cache(maxsize=20)
def fetch_data(target, day, verbose=False):
    ttype = &quot;JD&quot; if isinstance(day, float) else &quot;CAL&quot;
    cmd = f&quot;&quot;&quot;
!<span class="math-container">$$SOF
{base_cmd}
COMMAND={target}
TLIST_TYPE={ttype}
TLIST='{day}'
!$$</span>EOF
&quot;&quot;&quot;

    if verbose:
        print(cmd)
    req = requests.post(url, data={'format': 'text'}, files={'input': ('cmd', cmd)})
    m = re.search(r&quot;API VERSION:\s*(\S*)&quot;, req.text)
    if m is None:
        print(&quot;Malformed command file, aborting&quot;)
        print(cmd)
        print(req.text)
        return None

    version = m.group(1)
    if version != api_version:
        print(f&quot;Warning: API version is {version}, but this script is for {api_version}&quot;)

    m = re.search(r&quot;(?s)\\\$\\\$SOE(.*)\\\$\\\$EOE&quot;, req.text)
    if m is None:
        print(&quot;NO EPHEMERIS&quot;)
        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&amp;OBJ_DATA=YES&amp;MAKE_EPHEM=NO&amp;COMMAND=MB
# Sun: 10
# Moon: 301

def horizons_elong(day):
    &quot;&quot;&quot; Fetch Sun &amp; Moon ecliptic longitude from Horizons
        and thence compute the Moon's elongation
    &quot;&quot;&quot;
    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(&quot;Sun  &quot;, sun_pos)
    print(&quot;Moon &quot;, moon_pos)
    print(&quot;elong&quot;, elong, &quot;\n&quot;)
    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) &amp; (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(&quot;delta&quot;, delta)

    goal = (start_elong + delta) % 360
    print(&quot;goal &quot;, goal, &quot;\n&quot;)

    # 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(&quot;\nRound&quot;, i)

        x = interpolate(0, x0, y0, x1, y1)
        y = func(x)
        print(f&quot;{x=}\n{y=}&quot;)

        if abs(y) &lt; 5e-7:
            break

        # Replace the point which is furthest from the X axis
        if abs(y1) &gt; abs(y0):
            x1, y1 = x, y
        else:
            x0, y0 = x, y

    d = x - start_jd
    print(f&quot;\n{d} days&quot;)
    print(timedelta(days=d))
</code></pre>
<p>Dates may be given in most of the <a href="https://ssd.jpl.nasa.gov/horizons/manual.html#time" rel="nofollow noreferrer">time formats that Horizons accepts</a>, 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.</p>
<p>If you give a standard Julian day number, do <em>not</em> give the <code>JD</code> prefix, just enter the plain number.</p>
<p>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 <span class="math-container">$-999999$</span> (BC 7451-Feb-25).</p>
<p><code>delta</code> can be any numeric value. It can also be a numeric expression, eg <code>5*360 + 90</code>. If <code>delta</code> is zero, it just prints the Sun &amp; Moon longitudes and Moon elongation for that date.</p>
<p>This program uses the current value of the mean synodic month length, <span class="math-container">$29.5305889$</span>. The mean length is actually increasing, from ~<span class="math-container">$29.53056$</span> in 10,000 BC to ~<span class="math-container">$29.53061$</span> in 10,000 AD.</p>
<hr />
<p>Here's the output of the program for the data given in the question.</p>
<pre><code> 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
</code></pre>
<p>And here's a link to a <a href="https://sagecell.sagemath.org/?z=eJyVV21z2kYQ_q5fsVHSIKVElnDS2kzVlti4IS3gGpq2E2eYAx1wqTipulOCyvi_d_dOYDBtJ9XYlm53b9_u2d2z67pwxfVsCaNSwlPoZ5kEPktFrsUM0kwuhC4TDgnTDOZFtoLXWSH-yqRyAB-dwSxb5aXmoJfc7G4o4LSPaYGqvHzJFPdhnhXAYCE-cgm_jC9IH0dzWqx4YDR1UpU1YS5kYjTt8eHTEjcRcU8vaoXZkskFT2Ba7VQnPNUscIzGXwuhNdKQfd2H1o2QC2iFrVYQfo0_RuQGY0ZF15VeZvKR47qu44hVnhUaCt7E3z9LrrRyTODzUs50lqUKapG0KCczNltyyyeXjb81m76NP45TFinE4C61zlX75ESpJPiQp4FkigWL7OMJy8XJss7rZC5SHiDFdfDP5CMvFAWM26MgRP-mmM_JbJUQBf3td37sTrrXr7v9-PfuyLnoDsbdm_hlGH5_en7u_PxLZzDujXvdUXwaORednyZXw5t-Zxy_Go5fO53BD9v1ZfcH52L0drskVcNXbyaXnXEnHgxtZr7fBeyt2FqJv3jcCn0n4XOYE4QmBBJPs2LBdRPTUTUBvZ9misdXLEUUtC1mdJVz8v7NpQtiDkIJqTSTM-6ZPfM0Y9rHw1YcXHTZNbtsxHNy5NGTJ6PhlbPZZuLOuRj2-53BZbyxxu-c8U-90Xgy_v26izSyV5PixgZt3DVIRxd1mLhIPfpR-2qdpCcvhNQeGvANCdGAHmwxEeSZ0h4ebNPURrxpIMRXTDfa0NB8rRt3hOaUK-QIiRWCDK-ByhpNisW_s0pXRmWgOCtmS69wO9c9eNu9GfWGg_ateubdjp75rkFiQFr9rbMrTBsMMnnkrttnKXmCdYGVuWJYT-RGE9gUQYk14Pr_FuA96dCcDV6XhTQWbcLucbkKFkVW5l7k72XS8B7FsAfhh57O3V9ZIdGjNlDU200Y16b-xhROSyykJdLUrMCeRFxqJZs9vXcY0T_m0vtO-bdPbhEsXS94Zj67w-7_y-ZgCKa4uje90VHmPiNN_46r482W3oAvGvAMTsM6KtN597P8Lmq_t4DMPiGD-IHKU4Fbm9DYYtX4gRLvwvdN847e16Xl0erFe1T_eNfMYZolFfQuFdJGnMNndSpqUt9Z1McUydNdx8D28XSvMQ2GT7cV2n9FFkrZhijELxoYbYw1ckwb2fVA0-mpHdQ9w_3MMXU8oehhdqhgh_nvabU1ZRsOw9L-kDRBlXKCxU7d577LRaHpcHURo6JjEYzKyrxr2fMyhlDG28k_3yr34Qs4_Sp09qrSmvf3SC7FDu7OowOeyQjytroPmMYyMs27Ce6tdA9wUtM_JISJPmcSVCWzBFO7yqReQsrlQi-dmjhZkUQMrfPg5Wn48uzsnLZd0ehe13O4Qh_xsPSSafDWTah8qrDMjvFUYHU8xk8E9GKJfMxlFfp4rN4aU1ZFdqag57zIsxTHqodjwUrh24pYXCBoY1xhIqtwW29EiyFsPyzKdegcrOBL8Cqz08d6Q9P4vcbvE9SAw86YZzNtfFkxIT0cUoWOG60wDJ-_YfJ5-BVEZ-3otH16FkTRC6w-M-5jdJKVOpuUOV0JDmbfY-jNweihfORlwdMKZLnihZg1QRcc8yU03UkYvLkkzpQXdmoW1X1EVkNc17NZ2fPk6xnHLvmWpSXvFkVWHO8x70BpbKde3WIMaVKDwC4-0LR9UI61nV2eKdh_TPW-VqNoz88JYdqend1fM83qALSG4tZJrfcsMkY3KW_PYzxGK3FUQq6RRg303qK-PobRUsw1doOsSITEU1J0lcXOYC3gKYS1j3Nz7_PQ6yNAedEZgehBllASey5iiTTVTuEKZbe2e1JogVbwGiGwd3JDRjzG9ylDaNr0nNB2wud-7e2C0DyncfiJFQk2OfSktGBBLMcW4RFqaL0wRCqeihgmIER6s_6Kdlm54XgBt83xwDkauQILEgq6cXvn_tGcvJU3WSkTTLaotRk_0Np-GYdHZbwTrXaOPRyIc3ezju9u5aaKd6O-hiCbKg9byzfwkj__-t4neqZYS3_cC1NsecpmNrg8Q8X4f4XAHkU3irJAqtJ2fJDAb8DWQh2Zinz41n6F_qE5Gw5lHd87Dt1hH8iZ8LdyFmO0olFQH_4egueY101yR1NEufvFsfvfgmakihPf_xv_cQlW&amp;lang=python" rel="nofollow noreferrer">live version</a> 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 <a href="https://requests.readthedocs.io/en/latest/" rel="nofollow noreferrer">requests</a> library to fetch data from Horizons.</p>

CC BY-SA 4.0
2022-09-16T06:10:21.807Z
Thanks looking great. I have  questions regarding to other celestial bodies. I realized planets before Earth (Mercury and Venus), has different maximum delta value depends on date. And planets after Earth (Mars above) doesn't has limit. Are those two types of calculation differ from Moon itself? Would like to learn more if you can  provide website references. Thank you.
Arcana
CC BY-SA 4.0
2022-09-16T08:53:14.413Z
@Arcana Yes, it's a bit different calculating elongation of planets because of [apparent retrograde motion](https://en.wikipedia.org/wiki/Apparent_retrograde_motion).
PM 2Ring
CC BY-SA 4.0
2022-09-17T11:09:14.850Z
Is it very complex to calculate other planets?
Arcana
CC BY-SA 4.0
2022-09-17T18:58:33.033Z
@Arcana No, it's just a bit messy. My technique assumes the elongation is approximately linear. I recommend using a graph to find the stationary points.
PM 2Ring
CC BY-SA 4.0
2022-07-07T03:51:52.090Z
2022-07-07T03:51:52.090Z
PM 2Ring
<p>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).</p>
<p>But I'm not sure how the calculation done behind the scene.</p>
<p><strong>Here's the settings:</strong></p>
<p><strong>Planets:</strong> Sun and Moon (Geocentric Ecliptic)</p>
<p><strong>Synodic Degree:</strong> 180</p>
<p><strong>Begin:</strong> 12:00AM, 4 January 2022 GMT-5</p>
<p><strong>End result:</strong> 1:23PM, 19 January 2022 GMT-5</p>
<p>How to calculate <code>end result datetime</code> with synodic degree between Sun and Moon?</p>

CC BY-SA 4.0
2022-07-02T15:32:36.597Z
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](https://en.wikipedia.org/wiki/VSOP_model), or the [SPICE](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/info/intrdctn.html) Toolkit from JPL / NASA.
PM 2Ring
CC BY-SA 4.0
2022-07-03T16:15:56.860Z
@PM2Ring Can it be done using https://www.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.
Arcana
CC BY-SA 4.0
2022-07-04T02:11:43.347Z
Yes I had the stuffs I need practically. But I don't know how to compute the equation theoretically.
Arcana
CC BY-SA 4.0
2022-07-04T03:16:07.633Z
I show how to compute the ecliptic longitude (and right ascension) of the Sun to reasonable accuracy in [this answer](https://astronomy.stackexchange.com/a/49546/16685). Computing the Moon's position is *much* harder. As I said there, [lunar theory](https://en.wikipedia.org/wiki/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.
PM 2Ring
CC BY-SA 4.0
2022-07-04T08:17:36.907Z
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?
Arcana
CC BY-SA 4.0
2022-07-04T14:37:26.690Z
Let us [continue this discussion in chat](https://chat.stackexchange.com/rooms/137517/discussion-between-pm-2ring-and-arcana).
PM 2Ring
<observational-astronomy><orbit><orbital-elements><synodic-period>
Time taken calculation to reach synodic degrees between planets
Arcana
# Time taken calculation to reach synodic degrees between planets

(Asked by: username_0 on 2022-07-01)

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?

**Tags:** observational-astronomy, orbit, orbital-elements, synodic-period
> **Comments:**
> **username_1**: 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](https://en.wikipedia.org/wiki/VSOP_model), or the [SPICE](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/info/intrdctn.html) Toolkit from JPL / NASA.  
> **username_0**: @PM2Ring Can it be done using https://www.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.  
> **username_0**: Yes I had the stuffs I need practically. But I don't know how to compute the equation theoretically.  
> **username_1**: I show how to compute the ecliptic longitude (and right ascension) of the Sun to reasonable accuracy in [this answer](https://astronomy.stackexchange.com/a/49546/16685). Computing the Moon's position is *much* harder. As I said there, [lunar theory](https://en.wikipedia.org/wiki/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.  
> **username_0**: 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?  
> **username_1**: Let us [continue this discussion in chat](https://chat.stackexchange.com/rooms/137517/discussion-between-pm-2ring-and-arcana).  

---

## 1 answers

### Answer by username_1 on 2022-07-07. Score: 2 | Accepted answer

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.

> **Comments:**
> **username_0**: Thanks looking great. I have  questions regarding to other celestial bodies. I realized planets before Earth (Mercury and Venus), has different maximum delta value depends on date. And planets after Earth (Mars above) doesn't has limit. Are those two types of calculation differ from Moon itself? Would like to learn more if you can  provide website references. Thank you.  
> **username_1**: @Arcana Yes, it's a bit different calculating elongation of planets because of [apparent retrograde motion](https://en.wikipedia.org/wiki/Apparent_retrograde_motion).  
> **username_0**: Is it very complex to calculate other planets?  
> **username_1**: @Arcana No, it's just a bit messy. My technique assumes the elongation is approximately linear. I recommend using a graph to find the stationary points.  

---

