find_events on highly elliptical orbit satellite inconsistent

[ecdaero]

Hi! Thanks for the work you've done on Skyfield. I'm using it to find rise/set times of objects, and recently came across this corner case(?) where the find_events method is giving incomplete results.

The track of the object that I noticed this on:

I wrote a small program to demo this behavior:

#!/usr/bin/env python3

import datetime as dt
from datetime import datetime

from skyfield.api import EarthSatellite, Time, load, load_file, wgs84
from skyfield.toposlib import GeographicPosition
from skyfield.vectorlib import VectorSum


def main():
    """Attempt to find rise/set events for an example pass."""
    planets = load_file("de421.bsp")

    ts = load.timescale(builtin=True)

    site_latitude = 33
    site_longitude = -118
    site_altitude = 30
    min_elevation_angle = 10

    site = wgs84.latlon(site_latitude, site_longitude, elevation_m=site_altitude)

    # This object has a highly elliptical orbit:
    sat_target = EarthSatellite(
        "1 17328U 87008A   24305.23833080  .00000500  00000-0  22722-2 0  9997",
        "2 17328  63.0833   7.8955 7415190 253.9972  19.6892  2.00632911273252",
        "MOLNIYA 3-31",
        ts,
    )

    earth = planets["earth"]
    groundstation = earth + site
    satellite = earth + sat_target

    # This is an example pass where we noticed we weren't getting the right rise/sets
    obs_start_time = datetime(2024, 10, 31, 5, 28, 6, 0, tzinfo=dt.timezone.utc)
    obs_end_time = datetime(2024, 10, 31, 5, 40, 6, 0, tzinfo=dt.timezone.utc)
    t0 = ts.from_datetime(obs_start_time)
    t1 = ts.from_datetime(obs_end_time)

    event_times, events = find_events(site, groundstation, sat_target, satellite, t0, t1, min_elevation_angle)

    print("This should have found a 'rise':")
    print(events, event_times)
    print("...but it has not found a rise event.\n")

    # However, this *does* work:
    obs_start_time = datetime(2024, 10, 31, 2, 28, 6, 0, tzinfo=dt.timezone.utc)
    obs_end_time = datetime(2024, 10, 31, 7, 40, 6, 0, tzinfo=dt.timezone.utc)
    t0 = ts.from_datetime(obs_start_time)
    t1 = ts.from_datetime(obs_end_time)

    event_times, events = find_events(site, groundstation, sat_target, satellite, t0, t1, min_elevation_angle)

    print("Expand the boundary, and it finds events:")
    print(events, event_times)
    print("...including a rise and culmination.\n")

    obs_long_end_time = datetime(2024, 10, 31, 16, 40, 6, 0, tzinfo=dt.timezone.utc)
    t2 = ts.from_datetime(obs_long_end_time)

    event_times, events = find_events(site, groundstation, sat_target, satellite, t0, t2, min_elevation_angle)

    print("Expand it further, and it should have found a rise, culmination, and set:")
    print(events, event_times)
    print("...but it found a rise, two culminates, and no set.\n")

    # note: adjust the range fidelity to see more detail on the track
    track_time = ts.utc(2024, 10, 31, 5, range(20, 710, 4))
    print_track(groundstation, satellite, track_time)


def find_events(
    site: GeographicPosition,
    site_vector: VectorSum,
    satellite: EarthSatellite,
    sat_vector: VectorSum,
    t0: Time,
    t1: Time,
    min_elevation: float,
) -> tuple:
    alt_initial, az_initial, ssb_dist = site_vector.at(t0).observe(sat_vector).apparent().altaz()
    alt_final, az_final, ssb_dist = site_vector.at(t1).observe(sat_vector).apparent().altaz()
    print(f"Start/End (Az/El): {az_initial.degrees:.5f} / {alt_initial.degrees:.5f} ----> {az_final.degrees:.5f} / {alt_final.degrees:.5f}")

    # This fails to find the "rise" above 10 degrees elevation.
    event_times, events = satellite.find_events(site, t0, t1, altitude_degrees=min_elevation)

    return event_times, events


