Test with JWST TLE data

[mworion]

@brandon-rhodes,

we already had the situation with an downcoding starling satellite which resulted in runtime warnings in skyfield. Right now I tried to test the JWST telescope TLE data with skyfield and got the same issue. Where am I wrong in doing so?

If I try heavens above (, the display the path based??? on the give TLE parameters.

Attached a short script for testing:

from skyfield.api import load, wgs84, EarthSatellite

def findSatUp(sat, loc, tStart, tEnd, alt):
    t, events = sat.find_events(loc, tStart, tEnd, altitude_degrees=alt)
    return t, events

def main():
    ts = load.timescale()
    tStart = ts.now()
    tEnd = ts.tt_jd(tStart.tt + 12 / 86400)
    loc = wgs84.latlon(latitude_degrees=49, longitude_degrees=-11)
    tle = ["JWST",
           "1 50463U 21130A   21362.00000000  .00000000  00000-0  00000-0 0  9999",
           "2 50463   4.6198  89.0659 9884983 192.3200  17.4027  0.01958082    27"]
    sat = EarthSatellite(tle[1], tle[2],  name=tle[0])
    t, events = findSatUp(sat, loc, tStart, tEnd, 10)
    print(t, events)

if __name__ == "__main__":
    main()

Many thanks in advance for looking at it.

Michel

[brandon-rhodes]

I have tried out your script; thank you for providing one that could immediately be run and tried out. Is this the error you are asking about (so that I can be sure I am looking at the right one)?

/home/brandon/skyfield/skyfield/searchlib.py:184: RuntimeWarning: invalid value encountered in less
  indices = flatnonzero(dsd < 0)

[mworion]

I See the following:

/Users/mw/PycharmProjects/Envs/mw4/lib/python3.8/site-packages/skyfield/functions.py:85: RuntimeWarning: invalid value encountered in remainder
  phi = arctan2(y, x) % tau

Perhaps this might be related to zero's in fields 9, 10, 11 (First derivative of mean motion; the ballistic coefficient, Second derivative of mean motion, the drag term, or radiation pressure coefficient) of the two line elements. It may also related how SGP4 treats these elements. Still my knowledge about this is still limited.

Michel

Edit: I made numpy raise the full stack:

Traceback (most recent call last):
  File "/Users/mw/PycharmProjects/MountWizzard4/gists/brandon_runtime_warning2.py", line 26, in <module>
    main()
  File "/Users/mw/PycharmProjects/MountWizzard4/gists/brandon_runtime_warning2.py", line 21, in main
    t, events = findSatUp(sat, loc, tStart, tEnd, 10)
  File "/Users/mw/PycharmProjects/MountWizzard4/gists/brandon_runtime_warning2.py", line 8, in findSatUp
    t, events = sat.find_events(loc, tStart, tEnd, altitude_degrees=alt)
  File "/Users/mw/PycharmProjects/Envs/mw4/lib/python3.8/site-packages/skyfield/sgp4lib.py", line 250, in find_events
    tmax, altitude = find_maxima(t0, t1, altitude_at, half_second, 12)
  File "/Users/mw/PycharmProjects/Envs/mw4/lib/python3.8/site-packages/skyfield/searchlib.py", line 141, in find_maxima
    y = f(t)
  File "/Users/mw/PycharmProjects/Envs/mw4/lib/python3.8/site-packages/skyfield/sgp4lib.py", line 247, in altitude_at
    return at(t).altaz()[0].degrees
  File "/Users/mw/PycharmProjects/Envs/mw4/lib/python3.8/site-packages/skyfield/positionlib.py", line 609, in altaz
    return _to_altaz(self, temperature_C, pressure_mbar)
  File "/Users/mw/PycharmProjects/Envs/mw4/lib/python3.8/site-packages/skyfield/positionlib.py", line 835, in _to_altaz
    r_au, alt, az = to_spherical(position_au)
  File "/Users/mw/PycharmProjects/Envs/mw4/lib/python3.8/site-packages/skyfield/functions.py", line 85, in to_spherical
    phi = arctan2(y, x) % tau
FloatingPointError: invalid value encountered in remainder

[brandon-rhodes]

Alas, on my machine your script still prints:

/home/brandon/skyfield/skyfield/searchlib.py:184: RuntimeWarning: invalid value encountered in less
  indices = flatnonzero(dsd < 0)
<Time tt=[]> []

— instead of the error that you are seeing. Thank you for the full traceback from your machine. I think your guess might be correct that this satellite's elements are generating invalid values. One thing you might try is generating positions over the whole period you want to study, and check to see if any of them are nan.

If I could get the error to show up here on my own laptop, the first thing I would probably try would be print() statements in functions.py right above the line where the error is being triggered, to see what y and x look like. But, of course, you might be less happy about editing the source of an installed library? I once made such an edit and for the following year had a debug message occasionally printing to the screen because I forgot which library I had edited.

[davidmikolas]

perhaps helpful, perhaps not Why can n2yo somehow interpret this TLE but Skyfield can't? Why does it return nans and zeros for position?

Position relative to the geocenter +/- 2 days from epcoch returns a mix of nans and zeros.

I'm using Skyfield version 1.41

Somehow n2yo's implementation can still propagate this TLE. https://www.n2yo.com/satellite/?s=50463 minutes before the time of this post

import numpy as np
import matplotlib.pyplot as plt
from skyfield.api import load, EarthSatellite, Loader, wgs84
import datetime

loaddata = Loader('~/Documents/fishing/SkyData')  # avoids multiple copies of large files
ts = loaddata.timescale() # include builtin=True if you want to use older files (you may miss some leap-seconds)
eph = loaddata('de421.bsp') # a small one, fine for this

L1, L2 = """1 50463U 21130A   21362.00000000  .00000000  00000-0  00000-0 0  9999
2 50463   4.6198  89.0659 9884983 192.3200  17.4027  0.01958082    27""".splitlines()

sat = EarthSatellite(L1, L2)

days = np.arange(17)/4. - 2.  # +/- 2 days in 6 hour increments

times = ts.tt_jd(sat.epoch.tt + days)

epoch_string = datetime.datetime(*sat.epoch.utc[:5]).isoformat().replace('T', ' ') + 'UTC'

print('+/- 2 days around epoch of: ', epoch_string)

g = sat.at(times) # geocentric position object

p = g.position.km  # geocentric positions (km)

d = g.distance().km

if True:
    fig, ax = plt.subplots(1, 1)
    ax.plot(days, d)
    ax.plot(days, d, 'ok')
    ax.set_xlim(days.min() - 0.1, days.max()+0.1)
    ax.set_xlabel('days relative to TLE epoch')
    ax.set_ylabel('distance to geocenter (km)')
    ax.set_xlabel('days relative to TLE epoch')
    plt.suptitle('JWST TLE epoch: ' + epoch_string)
    plt.show()

[brandon-rhodes]

@davidmikolas — Could you verify whether those n2yo positions are accurate, and are providing a useful position for the satellite — maybe whether they are at least showing the correct lat/lon subpoint on their map? The Python wrapper for the raw C sgp4 routine currently does this:

    if (satrec.error && satrec.error < 6)
        r[0] = r[1] = r[2] = v[0] = v[1] = v[2] = NAN;

— to avoid having nonsense values returned. If you will print(g.message) in your sample script, you will see that all of those positions are failing to compute because the perturbed eccentricity is outside the range 0.0 to 1.0. Maybe the n2yo web site runs the C version of SGP4 but then ignores error conditions?

[brandon-rhodes]

@davidmikolas — Oh! I should have pointed out that in the case of that message, satrec.error equals 3, so the if statement triggers and runs its statement body.

[davidmikolas]

@brandon-rhodes oh the n2yo plot is of course nonsense which is why the Space SE question only asks " Why can n2yo somehow interpret this TLE..." then goes on to explain why this TLE shouldn't mean anything.

The answer is of course "Because they're not using Skyfield (yet!)" :-)

