Light Convergence Failure

[meawoppl]

Heyo @brandon-rhodes thanks for skyfield!

We do a lot of satellite math using your package, and I bumped into this gem while working with it. Here is a reproduction case:

def test_convergence_demo():
    tle_0 = """STARLINK-1885"""
    tle_1 = """1 47168U 20088AY  22221.92842503  .01483002  11676-2  34717-2 0  9991"""
    tle_2 = """2 47168  53.0385 305.4734 0002390 239.0067 121.2198 15.96617970 96526"""

    load
    satellite = EarthSatellite(tle_1, tle_2, tle_0)
    capp = wgs84.latlon(37.76115,-122.418056,12)
    ts = load.timescale()
    time = ts.tt_jd(2459831.9170483043)

    _planets = load('de421.bsp')

    # ssb = Solar System Barycentric
    location_ssb = _planets["earth"] + capp
    target_ssb = _planets["earth"] + satellite

    # NOTE(meawoppl) - altaz accepts temp/humidity correction...
    location_ssb.at(time).observe(target_ssb).apparent().altaz()

This code takes a road less traveled and reports the following:

observer = <Barycentric BCRS position and velocity at date t center=0 target=WGS84 latitude +37.7612 N longitude -122.4181 E elevation 12.0 m>
target = <VectorSum of 3 vectors:
 'de421.bsp' segment 0 SOLAR SYSTEM BARYCENTER -> 3 EARTH BARYCENTER
 'de421.bsp' segment 3 EARTH BARYCENTER -> 399 EARTH
 EarthSatellite 399 EARTH -> STARLINK-1885 catalog #47168 epoch 2022-08-09 22:16:56 UTC>

    def _correct_for_light_travel_time(observer, target):
        """Return a light-time corrected astrometric position and velocity.

        Given an `observer` that is a `Barycentric` position somewhere in
        the solar system, compute where in the sky they will see the body
        `target`, by computing the light-time between them and figuring out
        where `target` was back when the light was leaving it that is now
        reaching the eyes or instruments of the `observer`.

        """
        t = observer.t
        ts = t.ts
        whole = t.whole
        tdb_fraction = t.tdb_fraction

        cposition = observer.position.au
        cvelocity = observer.velocity.au_per_d

        tposition, tvelocity, gcrs_position, message = target._at(t)

        distance = length_of(tposition - cposition)
        light_time0 = 0.0
        for i in range(10):
            light_time = distance / C_AUDAY
            delta = light_time - light_time0
            if max(abs(delta)) < 1e-12:
                break

            # We assume a light travel time of at most a couple of days.  A
            # longer light travel time would best be split into a whole and
            # fraction, for adding to the whole and fraction of TDB.
            t2 = ts.tdb_jd(whole, tdb_fraction - light_time)

            tposition, tvelocity, gcrs_position, message = target._at(t2)
            distance = length_of(tposition - cposition)
            light_time0 = light_time
        else:
>           raise ValueError('light-travel time failed to converge')
E           ValueError: light-travel time failed to converge

This case was hard to find... I suspect we have called an equivalent of this subroutine maybe... a few million times without experiencing it....

Happy to dig in more if there isn't something obvious that springs to mind.

[meawoppl]

Also, if you ever find yourself at or near that lat/long, I owe you at least a couple drinks :beers:

[rmathar]

A shot in the dark: This might be one of the rare cases where the solutions are oscillatory/bi-stable? Probably this happens when there is more than one position with almost the same light delay, perhaps when the satellite is right above the observer's head at that point in time [or at the antipodal point, one earth diameter plus 300 km of altitude corresponding to 15.9 rounds per day of the LTE]. The positions of constant time delay are sets of spheres centered at the observer, and if the satellite just scrapes past one of these (almost tangentially) there might be 2 close positions...

Printing the position of the satellite for each of the 10 iterations ought to reveal that kind of quiver/oscillation. The usual recipe in the software would be not to take an absolute value as the convergence criterion but to leave the loop if subsequent updates of the delta don't decrease anymore.

[meawoppl]

That was my suspicion. Again this is a 1/1000000 case or so, but I took the time to round it up and make a repro as it might be helpful in testing/triage for others.

