0

I am attempting to write a python script that will use the SkyField and SciPy libraries to find quintuple planetary conjunctions... Specifically I am looking for the date when the 5 visible planets were all in conjunction within the constellation of Aries. This occurrence should be exceptionally rare and I just need something to find if and when it happened in the last 100K years or so...

There is this script for PyEphem and there is the SkyField solution here to find conjunctions.

I was able to modify the skyfield solution to find conjunctions for the last 5000 years:

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()

This outputs 9 dates... this seems about right for something so rare.

Can anyone confirm I am doing this right?

How do I import a star? Do I have to create it or is there an ephemeris for that? Is there an ephemeris for before 3000 BC?

Joshua Besneatte
  • 164
  • 7
  • 1
    You might want to ask on https://astronomy.stackexchange.com/ about what general techniques might make it possible to project planetary orbits that far back and at what confidence levels. If the uncertainty at 100k years is great enough, it might not be possible to know if the planets were in conjunction or not. Only once you know what math technique will predict planet positions that far back will you be in a position to select an astronomy library that uses the correct techniques. – Brandon Rhodes Dec 14 '19 at 23:06
  • perfect thanks! I will move my question and current research over there. thanks! – Joshua Besneatte Dec 15 '19 at 21:13

0 Answers0