set(g.message) returns {'perturbed eccentricity is outside the range 0.0 to 1.0'} so while some positions are nans and some are zeros, the error message thinks both may be due to the same situation.

All nans might be better, returning real numbers that might be used in downstream calculations (e.g. subtract 6378 km from the distance then do something else to it) might never catch that there's a problem, whereas returning nans would certainy get the user's attention that something is amiss.

In the mean time please feel free to post an answer there about how much more care Skyfield takes to make sure we don't get crazy results and how n2yo might not.

[cruisen]

Maybe I am missing a point here, but as far as I know TLE can only be "regularly fitted" in order to sort of represent trajectories of L2 (or for that matter any Lagrange point) Halo or Lissajous "orbit". TLE are not applicable.

But for example for JWST, a better source is: https://ssd.jpl.nasa.gov/horizons/app.html#/

And as said, I might be missing the point here, and if I did, please comment, so I can turn this into a learning opportunity.

More on that matter: https://github.com/skyfielders/python-skyfield/issues/572

[davidmikolas]

@cruisen you've got it right, TLEs are not for things that aren't bound to Earth and orbit it. However sometimes Project Pluto generates TLEs for extreme/unusual situations as simply a way to get an idea roughly where something is, and the n2yo site gathers TLEs from several places besides official sources.

