0

Expanding on the solution provided here: Determine coordinates at conjunction times ... I coded the following to give me all the conjunctions of 5 planets (and the sun) for any given year (within the ephemeris, of course), sorted by date.

CORRECTED TEXT:

My question is how can the results be improved?. The results below are encouraging, but ...

  • Normally one can expect that a yearly scan requires 13 monthly 'starting points', e.g. 1st Jan 2019 to 1st Jan 2020 inclusive. However 14 'starting points' including 1st. Feb 2020 are required (see below). Why is this?

  • Using monthly search 'starting points' is an arbitrary choice. It appears to work well with slowly moving celestial objects, however Mercury dances around like a yo-yo and can cause multiple conjunctions within one month. Switching to weekly 'starting points' with t = ts.utc(yy, 0, range(0, 58*7, 7)) does not appear to help. Why is this?

  • Comparing with USNO data looks good. However, to pick one discrepancy in 2020: although "Feb 26 02h Mercury in inferior conjunction" is detected, "Jan 10 15h Mercury in superior conjunction" is not. Why is this?

Note: The scipy.optimize.brentq search can go either way - forwards or backwards - so it is normal to expect that incorrect (previous/next) years have to be filtered out of the results.

Run the following in Python:

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

# find conjunctions of 5 objects within a given year
# INITIALLY RUN: pip install scipy

import sys
import scipy.optimize
from skyfield.api import load, pi, tau

ts = load.timescale()
eph = load('de421.bsp')
sun = eph['sun']
earth   = eph['earth']
mercury = eph['mercury']
venus   = eph['venus']
jupiter = eph['jupiter barycenter']
saturn  = eph['saturn barycenter']
mars    = eph['mars']

objects = [sun, mercury, venus, mars, jupiter, saturn]
object_name = ['sun', 'mercury', 'venus', 'mars', 'jupiter', 'saturn']
conjunctions_per_year = []

def f(jd,j,k):
    "Compute how far away in longitude the two celestial objects are."
    t = ts.tt(jd=jd)
    e = earth.at(t)
    lat, lon, distance = e.observe(objects[j]).ecliptic_latlon()
    sl = lon.radians
    lat, lon, distance = e.observe(objects[k]).ecliptic_latlon()
    vl = lon.radians
    relative_lon = (vl - sl + pi) % tau - pi
    return relative_lon

## MAIN ##
s = input("""
  Enter year as 'YYYY': """)
okay = False
if len(s) == 4:
    if s.isnumeric():
        yy = int(s)
        if 1900 <= yy <= 2050: okay = True
if not okay:
    print("ERROR: Please pick a year between 1900 and 2050")
    sys.exit(0)

# Process monthly starting points spanning the chosen year
t = ts.utc(yy, range(13))
print("Found conjunctions:")    # let's assume we find some

# Where in the sky were the two celestial objects on those dates?
e = earth.at(t)

for j in range(len(objects)):
    lat, lon, distance = e.observe(objects[j]).ecliptic_latlon()
    sl = lon.radians

    for k in range(j+1, len(objects)):
        lat, lon, distance = e.observe(objects[k]).ecliptic_latlon()
        vl = lon.radians
# Where was object A relative to object B?  Compute their difference in
# longitude, wrapping the value into the range [-pi, pi) to avoid
# the discontinuity when one or the other object reaches 360 degrees
# and flips back to 0 degrees.
        relative_lon = (vl - sl + pi) % tau - pi

# Find where object B passed from being ahead of object A to being behind:
        conjunctions = (relative_lon >= 0)[:-1] & (relative_lon < 0)[1:]

# For each month that included a conjunction,
# ask SciPy exactly when the conjunction occurred.
        for i in conjunctions.nonzero()[0]:
            t0 = t[i]
            t1 = t[i + 1]
            #print("Starting search at", t0.utc_jpl())
            jd_conjunction = scipy.optimize.brentq(f, t[i].tt, t[i+1].tt, args=(j,k))
            # append result as tuple to a list
            conjunctions_per_year.append((jd_conjunction, j, k))

conjunctions_per_year.sort()    # sort tuples in-place by date

