Invalid results with script that uses skyfield

I'm exploring the possibilities of the magnificent software Skyfield by Brandon Rhodes. I've made a script to calculate conjunctions in Right Ascension between random objects. I use the following script:

from skyfield import almanac
from skyfield.searchlib import find_maxima, find_minima, find_discrete
from skyfield.api import Star, load
from datetime import datetime, date,timedelta
import pytz

planets = load('de430t.bsp')
earth = planets['earth']

x = [
['Aldebaran',[4, 35, 55.2],[16, 30, 33]],
['Regulus',[10, 8, 22.3],[11, 58, 2]],
['Pollux',[7, 45, 18.9],[28, 1, 34]],
['Antares',[16, 29, 24.4],[-26, 25, 55]],
]

ts = load.timescale(builtin=True)
t = ts.now()
tzn = 'Europe/Amsterdam'

tz = pytz.timezone(tzn)
now = datetime(2020, 1, 1, 0, 0, 0)
t0 = ts.utc(tz.localize(now))
t1 = ts.utc(tz.localize(now) + timedelta(days=+365))

def difference(t):
    e = earth.at(t)
    ra11, dec11, distance = e.observe(object).radec()
    ra12, dec12, distance2 = e.observe(planets[target1]).radec()
    diff = ra11.hours - ra12.hours
    return diff >= 0

difference.rough_period = 1.0

for count in range (len(x)):
    object = Star(ra_hours=(x[count][1][0], x[count][1][1], x[count][1][2]),dec_degrees=(x[count][2][0], x[count][2][1], x[count][2][2]))
    target1 = 'venus'

    t, b = find_discrete(t0, t1, difference)
    if len(t) > 0:
        print (f"{x[count][0]} and {target1}")

        for i, (a, b) in enumerate(zip(t, b)):
            e = earth.at(a)
            ra11, dec11, distance = e.observe(object).radec('date')
            ra12, dec12, distance2 = e.observe(planets[target1]).radec('date')
            print (f"Diff: {ra11.hours - ra12.hours}, ra_{x[count][0]}: {ra11.hours}, ra_{target1}: {ra12.hours}")
            print(f"{a.utc_iso()},{dec11._degrees - dec12._degrees}")
            print ("")

I believe this is calculating the time instance that both objects have the same RA.

Unfortunately I get these results:

Aldebaran and venus
Diff: 4.600705365706519, ra_Aldebaran: 4.6176105612629375, ra_venus: 0.016905195556418524
2020-02-07T21:20:06Z,16.962942031863825

Diff: -0.0014205156698450239, ra_Aldebaran: 4.617748605529588, ra_venus: 4.619169121199433
2020-04-17T20:25:49Z,-10.136618155596008

Diff: -0.000670093218860579, ra_Aldebaran: 4.617892691238783, ra_venus: 4.618562784457644
2020-06-08T07:56:08Z,-4.921187478324768

Diff: -0.0001286749609095139, ra_Aldebaran: 4.618000948409605, ra_venus: 4.618129623370515
2020-07-12T06:44:16Z,-0.962286810883981

Regulus and venus
Diff: 10.140247344361857, ra_Regulus: 10.157152539918275, ra_venus: 0.016905195556418524
2020-02-07T21:20:06Z,12.28436748615726

Diff: 5.852858068422506e-06, ra_Regulus: 10.157702949500333, ra_venus: 10.157697096642265
2020-10-02T23:42:45Z,0.0903429562716429

Pollux and venus
Diff: 7.758742204719277, ra_Pollux: 7.7756474002756955, ra_venus: 0.016905195556418524
2020-02-07T21:20:06Z,28.391501492522927

Diff: 0.001226225278400328, ra_Pollux: 7.776229220287632, ra_venus: 7.775002995009232
2020-09-01T16:39:22Z,8.682000412217121