See answers to https://space.stackexchange.com/questions/24314/has-a-tle-ever-been-issued-for-a-spacecraft-trajectory-not-bound-to-earth-orbit

Also https://space.stackexchange.com/questions/57771/does-the-james-webb-telescope-have-a-two-line-element-set

It turns out that space-track.org has two official TLEs for 2021-130A issued at epochs 21359.67 and 21362.00 so as far as I can tell, these appears to be good candidates for TLEs of objects not strictly bound to Earth. During this period JWST was traveling along a trajectory close to the 3-body stable manifold associated with Sun-Earth L2. This is not a situation where a TLE would be reliable, but it seems the TLEs were issues as "helpful hints" to give folks at least some handle on where it was during this time period, since Horizons doesn't update their ephemeris in real time.

The folks at JPL do act fairly quickly, and once sufficiently reliable tracking data is available through the Deep Space Network they will do some numerical orbit propagation that is fitted to DSN tracking data points, then update their ephemeris accordingly. It will take a while but you can see how the comments and answers to https://space.stackexchange.com/questions/57822/is-jwst-halo-orbit-prograde-or-retrograde-and-why evolved over time as JPL expanded their coverage of their JWST ephemeris.

This is a really, really, really interesting situation! I will update here if I find out more about the issuing of those two TLEs for a spacecraft not bound to Earth.

update: https://space.stackexchange.com/questions/58319/did-18-spcs-issue-two-tles-for-jwst-when-it-was-not-bound-to-earth-nor-orbiting

[brandon-rhodes]

Please feel free to continue discussion here around interesting topics like L2 orbits, I am enjoying following all the links! But I am going to go ahead and mark this issue as closed, as I think the circumstances of the original question are clear: the SGP4 module returns NaN values when an error is raised, and older versions of NumPy print various warnings to the screen when Skyfield's position logic tries doing math with these NaN values.

I wish NumPy would just be quiet, and let NaN inputs produce NaN outputs as Nature intended. Indeed, the most recent NumPy versions seem far quieter, so the next time this is asked I might simply ask folks to try updating NumPy to the latest version.

Also, I'll be quicker next time to point users to a check of the error messages being produced by the position.

It feels like the NumPy warning/error settings isn't something Skyfield should be reaching in and adjusting on the user's behalf, but I'm certainly open to updating the documentation, if someone has figured out how a script can ask NumPy to quiet down before they attempt a computation that might involves some dates on which positions come out as NaN.

Again, feel free to comment further even as this is marked Closed, especially if you think there's action here for Skyfield to take, that I've missed. Thanks!