I suspect the "right" thing to do here is roughly:

• have a convergence criteria, define behavior for the fail case
• then use a more robust optimizer?

[brandon-rhodes]

I think the problem is that the Starlink satellite is moving at super-relativistic speed, and therefore images of the satellite are arriving at the viewer in reverse order—by the time its distant image has reached you, it's already close enough to Earth to have delivered several more recent images, so that the observer will actually see it appear to travel backwards in time as, looking up, you see the satellite appearing to move away as more and more distant images finally reach you.

To wit:

tle_0 = """STARLINK-1885"""
tle_1 = """1 47168U 20088AY  22221.92842503  .01483002  11676-2  34717-2 0  9991"""
tle_2 = """2 47168  53.0385 305.4734 0002390 239.0067 121.2198 15.96617970 96526"""

from skyfield.api import EarthSatellite, load
from skyfield.functions import length_of

c_km_per_s = 299792.458

satellite = EarthSatellite(tle_1, tle_2, tle_0)
ts = load.timescale()
t = ts.tt_jd(2459831.9170483043)
xyz_previous = satellite.at(t).xyz
one_second = 1/24.0/3600.0
for i in range(10):
    t = t + one_second
    xyz = satellite.at(t).xyz
    range_km = length_of(xyz.km)
    moved_km = length_of(xyz.km - xyz_previous.km)
    print(
        t.utc_strftime(),
        ' Range from Earth (km): {:.1f}'.format(range_km),
        ' Speed (c): {:.3f}'.format(moved_km / c_km_per_s),
    )
    xyz_previous = xyz

Which produces:

2022-09-09 09:59:25 UTC  Range from Earth (km): 1685123.5  Speed (c): 8.761
2022-09-09 09:59:26 UTC  Range from Earth (km): 1613951.1  Speed (c): 8.486
2022-09-09 09:59:27 UTC  Range from Earth (km): 1547649.9  Speed (c): 8.935
2022-09-09 09:59:28 UTC  Range from Earth (km): 1485763.7  Speed (c): 2.711
2022-09-09 09:59:29 UTC  Range from Earth (km): 1427889.8  Speed (c): 5.380
2022-09-09 09:59:30 UTC  Range from Earth (km): 1373671.8  Speed (c): 8.900
2022-09-09 09:59:31 UTC  Range from Earth (km): 1322793.0  Speed (c): 6.029
2022-09-09 09:59:32 UTC  Range from Earth (km): 1274971.2  Speed (c): 3.727
2022-09-09 09:59:33 UTC  Range from Earth (km): 1229954.1  Speed (c): 2.316
2022-09-09 09:59:34 UTC  Range from Earth (km): 1187515.8  Speed (c): 6.507

Check my math here—scripts that I write after midnight have been known to have impressive flaws—but it looks like the satellite's speed second-to-second is varying between about 2× the speed of light and 8× the speed of light. Under such conditions a convergence won't be found.

‘Luke, at that speed do you think you'll be able to pull out in time?’

[rmathar]

So the problem appears to be that the program is trying to observe a satellite on 2022-09-09 which decayed in 2022-08-14 and was for that end-of-life starting an extreme (non-physical) elliptic orbit approximately 2 days earlier while actually crashing into the atmosphere: https://isstracker.pl/en/satelity/47168 . The extreme LTE parameter set to 0.0148 (on 2022-08-10) for the change in the n-parameter is in some way indicating that sort of problem.

[meawoppl]

That appears to be quite non-physical! As an optical imaging group, I think our chances of catching super-relativistic objects on camera are pretty low :wink:

I appreciate the time/triage effort you spent here!

Might it make sense to check the value of tvelocity for this then balk explicitly? Happy to open a PR if you think it would be helpful.

[brandon-rhodes]

Might it make sense to check the value of tvelocity for this then balk explicitly? Happy to open a PR if you think it would be helpful.

Since the caller might be asking about n dates, I actually would like to back away from the ValueError that Skyfield currently raises, and probably just drop back to returning nan if the result doesn't converge, so that if there are multiple dates involved maybe some of them could still return good results. My guess is that for now a pull request might not be worth it, as this is the only time this problem has ever come up?