Your conjection-finding code looks very good, and I suspect its results are far better than those of the Wikipedia — looking at the version history, it’s not clear where their numbers even came from, and unattributed astronomy calculations can’t easily be double-checked without knowing from which ephemeris and software they derived them.
I attach below a slightly improved version of your solver. Here are the tweaks I recommend:
- Pass
epoch='date' when computing coordinates, since conjunctions and oppositions are traditionally defined according to the celestial sphere of the year in which they happen, not the standard J2000 celestial sphere.
- Before doing any math on the dates, force them into a range of zero to a full circle (360 degrees or 24 hours). Otherwise, you will see the result of the subtraction jump by ±360°/24h whenever one of the right ascensions or longitudes happens to cross 0°/0h and change signs. This will give you “phantom conjunctions” where the planets are not really swapping places but merely switching the sign of the angle they’re returning. (For example: jumping from -69° to 291° is really no motion in the sky at all.)
- Note that both of our scripts should also find when Jupiter and Saturn are across the sky from each other, since the sign of their difference should also flip at that point.
- In case any sources you track down cite the moment of closest approach or the angle between the planets at that moment, I've added it in.
Here’s the output, very close to yours, and again not agreeing well with those old unexplained uncredited Wikipedia numbers:
Great conjunction in ecliptic longitude:
['A.D. 2020-Dec-21 18:20:37.5144 UT']
Great conjunction in right ascension:
['A.D. 2020-Dec-21 13:32:04.9486 UT']
Great conjunction at closest approach:
A.D. 2020-Dec-21 18:21:00.2161 UT - 0.1018 degrees
And the script:
from skyfield.api import load
from skyfield.searchlib import find_discrete, find_minima
planets = load('de421.bsp')
sun = planets['sun']
earth = planets['earth']
jupiter = planets['jupiter barycenter']
saturn = planets['saturn barycenter']
ts = load.timescale(builtin=True)
def longitude_difference(t):
e = earth.at(t)
j = e.observe(jupiter).apparent()
s = e.observe(saturn).apparent()
_, lon1, _ = s.ecliptic_latlon(epoch='date')
_, lon2, _ = j.ecliptic_latlon(epoch='date')
return (lon1.degrees - lon2.degrees) % 360.0 > 180.0
def ra_difference(t):
e = earth.at(t)
j = e.observe(jupiter).apparent()
s = e.observe(saturn).apparent()
jRa, _, _ = j.radec(epoch='date')
sRa, _, _ = s.radec(epoch='date')
return (sRa.hours - jRa.hours) % 24.0 > 12.0
def separation(t):
e = earth.at(t)
j = e.observe(jupiter).apparent()
s = e.observe(saturn).apparent()
return j.separation_from(s).degrees
longitude_difference.step_days = 30.0
ra_difference.step_days = 30.0
separation.step_days = 30.0
print()
print("Great conjunction in ecliptic longitude:")
t, b = find_discrete(ts.utc(2020), ts.utc(2021), longitude_difference)
for ti in t:
print(t.utc_jpl())
print()
print("Great conjunction in right ascension:")
t, b = find_discrete(ts.utc(2020), ts.utc(2021), ra_difference)
for ti in t:
print(t.utc_jpl())
print()
print("Great conjunction at closest approach:")
t, b = find_minima(ts.utc(2020, 6), ts.utc(2021), separation)
for ti, bi in zip(t, b):
print('{} - {:.4f} degrees'.format(ti.utc_jpl(), bi))