Antares and venus
Diff: 16.493491164600684, ra_Antares: 16.510396360157102, ra_venus: 0.016905195556418524
2020-02-07T21:20:06Z,-26.059118330110437

Diff: 0.000832014094154232, ra_Antares: 16.51126040187071, ra_venus: 16.510428387776557
2020-12-23T00:34:39Z,-5.652225571050259

The line starting with "Diff" is a line to monitor the validity of the output. Diff stands for the calculated difference in RA. It should be close to zero. The other two values are the Right Ascensions of both objects. They should be closely similar. The second line is the result I want and it is the calculated time and the distance in degrees between the objects. As you can see for some reason for each set of objects I get an invalid result for the time instance: 2020-02-07T21:20:06Z and the difference value for that instance is certainly not close to zero. If I change the object venus to moon it is getting worse because every second result is invalid. I checked the other results against the Skychart / Cartes du Ciel software and those checkout.

I can't figure out what's wrong here. Can someone help me please?

Good question! I should add a new section to https://rhodesmill.org/skyfield/searches.html explaining this common behavior seen when subtracting two longitudes or right ascensions. The key to unraveling the mystery is to watch what happens to the angle difference at one of the moments that is showing up in your output as a phantom conjunction. I’ve attached a script which prints this for the very first event you print, between Venus and Aldebaran:

2020-02-07 Difference: -19.33880215224481 Venus RA: 23.93746881891146
2020-02-08 Difference: 4.5908654129343995 Venus RA: 0.007801253732248779

The angle difference jumps between -19.3 and 4.6, which should immediately strike us as suspicious, since those are simply two different names for exactly the same angle! You can confirm this by adding 24.0 to -19.3 and you will come out with an angle very close to 4.6 (give or take the bit of actual motion that Venus accomplished over one day).

Why would the result jump between two aliases for exactly the same angular difference in the sky?

The answer is in the second fact printed above: the RA of Venus. The discontinuity happens at exactly the moment that Venus happens to cross 0h right ascension! Even though 23.93746881891146 and 0.007801253732248779 are very nearly the same angle, they differ by 24.0 because they straddle the place in the sky where we reset how we name right ascensions.

My script below also displays a plot clarifying the situation:

You can see, in the top plot, that it’s at the exact moment that Venus resets its RA back to zero that the RA difference makes a leap of 24.0 hours from one to another alias for the same four-and-a-half-hour difference in RA.

The solution?

Angular differences need to be constrained to a range like [-12h, +12h) to force the choice of one preferred alias for each possible angular value. In Python, as you can see in the script below, the typical maneuver is:

(ra1.hours - ra2.hours + 12.0) % 24.0 - 12.0

This is shown in the second plot above as the “Improved difference” and not only does it correctly hide the February 7 discontinuity, no longer making it look like an event, but it now correctly identifies an opposition between Venus and Aldebaran out near the end of 2020 (at the right edge of the plot) that had previously appeared as merely the difference exceeding -12.0 but that now shines forth as a crucial moment for the angular difference.

Finally, this script checks for oppositions and filters them out of your search result. You will also find a few possible tweaks to your Python code that you might consider as you continue coding in Python. Here goes:

from skyfield.searchlib import find_maxima, find_minima, find_discrete
from skyfield.api import Star, load
from datetime import datetime, date,timedelta
import pytz

planets = load('de430t.bsp')
earth = planets['earth']

stars = [
    ['Aldebaran', (4, 35, 55.2), (16, 30, 33)],
    ['Regulus', (10, 8, 22.3), (11, 58, 2)],
    ['Pollux', (7, 45, 18.9), (28, 1, 34)],
    ['Antares', (16, 29, 24.4), (-26, 25, 55)],
]

ts = load.timescale(builtin=True)
t0 = ts.utc(2020, 1, 1)
t1 = ts.utc(2021, 1, 1)

# Exploring the first bad result for Venus and Aldebaran.

