How do I compute correct planet oppositions?

Question

I'm having a hard time to find a way to compute the exact moment when a planet comes in opposition. By definition, planet is said to be in opposition when, from the point of view of Earth , it is at the opposite of the Sun , so my assumption was that it occurs at the moment when the angle was at its maximum value.

To find this, I wrote the following code with Skyfield:

def get_skf_objects():
    return get_loader()('de421.bsp')

def get_angle(time: Time):
    sun_pos = earth.at(time).observe(get_skf_objects()['sun'])
    aster_pos = earth.at(time).observe(get_skf_objects()[aster.skyfield_name]) # e.g. aster.skyfield_name = 'MARS'
    degrees = sun_pos.separation_from(aster_pos).degrees
    return degrees

get_angle.rough_period = 1.0

times, angles = _find_maxima(start_time, end_time, get_angle, epsilon=1./3600/24)

For some reasons I don't really get, this produces values for the times when a opposition actually happens. But my real problem is that, compared to the real opposition times, the hour I get is always false. For instance:

• Mars opposition of July 27, 2018 is given to be at 2:15 UTC, but in reality, it actually was at 5:12 UTC.
• Mars opposition of October 13, 2020 will happen at 23:59 UTC, but Skyfield tells me it will happen on October 14th at 01:55 UTC.

And so on…

I added print calls on the degrees variable in get_angle and on the angles after the call of _find_maxima() and got this result for October 14th:

[177.00547411 177.00728148 177.00435833 176.99671953 176.98440208
 176.96746426 176.94598429 176.92005866 176.88980009 176.85533528
 176.81680257 176.77434952]
[177.00547411 177.00615453 177.00667867 177.00704644 177.00725779
 177.00731269 177.00721115 177.00695318 177.00653883 177.00596817
 177.00524129 177.00435833]
[177.00725779 177.00727941 177.00729586 177.00730713 177.00731323
 177.00731417 177.00730993 177.00730051 177.00728593 177.00726617
 177.00724125 177.00721115]
[177.00731323 177.00731379 177.00731417 177.00731438 177.00731443
 177.0073143  177.00731399 177.00731352 177.00731288 177.00731207
 177.00731108 177.00730993]
[177.00731438 177.0073144  177.00731442 177.00731443 177.00731443
 177.00731443 177.00731442 177.00731441 177.00731439 177.00731436
 177.00731433 177.0073143 ]
[177.00731443 177.00731443 177.00731443 177.00731443 177.00731443
 177.00731443 177.00731443 177.00731443 177.00731443 177.00731443
 177.00731443 177.00731443]
[177.00731443 177.00731443 177.00731443 177.00731443 177.00731443
 177.00731443 177.00731443 177.00731443 177.00731443 177.00731443
 177.00731443 177.00731443]


[177.00731443]

What is happening?

Answer 1

You are on the right track, but there are a few problems with your approach.

1. Your source for the exact time of opposition (did you mention the source? if so then I haven't yet found it in the text of your question) appears to be an incorrect one. The exact time of the 2020 Mars opposition is 23:25 UT, as you can verify at this Google Books link.

2. Your guess that opposition would be defined using the raw separation angle is close but not exact. According to this article on the Wikipedia, opposition is the moment when the Sun's and the planet's ecliptic longitude are 180° away from each other. Among other things, this definition means that in the old days the moment of opposition could be found by subtracting two angles from a table of planetary longitudes, whereas a definition involving raw angular separation would have required spherical trigonometry between a pair of lat-lon (or RA-dec) coordinates.

3. The definition, as it happens, uses apparent positions, so you need to call .apparent() to transform the coordinates you receive.

Putting those adjustments together, see if the following script gives you the correct answer:

from skyfield.api import load, tau, pi
from skyfield.almanac import find_discrete

planets = load('de421.bsp')
sun = planets['sun']
earth = planets['earth']
mars = planets['mars']

ts = load.timescale(builtin=True)

def longitude_difference(t):
    e = earth.at(t)
    s = e.observe(sun).apparent()
    m = e.observe(mars).apparent()
    _, lon1, _ = s.ecliptic_latlon()
    _, lon2, _ = m.ecliptic_latlon()
    return (lon1.degrees - lon2.degrees) - 180 > 0

longitude_difference.rough_period = 300.0

t, b = find_discrete(ts.utc(2020), ts.utc(2021), longitude_difference)
for ti in t:
    print(t.utc_jpl())