def print_track(groundstation: VectorSum, satellite: EarthSatellite, track_time: Time) -> None:
    """Print a minute-by-minute track of a satellite."""
    print("Propagated track:")
    for obs_time in track_time:
        ssb_alt, ssb_az, ssb_dist = groundstation.at(obs_time).observe(satellite).apparent().altaz()
        print(f"{obs_time.utc_strftime('%Y-%m-%d %H:%M')}: {ssb_az.degrees:.5f} / {ssb_alt.degrees:.5f}")


if __name__ == "__main__":
    main()

The output, trimmed a bit for brevity:

Start/End (Az/El): 187.17892 / 3.56302 ----> 181.04224 / 31.75064
This should have found a 'rise':
[] <Time tt=[]>
...but it has not found a rise event.

Start/End (Az/El): 8.53929 / -5.12598 ----> 342.28175 / 80.02753
Expand the boundary, and it finds events:
[0 1] <Time tt=[2460614.73033378 2460614.79156628]>
...including a rise and culmination.

Start/End (Az/El): 8.53929 / -5.12598 ----> 153.81261 / 25.31236
Expand it further, and it should have found a rise, culmination, and set:
[0 1 1] <Time tt=[2460614.73033391 2460614.79156647 2460615.17820007]>
...but it found a rise, two culminates, and no set.

Propagated track:
2024-10-31 05:20: 192.82956 / -19.59808
2024-10-31 05:24: 189.81945 / -7.95510
2024-10-31 05:28: 187.23907 / 3.29044
2024-10-31 05:32: 184.96974 / 13.76284
2024-10-31 05:36: 182.93587 / 23.21467
2024-10-31 05:40: 181.08645 / 31.55621
2024-10-31 05:44: 179.38535 / 38.82250
2024-10-31 05:48: 177.80579 / 45.11913
2024-10-31 05:52: 176.32707 / 50.57753
2024-10-31 05:56: 174.93248 / 55.32799
2024-10-31 06:00: 173.60790 / 59.48724
2024-10-31 06:04: 172.34077 / 63.15455
2024-10-31 06:08: 171.11921 / 66.41212
2024-10-31 06:12: 169.93118 / 69.32710
2024-10-31 06:16: 168.76358 / 71.95414
2024-10-31 06:20: 167.60085 / 74.33764
2024-10-31 06:24: 166.42309 / 76.51380
2024-10-31 06:28: 165.20258 / 78.51223
2024-10-31 06:32: 163.89761 / 80.35721
2024-10-31 06:36: 162.43976 / 82.06872
2024-10-31 06:40: 160.70536 / 83.66305
2024-10-31 06:44: 158.44144 / 85.15312
2024-10-31 06:48: 155.03239 / 86.54776
2024-10-31 06:52: 148.50069 / 87.84690
2024-10-31 06:56: 128.13765 / 89.00029
2024-10-31 07:00: 39.69036 / 89.30654
2024-10-31 07:04: 3.38661 / 88.37322
2024-10-31 07:08: 354.56943 / 87.32850
2024-10-31 07:12: 350.60919 / 86.30500
2024-10-31 07:16: 348.25847 / 85.31556
2024-10-31 07:20: 346.63795 / 84.36104
2024-10-31 07:24: 345.41435 / 83.43988
2024-10-31 07:28: 344.43383 / 82.55000
2024-10-31 07:32: 343.61550 / 81.68929
2024-10-31 07:36: 342.91262 / 80.85572
2024-10-31 07:40: 342.29623 / 80.04743
2024-10-31 07:44: 341.74737 / 79.26273
2024-10-31 07:48: 341.25300 / 78.50008
2024-10-31 07:52: 340.80385 / 77.75808
2024-10-31 07:56: 340.39309 / 77.03548
2024-10-31 08:00: 340.01557 / 76.33114
2024-10-31 08:04: 339.66727 / 75.64402
2024-10-31 08:08: 339.34503 / 74.97318
2024-10-31 08:12: 339.04628 / 74.31778
2024-10-31 08:16: 338.76892 / 73.67704
...
2024-10-31 12:32: 339.18931 / 51.25362
2024-10-31 12:36: 339.32071 / 51.16099
2024-10-31 12:40: 339.45316 / 51.07786
2024-10-31 12:44: 339.58656 / 51.00441
2024-10-31 12:48: 339.72082 / 50.94083
2024-10-31 12:52: 339.85584 / 50.88734
2024-10-31 12:56: 339.99152 / 50.84416
2024-10-31 13:00: 340.12775 / 50.81152
2024-10-31 13:04: 340.26446 / 50.78967
2024-10-31 13:08: 340.40153 / 50.77889
2024-10-31 13:12: 340.53889 / 50.77946
2024-10-31 13:16: 340.67642 / 50.79169
2024-10-31 13:20: 340.81404 / 50.81590
2024-10-31 13:24: 340.95166 / 50.85244
2024-10-31 13:28: 341.08918 / 50.90169
2024-10-31 13:32: 341.22653 / 50.96405
2024-10-31 13:36: 341.36361 / 51.03994
2024-10-31 13:40: 341.50034 / 51.12985
2024-10-31 13:44: 341.63665 / 51.23427
2024-10-31 13:48: 341.77246 / 51.35374
2024-10-31 13:52: 341.90769 / 51.48884
...
2024-10-31 15:48: 346.96691 / 69.50775
2024-10-31 15:52: 347.66689 / 71.36404
2024-10-31 15:56: 348.66738 / 73.45936
2024-10-31 16:00: 350.20602 / 75.84171
2024-10-31 16:04: 352.83769 / 78.56846
2024-10-31 16:08: 358.20033 / 81.69377
2024-10-31 16:12: 13.60105 / 85.15166
2024-10-31 16:16: 81.85831 / 87.11603
2024-10-31 16:20: 134.10787 / 83.04728
2024-10-31 16:24: 145.82440 / 76.28422
2024-10-31 16:28: 150.17777 / 67.48542
2024-10-31 16:32: 152.25411 / 56.25051
2024-10-31 16:36: 153.31592 / 42.23893
2024-10-31 16:40: 153.80583 / 25.74685
2024-10-31 16:44: 153.90416 / 8.11710
2024-10-31 16:48: 153.67806 / -8.82336