star = Star(ra_hours=stars[0][1], dec_degrees=stars[0][2])

for utc in (2020, 2, 7), (2020, 2, 8):
    t = ts.utc(*utc)
    ra1, _, _ = planets['earth'].at(t).observe(star).radec()
    ra2, _, _ = planets['earth'].at(t).observe(planets['venus']).radec()
    print(t.utc_strftime('%Y-%m-%d'),
          'Difference:', ra1.hours - ra2.hours,
          'Venus RA:', ra2.hours)

# Plot showing how to protect an angular difference against discontinuity.

import matplotlib.pyplot as plt

t = ts.utc(2020, 1, range(365))
e = planets['earth'].at(t)
star = Star(ra_hours=stars[0][1], dec_degrees=stars[0][2])
ra1, _, _, = e.observe(star).radec()
ra2, _, _, = e.observe(planets['venus']).radec()

fig, (ax, ax2) = plt.subplots(2, 1)
ax.plot(t.J, ra2.hours, label='Venus RA')
ax.plot(t.J, ra1.hours - ra2.hours, label='RA difference')
ax.set(xlabel='Year', ylabel='Hours')
ax.grid()
ax.legend()
ax2.plot(t.J, ra2.hours, label='Venus RA')
ax2.plot(t.J, (ra1.hours - ra2.hours + 12.0) % 24.0 - 12.0,
         label='Improved difference')
ax2.set(xlabel='Year', ylabel='Hours')
ax2.grid()
ax2.legend()
fig.savefig('tmp.png')

# So we need to force the difference into the range [-12 hours .. +12 hours]

def difference(t):
    e = earth.at(t)
    ra11, dec11, distance = e.observe(object).radec()
    ra12, dec12, distance2 = e.observe(planets[target1]).radec()
    diff = (ra11.hours - ra12.hours + 12.0) % 24.0 - 12.0
    return diff >= 0

difference.rough_period = 1.0

for name, ra_hms, dec_dms in stars:
    object = Star(ra_hours=ra_hms, dec_degrees=dec_dms)
    target1 = 'venus'

    t, b = find_discrete(t0, t1, difference)
    if len(t) == 0:
        break

    print (f"{name} and {target1}")

    for a, b in zip(t, b):
        e = earth.at(a)
        ra1, dec1, _ = e.observe(object).radec('date')
        ra2, dec2, _ = e.observe(planets[target1]).radec('date')
        if abs(ra1.hours - ra2.hours) > 6.0:
            continue  # ignore oppositions
        print(f"Diff: {ra1.hours - ra2.hours:.4f}, ra_{name}: {ra1.hours}, ra_{target1}: {ra2.hours}")
        print(f"{a.utc_iso()},{dec1._degrees - dec2._degrees}")
        print()

The events now printed are:

Aldebaran and venus
Diff: -0.0014, ra_Aldebaran: 4.617748605529591, ra_venus: 4.619169121320681
2020-04-17T20:25:49Z,-10.136618155920797

Diff: -0.0007, ra_Aldebaran: 4.617892691238804, ra_venus: 4.618562784294269
2020-06-08T07:56:08Z,-4.921187477025711

Diff: -0.0001, ra_Aldebaran: 4.618000948409604, ra_venus: 4.618129623302464
2020-07-12T06:44:16Z,-0.9622868107603999

Regulus and venus
Diff: 0.0000, ra_Regulus: 10.157702949500333, ra_venus: 10.157697096607553
2020-10-02T23:42:45Z,0.09034295610918264

Pollux and venus
Diff: 0.0012, ra_Pollux: 7.776229220287636, ra_venus: 7.775002995218732
2020-09-01T16:39:22Z,8.68200041253418

Antares and venus
Diff: 0.0008, ra_Antares: 16.511260401870718, ra_venus: 16.51042838802181
2020-12-23T00:34:39Z,-5.652225570418828

Which I believe addresses and corrects your problem!