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!
- What does cls() function do inside a class method?
- Load Registered Component in Azure ML for Pipeline using Python sdk v2
- Can I add message to the tqdm progressbar?
- How can I preserve the previous value to find the row that is greater than it?
- What's a correct and good way to implement __hash__()?
- Python Harvesters Image Acquisition GigeCam
- Why do I get a recursion error when the depth of the expected recursion should be way less than 999?
- How to (intermittently) skip certain cells when running IPython notebook?
- How to give jupyter cell standard input in python?
- Connect EtherCAT Device with pysoem
- Regex to substitute the next two words after a matching point
- What is the time complexity of heapq.nlargest?
- How to check if a given number is a power of two?
- How to move Jupyter notebook cells up/down using keyboard shortcut?
- cannot override sys.excepthook
- Getting Python's unittest results in a tearDown() method
- Is redis in Python asynchronous?
- How to tell Click to always show option defaults
- flask sqlalchemy with circular imports with db models
- .join with a query
- Python's argparse to show program's version with prog and version string formatting
- classifiers in scikit-learn that handle nan/null
- How does improve detection of line in images by open cv?
- Performance Advantages to Iterators?
- Problem with selecting a specific web element with Playwright in Python
- Thread long running python task, doesn't work inside Docker
- How to add a column with JSON representation of rows in Polars DataFrame?
- Pandas insert empty row at 0th position
- What should HTTP 201 response body be when responding to a POST request with large data?
- Can Python's unittest test in parallel, like nose can?