Invalid value encountered when propagating satellite and computing latitude/longitude

[AshaJain21]

Hello, I am propagating reentered satellites from their latest TLE to their predicted reentry time, and for several satellites I get the error: invalid value encountered in double_scalars lon = (arctan2(y, x) - pi) % tau - pi. I believe this error may be because aractan2 is returning a nan value (possibly because x and y are nonsensical) but I am unsure why x and y are invalid or how to fix it. I have checked that the TLE's are within 5 days of the predicted reentry time.

For example, the error occurs when propagating Starlink-2638 (catalog #48477 epoch 2021-06-02 21:34:18 UTC) to June 3rd 2021 12:19:00 UTC, using the following python code:

ts = load.timescale()
line_1 = '1 48477U 21040BB  21153.89881469  .33258120  12317-4  10940-2 0  9990'
line_2 = '2 48477  53.0189  95.7131 0009567 170.8802 189.9362 16.43956428  5244'
satellite = EarthSatellite(line_1, line_2, 'Starlink-2638', ts)
print(satellite)
t = ts.utc(2021, 6, 3, 12, 19, 0)
days = t - satellite.epoch
print('{:.3f} days away from epoch'.format(days))
geocentric = satellite.at(t)
lat_true, lon_true = wgs84.latlon_of(geocentric)
print("True lat: " + str(lat_true.degrees))
print("True lon: " + str(lon_true.degrees))

The output of the code above:

Starlink-2638 catalog #48477 epoch 2021-06-02 21:34:18 UTC
0.614 days away from epoch
True lat: nan
True lon: nan

The full error is:

/Users/ashajain/anaconda3/lib/python3.8/site-packages/skyfield/toposlib.py:195: RuntimeWarning: invalid value encountered in double_scalars
 lon = (arctan2(y, x) - pi) % tau - pi

Thanks for the help!

[brandon-rhodes]

1. I suggest that you plot the satellite's altitude over the period of time in question. My guess is that the nan values are (correctly) denoting a time at which the satellite no longer has a valid position.
2. But the error message is indeed annoying. Skyfield should hopefully avoid such messages in cases where nan is a meaningful result. I'll go see if the message can be turned off.

[brandon-rhodes]

An update: I can't get that error to appear here on my machine. (I thought for a moment I'd seen it, but it was an unrelated warning.) On my version of NumPy, arctan(nan, nan) comes out silently and cleanly as nan without any unnecessary drama. Even if:

• I set the PYTHONWARNINGS=d environment variable.
• And do import numpy followed by numpy.seterr('raise').

@AshaJain21 — Let's compare NumPy versions. Is yours newer or older than numpy==1.15.4?

[AshaJain21]

Hi Brandon, thanks for the quick replies. I have numpy 1.20.1 according to pip list in my juypter notebook.TLEtoLatLong.pdf In testing some more, I get a silent failing sometimes when I initialize an EarthSatellite from a tle directly but always get the warning explicitly when I initialize an EarthSatellite with an explicit definition of SatRec() as shown below:

satrec = Satrec()
    satrec.sgp4init(
        WGS72,           # gravity model
        'i',             # 'a' = old AFSPC mode, 'i' = improved mode
        int(satellite_info['norad'].values[0]),               # satnum: Satellite number
        epoch_difference_days,       # epoch: days since 1949 December 31 00:00 UT
        satellite_info['bstar'].values[0], # bstar: drag coefficient (/earth radii)
        satellite_info['dn_o2'].values[0], # ndot: ballistic coefficient (revs/day)
        satellite_info['ddn_o6'].values[0], # nddot: second derivative of mean motion (revs/day^3)
        satellite_info['ecc'].values[0],       # ecco: eccentricity[TLEtoLatLong.pdf](https://github.com/skyfielders/python-skyfield/files/7864484/TLEtoLatLong.pdf)        satellite_info['argp'].values[0]* (np.pi/180), # argpo: argument of perigee (radians)
        satellite_info['inc'].values[0]* (np.pi/180), # inclo: inclination (radians)
        satellite_info['M'].values[0]* (np.pi/180), # mo: mean anomaly (radians)
        satellite_info['n'].values[0]* (2*np.pi) / (24*60), # no_kozai: mean motion (radians/minute)
        satellite_info['raan'].values[0]* (np.pi/180), # nodeo: right ascension of ascending node (radians)
    )
    satellite = EarthSatellite.from_satrec(satrec, ts)

Attached is a copy of the juypter notebook exported to a .pdf file (thus the file is not run-ready)

[brandon-rhodes]

I'm happy that you've discovered a reliable way to induce the warning to print! Could you provide a small script that does so? Without a satellite_info data structure I can't run you example code.

[AshaJain21]

Definitely, attached is a simple script. Note that I had to add the .txt extension to allow the file to attach, so please remove that extension. I was not able to test and run the file due to current issues on my computer, but I will upload the output once I am able to.

propogation_test.py.txt