for jdt, j, k in conjunctions_per_year:
    tt = ts.tt(jd=jdt)
    #if int(tt.utc_strftime("%Y")) != yy: continue  # filter out incorrect years
    print(" {:7}-{:7}: {}".format(object_name[j], object_name[k], tt.utc_jpl()))

The output generated for 2019 with 14 monthly 'starting points' in line 49 (Jan 1 2019 to Feb 1 2020) is:

  Enter year as 'YYYY': 2019
Found conjunctions:
 mercury-jupiter: A.D. 2018-Dec-21 17:37:00.2965 UTC
 sun    -saturn : A.D. 2019-Jan-02 05:49:31.2431 UTC
 mercury-saturn : A.D. 2019-Jan-13 13:31:05.0947 UTC
 venus  -jupiter: A.D. 2019-Jan-22 12:25:50.9449 UTC
 venus  -saturn : A.D. 2019-Feb-18 10:51:54.0668 UTC
 sun    -mercury: A.D. 2019-Mar-15 01:47:38.1785 UTC
 mercury-mars   : A.D. 2019-Jun-18 16:04:25.1024 UTC
 sun    -mercury: A.D. 2019-Jul-21 12:34:05.4668 UTC
 venus  -mars   : A.D. 2019-Aug-24 17:04:32.8511 UTC
 sun    -mars   : A.D. 2019-Sep-02 10:42:14.4417 UTC
 mercury-mars   : A.D. 2019-Sep-03 15:39:54.5854 UTC
 mercury-venus  : A.D. 2019-Sep-13 15:10:34.7771 UTC
 sun    -mercury: A.D. 2019-Nov-11 15:21:41.6804 UTC
 venus  -jupiter: A.D. 2019-Nov-24 13:33:28.2480 UTC
 venus  -saturn : A.D. 2019-Dec-11 10:04:50.4542 UTC
 sun    -jupiter: A.D. 2019-Dec-27 18:25:26.4797 UTC

Furthermore, with the expected 13 monthly 'starting points' (in line 49) the December 2019 conjunctions are excluded:

  Enter year as 'YYYY': 2019
Found conjunctions:
 mercury-jupiter: A.D. 2018-Dec-21 17:37:00.2965 UTC
 sun    -saturn : A.D. 2019-Jan-02 05:49:31.2431 UTC
 mercury-saturn : A.D. 2019-Jan-13 13:31:05.0947 UTC
 venus  -jupiter: A.D. 2019-Jan-22 12:25:50.9449 UTC
 venus  -saturn : A.D. 2019-Feb-18 10:51:54.0668 UTC
 sun    -mercury: A.D. 2019-Mar-15 01:47:38.1785 UTC
 mercury-mars   : A.D. 2019-Jun-18 16:04:25.1024 UTC
 sun    -mercury: A.D. 2019-Jul-21 12:34:05.4668 UTC
 venus  -mars   : A.D. 2019-Aug-24 17:04:32.8511 UTC
 sun    -mars   : A.D. 2019-Sep-02 10:42:14.4417 UTC
 mercury-mars   : A.D. 2019-Sep-03 15:39:54.5854 UTC
 mercury-venus  : A.D. 2019-Sep-13 15:10:34.7771 UTC
 sun    -mercury: A.D. 2019-Nov-11 15:21:41.6804 UTC
 venus  -jupiter: A.D. 2019-Nov-24 13:33:28.2480 UTC

Note: Uncomment the penultimate line to filter out incorrect years ('2018' in the example above).

Progman
  • 16,827
  • 6
  • 33
  • 48
Aendie
  • 321
  • 1
  • 13

1 Answers1

2

To answer my own question following further investigation...

Question #2:

Monthly search 'starting points' are adequate when considering any two celestial objects that have a maximum of one conjunction per month. As Sun-Mercury probably have the most frequent conjunctions ... they have typically 6 or 7 conjunctions per year so a monthly interval is sufficient (and a weekly interval brings no benefit).

Question #1:

