Как я могу использовать Skyfield или PyEphem для определения дат прошлых пятикратных соединений?

Я пытаюсь написать скрипт Python, который будет использовать библиотеки SkyField и SciPy для поиска пятикратных планетарных соединений... В частности, я ищу дату, когда все 5 видимых планет были соединены в созвездии Овна. Это явление должно быть исключительно редким, и мне просто нужно что-то выяснить, произошло ли это и когда это произошло за последние 100 тысяч лет или около того...

Существует этот скрипт для PyEphem и решение SkyField здесь, чтобы найти соединения.

Я смог изменить решение небесного поля, чтобы найти соединения за последние 5000 лет:

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

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

#aries = star['aries']

# Every month from year 2000 to 2050.
t = ts.utc(-2999, range(12 * 5000))

# Where in the sky were Venus and the Sun on those dates?
e = earth.at(t)

#lat, lon, distance = e.observe(aries).ecliptic_latlon()
#al = lon.radians

lat, lon, distance = e.observe(sun).ecliptic_latlon()
sl = lon.radians

lat, lon, distance = e.observe(venus).ecliptic_latlon()
vl = lon.radians

lat, lon, distance = e.observe(mercury).ecliptic_latlon()
ml = lon.radians

lat, lon, distance = e.observe(mars).ecliptic_latlon()
mal = lon.radians

lat, lon, distance = e.observe(jupiter).ecliptic_latlon()
jl = lon.radians

lat, lon, distance = e.observe(saturn).ecliptic_latlon()
sal = lon.radians

# Where was Mercury conjoined with Venus?  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 - ml + pi) % tau - pi
relative_lon2 = (mal - ml + pi) % tau - pi
relative_lon3 = (jl - ml + pi) % tau - pi
relative_lon4 = (sal - ml + pi) % tau - pi

# Find where Venus passed from being ahead of the Sun to being behind.
conjunctions = (relative_lon >= 0)[:-1] & (relative_lon < 0)[1:] & (relative_lon2 >= 0)[:-1] & (relative_lon2 < 0)[1:] & (relative_lon3 >= 0)[:-1] & (relative_lon3 < 0)[1:] & (relative_lon4 >= 0)[:-1] & (relative_lon4 < 0)[1:]

# For each month that included a conjunction, ask SciPy exactly when
# the conjunction occurred.

def f(jd):
    "Compute how far away in longitude Venus and Mercury are."
    t = ts.tt(jd=jd)
    e = earth.at(t)
    lat, lon, distance = e.observe(venus).ecliptic_latlon()
    vl = lon.radians
    lat, lon, distance = e.observe(mercury).ecliptic_latlon()
    ml = lon.radians
    relative_lon = (vl - ml + pi) % tau - pi
    return relative_lon

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)
    print("Found conjunction:", ts.tt(jd=jd_conjunction).utc_jpl())
    print()

Это выводит 9 дат... это кажется правильным для чего-то такого редкого.

Кто-нибудь может подтвердить, что я делаю это правильно?

Как импортировать звезду? Нужно ли мне его создавать или для этого есть эфемериды? Существуют ли эфемериды до 3000 г. до н.э.?