# Finding Quintuple Planetary Alignments with SkyField

I am working on a python script that will use the SkyField and SciPy libraries to find quintuple planetary conjunctions and their corresponding constellation location. Specifically I am looking for dates 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 13K years or so...

I found this SkyField solution here to find conjunctions.

I was able to modify the above solution to find quintuple conjunctions for the last 15000 years. At least I think I did. Here is my solution:

import scipy.optimize

efile = "de431t.bsp"; # ephemeris to use

# Define planets
earth = eph['earth barycenter']
venus = eph['venus barycenter']
mercury = eph['mercury barycenter']
mars = eph['mars barycenter']
jupiter = eph['jupiter barycenter']
saturn = eph['saturn barycenter']

# Every month from start year
t = ts.utc(-12999, range(12 * 15000))

print("\nCalculation plaetary locations. This may take a while...\n")

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

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

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

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

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

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

print("Looking for conjunctions...\n")

# When was Mercury conjoined with the other planets?  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 = (sl - ml + pi) % tau - pi

# Find where all planets are within a degrees of one another...
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()
lat, lon, distance = e.observe(mercury).ecliptic_latlon()
relative_lon = (vl - ml + pi) % tau - pi
return relative_lon

for i in conjunctions.nonzero():
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())
e = earth.at(ts.tt(jd=jd_conjunction))
print("In constellation:", constellation_at(vr))
print()


This seems to be working and the next step is to plug the dates into XePhem and see what they look like.

Can anyone confirm I am doing this right?

• I assume you don't want general programming advice (because this isn't Code Review) And you say it seems to be working (for example you could check for much more common alignments, eg three planets within 10 degrees.) And there is no reason to think that the underlying model is flawed. So I don't know if there can be much of an answer beyond "looks okay" – James K Dec 15 '19 at 22:06
• @JamesK I guess I was kind of looking for a "looks OK" answer.. Though the given answer is interesting, because I am concerned I might be missing some things... however plugging the dates into Xephem shows that am pretty close... but not exact... because it was code I was able to more or less figure out what's going on and make it work(ish), however I am by no means an astrophysicist or mathematician :/ – Joshua Besneatte Dec 17 '19 at 17:24

I computed all major planet conjunctions in DE431 to answer How to calculate conjunctions of 2 planets and you might be able to port what I did (using CSPICE) to skyfield.

Several of the more interesting conjunctions I found are here:

http://search.astro.barrycarter.info/table.html

including 5 and 6 planet conjunctions:

I note that the 6 planet conjunction is unique in DE431, but I refer to the 5 planet conjunction as "rare", so it's likely there's at least one other one.

I'm too lazy to look through my own results, but this might give you a start.

Wow, DE431 covers a wide range of dates.

The main weakness of your approach is that you only check every month. The inner planets, in particular, move rapidly between constellations. If Mercury is in Aries one the first day of one month, and Venus the next, your search will not find a match, even though they might both have spent most of the month in the constellation — Venus entering right after the month starts, and Mercury leaving only at month’s end.

Using a smaller step still leaves you vulnerable, it just decreases the fraction of events you will wind up missing.

I would, for each planet, ask Skyfield’s almanac.py search function (there's lots of examples in the file) for a list of the time periods when the planet is in the target location. You can do that with a 0% miss rate, I think, by choosing a high enough step size for each planet that it would have no chance to go through a constellation without your noticing — and you could even guarantee that you didn't miss anything by verifying that the planet never skips over a constellation at the step size you choose, and having the program error out if a planet does. Then I would try doing a union of those periods; finding the union of a list of time periods is a known problem with good solutions, last time I had to do it.

Here is how the almanac searcher works, and some examples:

• I think there might be a better solution but it might be difficult to script an implementation that never misses, using the (usually) slow rate at which the planets' orbits drift. An ephemeris could be sparsely polled, say six points within one period, repeated every 100 periods. Each sextuplet would be used to built a simple, low-order $\mathbf{x}(t)$ interpolator coefficient set and some smart and fast numpy could look for approximate co-linear configurations by inverting some matrix. The trend of the coefficients with time could themselves be parameterized. – uhoh Dec 17 '19 at 3:02
• I know it's easy to suggest that without thinking it through, but something like that might exist already, possibly in a different context (e.g. traffic light timing in Manhattan) and might have an existing name and implementation. – uhoh Dec 17 '19 at 3:04
• I was able to use the dates from my OP code to plug into Xephem and then move around in time until the full alignment was as close as possible. But, no, the dates were definitely not exact. I did find what I was looking for, which was a triple conjunction visible from India... There was only ONE instance I found using the huge ephemeris... However, now I am concerned I may have missed some... I will look into the almanac.py function and get back.... but since I need ALL the planets to conjoin... what if I just ran my search against a slower moving planet, say Mercury vs Saturn??? – Joshua Besneatte Dec 17 '19 at 17:21
• OR maybe I could make the search weekly??? I don't care if it's slow.... do you know how I would chage the timescale to weekly??? I tried changing the range but that didn't work... I am guessing load.timescale() needs to be altered? \ – Joshua Besneatte Dec 17 '19 at 17:32
• fyi I've just asked So, what, exactly, is a tropical period? – uhoh Mar 21 at 8:13