Each search 'starting point' seeks a conjunction except the last one, i.e. start #n is compared to start #n+1 (as a search ending point). So with 13 monthly 'starting points' (14 dates in total), the first one considers Jan 1 to Feb 1 and the last one considers Dec 1 to Jan 1 of next year. scipy.optimize.brentq only searches within the limits given to it.

Question #3:

I soon realized that the original code sample was only considering celestial objects passing in one direction only ... if it came from the other direction it would not be considered as a conjunction. Solving this now gives us all conjunctions per year!

However there was catch to this approach: every time two celestial objects are 180° apart in apparent celestial longitude, the relative longitude changes sign from negative to positive. However this is not a conjunction of planets ... instead they are in opposition ... the difference being that the change in relative longitude is a huge step, so they are easily detected. These cases can be filtered out ... on the other hand we can include them as a by-product. Next catch: a planet is considered to be in opposition if the other object is the sun: the alignment is planet-earth-sun. But planetA-earth-planetB is not considered as an opposition although these cases are also caught. We simply filter them out.

Excuse me ... due to numerous changes, it's simpler to include the whole code again:

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

# Find conjunctions of 5 celestial objects within a given year
#    ... and oppositions of the 3 superior planets

# INITIALLY RUN: pip install scipy skyfield
import sys
import scipy.optimize
from skyfield.api import load, pi, tau

ts = load.timescale()
eph = load('de421.bsp')
sun = eph['sun']
earth   = eph['earth']
mercury = eph['mercury']
venus   = eph['venus']
jupiter = eph['jupiter barycenter']
saturn  = eph['saturn barycenter']
mars    = eph['mars']

objects = [sun, mercury, venus, mars, jupiter, saturn]
object_name = ['sun', 'mercury', 'venus', 'mars', 'jupiter', 'saturn']
conjunctions_per_year = []
oppositions_per_year = []

def f(jd,j,k):
    # Compute how far away in longitude the two celestial objects are.
    t = ts.tt(jd=jd)
    e = earth.at(t)
    lat, lon, distance = e.observe(objects[j]).ecliptic_latlon()
    sl = lon.radians
    lat, lon, distance = e.observe(objects[k]).ecliptic_latlon()
    vl = lon.radians
    relative_lon = (vl - sl + pi) % tau - pi
    return relative_lon

## MAIN ##
s = input("""
 Enter year as 'YYYY': """)
okay = False
if len(s) == 4:
    if s.isnumeric():
        yy = int(s)
        if 1900 <= yy <= 2050: okay = True
if not okay:
    print("ERROR: Please pick a year between 1900 and 2050")
    sys.exit(0)

# Process monthly starting points spanning the chosen year
t = ts.utc(yy, range(1,14))     # monthly starting and end points

# Where in the sky were the two celestial objects on those dates?
e = earth.at(t)

for j in range(len(objects)):
    lat, lon, distance = e.observe(objects[j]).ecliptic_latlon()
    sl = lon.radians

    for k in range(j+1, len(objects)):
        opp_ndx = []        # collect indices for planets in opposition
        lat, lon, distance = e.observe(objects[k]).ecliptic_latlon()
        vl = lon.radians
# Where was object A relative to object B?  Compute their difference in
# longitude, wrapping the value into the range [-pi, pi) to avoid
# the discontinuity when one or the other object reaches 360 degrees
# and flips back to 0 degrees.
        relative_lon = (vl - sl + pi) % tau - pi

# Find where object B passed from being ahead of object A to being behind:
        conjunctionsInf = (relative_lon >= 0)[:-1] & (relative_lon < 0)[1:]
# Find where object A passed from being ahead of object B to being behind:
        conjunctionsSup = (relative_lon < 0)[:-1] & (relative_lon >= 0)[1:]
# ignore planets in opposition (only within conjunctionsSup) ...
        for i in conjunctionsSup.nonzero()[0]:
            if relative_lon[i+1] - relative_lon[i] > 5.0: opp_ndx.append(i)
# all conjunctions is the sum of both
        conjunctions = conjunctionsInf + conjunctionsSup

