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).
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.brentqonly 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:
By comparison we now have the following results for 2019 (which agrees well with USNO data):
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).