I would prefer to continue using Skyfield for finding rise and set events, and I don't think I'm doing anything wrong. Is this a bug/limitation in Skyfield? Are there intentional limitations here that I should work around?

Thanks again!

[nils10]

Hi,
I ran into a similar issue with pass prediction regularly failing for the recently launched Proba-3 satellite (that is also in a highly elliptical orbit). In case you want to have another example case to reproduce, some info below:

TLE:

PROBA3-estimation
1 99991U 24999A   24340.44722222  .00000010  00000-0  52554-3 0  0017
2 99991  59.0000 142.9829 8111000 188.0000 000.0000 01.22265304000001

Observer location: 53°N, 5°E
Min. elevation: 0.0 deg
Start time (t0): 2024-12-07 12:35 UTC
End time (t1): 2024-12-09 12:35 UTC

Calling the following code to find all passes/events between the start and end time (2 days) results in the following events being returned:
t, events = self.satellite.find_events(self.observer, t0, t1, altitude_degrees=0.0)

Passes/Events:

Pass 1 (correct):
2024-12-06 21:40:24.319836+00:00 Event: 0 | EL: 0.0 deg CORRECT
2024-12-07 00:31:40.045713+00:00 Event: 1 | EL: 21.1 deg CORRECT
2024-12-07 01:26:46.400018+00:00 Event: 2 | EL: -0.0 deg CORRECT
Pass 2 (correct):
2024-12-07 04:11:32.986707+00:00 Event: 0 | EL: 0.0 deg CORRECT
2024-12-07 06:22:31.528189+00:00 Event: 1 | EL: 10.1 deg CORRECT
2024-12-07 09:27:28.588537+00:00 Event: 2 | EL: -0.0 deg CORRECT
Pass 3 (NOT correct):
2024-12-08 00:26:41.619490+00:00 Event: 0 | EL: 0.0 deg CORRECT
2024-12-08 05:15:41.914492+00:00 Event: 1 | EL: 30.6 deg CORRECT
<MISSING: here the 3rd pass actually has LOS at ~12:10:55 UTC>
<MISSING: here a 4th pass actually has AOS at ~15:33:55 UTC>
2024-12-08 16:53:26.786992+00:00 Event: 1 | EL: 16.3 deg (TME and EL is correct, but this is of a 4th pass, not 3rd pass!)
<MISSING: here the 4th pass actually has LOS at ~17:05:39 UTC>
<MISSING: here a 5th pass actually has AOS at ~22:56:48 UTC>
2024-12-09 05:25:31.332216+00:00 Event: 1 | EL: 45.2 deg (TME and EL is correct, but this is of a 5th pass, not 3rd pass!)
2024-12-09 12:31:59.786383+00:00 Event: 1 | EL: 88.1 deg (TME and EL is correct, but this is of a 5th pass, not 3rd pass! Its correct that the 5th pass indeed has 2 TME's.

So in this example the prediction has missed that passes have ended/new passes have started. It would be great if this could be fixed as its hard to know in production when this happens and can have significant impact when missed!

I found a comment in the sgp4lib.py that might be a starting point to fixing this:
# Long-period satellites might rise each day not because of their own motion, but because the Earth rotates under them, so check position at least each quarter-day. We might need to tighten this even further if experience someday shows it missing a pass of a particular satellite.

[brandon-rhodes]

Thanks for this bug report!

I have just landed a fix for your first issue, the missing rising event. I'll release a new version of Skyfield tomorrow, hopefully, with the fix. It makes your first output look like:

Start/End (Az/El): 187.17904 / 3.56301 ----> 181.04240 / 31.75064
This should have found a 'rise':
[0] <Time tt=[2460614.73032955]>
...but it has not found a rise event.

So I think it's fixed!

By the way, your script might be my favorite bug report test script ever! Some test scripts simply print a wall of numbers, and I need to wade through the person's explanation to learn where I should be looking in the output. But your little script explains, complete with arrows, exactly what output I'm supposed to look at and why!

Your third case, I think, is the result of a slight error in your bracketing. Your full track shows the satellite still above 10° at your end time of 16, 40, 6, so you need to adjust to 16, 44, 6 if you want to catch the setting time. Adjusting your line to look like this:

obs_long_end_time = datetime(2024, 10, 31, 16, 44, 6, 0, tzinfo=dt.timezone.utc)

—produces the correct output:

Start/End (Az/El): 8.53933 / -5.12596 ----> 153.90247 / 7.67826
Expand it further, and it should have found a rise, culmination, and set:
[0 1 1 2] <Time tt=[2460614.7303336775 ... 2460615.1977260807] len=4>
...but it found a rise, two culminates, and no set.

The reason for the two culminations will, I think, be clear if you review your diagram: the satellite passes the zenith not once, but twice. Maybe you could graph its altitude alone to see the two culminations.

Since I'm hopeful the bug here is fixed, I'm going to close the issue. If you want to go ahead and try the fix yourself before the actual release, you can install the development version of this project like this:

pip install -U https://github.com/skyfielders/python-skyfield/archive/master.zip

Please feel free to comment further here if there's a subtlety I haven't caught and there's still a problem. Thanks again for the bug report!

[ecdaero]

Excellent, thank you! I've tested out the latest release, and things do seem to work as expected now. Thanks for your patience, I did indeed goof on the 2 culminations and the end time, but I'm glad it didn't cause too much grief.

I'll let you know if we encounter any other unusual behavior.

[brandon-rhodes]

I'm glad it now works for you! Out of curiosity — what's the grey snake peeking into your plot from the left edge? Is it the Sun or Moon?

[ecdaero]

I believe that one shows the sun, yeah. We've written an internal tool to plan tracks on our telescope ahead of time, and have a separate tool that shows the track when it gets scheduled on our telescope.