# For each month that included a conjunction,
# ask SciPy exactly when the conjunction occurred.
        for i in conjunctions.nonzero()[0]:
            t0 = t[i]
            t1 = t[i + 1]
            #print("Starting search at", t0.utc_jpl())
            jdt = scipy.optimize.brentq(f, t[i].tt, t[i+1].tt, args=(j,k))
            # A search starting on Jan 1 of the chosen year may well
            # find a conjunction in December of the previous year.
            tt = ts.tt(jd=jdt)
            if i in opp_ndx:    # planets in opposition?
                if j == 0:  # ignore if planet is not in opposition to the sun
                    # if j != 0 the planets are also theoretically "in opposition" in
                    # that their apparent geocentric celestial longitudes differ by 180°
                    oppositions_per_year.append((jdt, j, k))
            else:
                # append result as tuple to a list
                conjunctions_per_year.append((jdt, j, k))

conjunctions_per_year.sort()    # sort tuples in-place by date
oppositions_per_year.sort()     # sort tuples in-place by date

print("Found {} conjunctions:".format(len(conjunctions_per_year)))
for jdt, j, k in conjunctions_per_year:
    tt = ts.tt(jd=jdt)
    print(" {:7}-{:7}: {}".format(object_name[j], object_name[k], tt.utc_jpl()))

print("\nFound {} oppositions:".format(len(oppositions_per_year)))
for jdt, j, k in oppositions_per_year:
    tt = ts.tt(jd=jdt)
    print(" {:7}-{:7}: {}".format(object_name[j], object_name[k], tt.utc_jpl()))

By comparison we now have the following results for 2019 (which agrees well with USNO data):

 Enter year as 'YYYY': 2019
Found 22 conjunctions:
 sun    -saturn : A.D. 2019-Jan-02 05:49:31.2431 UTC
 mercury-saturn : A.D. 2019-Jan-13 13:31:05.0947 UTC
 venus  -jupiter: A.D. 2019-Jan-22 12:25:50.9449 UTC
 sun    -mercury: A.D. 2019-Jan-30 02:51:47.8429 UTC
 venus  -saturn : A.D. 2019-Feb-18 10:51:54.0668 UTC
 sun    -mercury: A.D. 2019-Mar-15 01:47:38.1785 UTC
 sun    -mercury: A.D. 2019-May-21 13:06:51.1057 UTC
 mercury-mars   : A.D. 2019-Jun-18 16:04:25.1024 UTC
 mercury-mars   : A.D. 2019-Jul-08 22:27:22.7258 UTC
 sun    -mercury: A.D. 2019-Jul-21 12:34:05.4668 UTC
 mercury-venus  : A.D. 2019-Jul-25 00:26:14.9056 UTC
 sun    -venus  : A.D. 2019-Aug-14 06:07:35.3800 UTC
 venus  -mars   : A.D. 2019-Aug-24 17:04:32.8511 UTC
 sun    -mars   : A.D. 2019-Sep-02 10:42:14.4417 UTC
 mercury-mars   : A.D. 2019-Sep-03 15:39:54.5854 UTC
 sun    -mercury: A.D. 2019-Sep-04 01:40:13.2892 UTC
 mercury-venus  : A.D. 2019-Sep-13 15:10:34.7771 UTC
 mercury-venus  : A.D. 2019-Oct-30 22:05:08.2419 UTC
 sun    -mercury: A.D. 2019-Nov-11 15:21:41.6804 UTC
 venus  -jupiter: A.D. 2019-Nov-24 13:33:28.2480 UTC
 venus  -saturn : A.D. 2019-Dec-11 10:04:50.4542 UTC
 sun    -jupiter: A.D. 2019-Dec-27 18:25:26.4797 UTC

Found 2 oppositions:
 sun    -jupiter: A.D. 2019-Jun-10 15:12:47.7731 UTC
 sun    -saturn : A.D. 2019-Jul-09 16:51:29.7713 UTC

Note that the original question which spawned this discussion (Determine coordinates at conjunction times) only seeks inferior conjunctions with Venus ... here I am looking for both Inferior and Superior conjunctions of all planets (and oppositions of the three superior planets).

Aendie
  • 321
  • 1
  • 13