Getting "RuntimeWarning: invalid value encountered in double_scalars" in keplerlib.py

[Bernmeister]

Running the following code:

#!/usr/bin/env python3

import datetime, io, skyfield.api, skyfield.constants, skyfield.data.mpc

print( "Skyfield:", skyfield.__version__ )
latitude = -33
longitude = 151
now = datetime.datetime.strptime( "2020-08-08", "%Y-%m-%d" )
data = "    CK19Y04a  2020 05 31.0420  0.251014  1.001333  177.2464  120.9277   45.8250  20200807  11.6  6.0  C/2019 Y4-A (ATLAS)                                      MPEC 2020-L06"
timeScale = skyfield.api.load.timescale( builtin = True )
topos = skyfield.api.Topos( latitude_degrees = latitude, longitude_degrees = longitude )
ephemeris = skyfield.api.load( "de421.bsp" )
sun = ephemeris[ "sun" ]
earth = ephemeris[ "earth" ]
with io.BytesIO( data.encode() ) as f:
    dataframe = skyfield.data.mpc.load_comets_dataframe( f ).set_index( "designation", drop = False )
    body = sun + skyfield.data.mpc.comet_orbit( dataframe.loc[ "C/2019 Y4-A (ATLAS)" ], timeScale, skyfield.constants.GM_SUN_Pitjeva_2005_km3_s2 )

I am getting the following result:

Skyfield: 1.26
/usr/local/lib/python3.5/dist-packages/skyfield/keplerlib.py:256: RuntimeWarning: invalid value encountered in double_scalars
  return 2 * arctan(((1+e)/(1-e))**.5 * tan(E/2))

I stepped into the function true_anomaly and the parameter E is 0.0...not sure where to go from here!

[brandon-rhodes]

I started commenting out parts of the expression, and happily E is not involved; if I run tan(E/2) by itself then no warning is printed.

Trying instead the first half of the expression, I get the warning only if **.5 is applied to the computed fraction, so I printed the value of the fraction and got:

    print((1+e)/(1-e))

-1501.375093773411

My guess is that this formula only applies for closed, circular or elliptical orbits, but not for hyperbolic ones. @JoshPaterson wrote keplerlib.py so he might have some insight, but I'll at least go do a quick web search to see if there's an alternative formula for a hyperbolic orbit.

[brandon-rhodes]

Well, that was pretty quick. I wound up copying and pasting a formula from the Wikipedia. Could you try installing Skyfield from master to see if it fixes your problem?

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

Thanks for any feedback you can provide!

[Bernmeister]

To test, I downloaded that zip file, opened up and dropped the new keplerlib.py into /usr/local/lib/python3.6/dist-packages/skyfield after I had renamed the existing keplerlib.py to keplerlib.py.ORIGINAL (wasn't sure how the pip would work nor how to back out). Assuming that's cool, I ran the test and it worked.

So I then ran the test over the full comet set from MPC, and got an exception out in mpc.py:

/usr/local/lib/python3.6/dist-packages/skyfield/data/mpc.py:207: RuntimeWarning: divide by zero encountered in double_scalars
  a = row.perihelion_distance_au / (1.0 - e)
/usr/local/lib/python3.6/dist-packages/skyfield/data/mpc.py:208: RuntimeWarning: invalid value encountered in double_scalars
  p = a * (1.0 - e*e)
/usr/local/lib/python3.6/dist-packages/skyfield/keplerlib.py:242: RuntimeWarning: invalid value encountered in double_scalars
  dE = dM/(1 - e*cos(E))
Traceback (most recent call last):
  File "magnitudeTest.py", line 17, in <module>
    body = sun + skyfield.data.mpc.comet_orbit( dataframe.loc[ "C/2015 A2 (PANSTARRS)" ], timeScale, skyfield.constants.GM_SUN_Pitjeva_2005_km3_s2 )
  File "/usr/local/lib/python3.6/dist-packages/skyfield/data/mpc.py", line 222, in comet_orbit
    row['designation'],
  File "/usr/local/lib/python3.6/dist-packages/skyfield/keplerlib.py", line 156, in _from_mean_anomaly
    E = eccentric_anomaly(eccentricity, M)
  File "/usr/local/lib/python3.6/dist-packages/skyfield/keplerlib.py", line 247, in eccentric_anomaly
    raise ValueError('Failed to converge')
ValueError: Failed to converge

However the comet data which caused the exception was this:

"    CK15A020  2015 08  1.8353  5.341055  1.000000  208.8369  258.5042  109.1696            10.5  4.0  C/2015 A2 (PANSTARRS)                                    MPC 93587"

So to test, use the original code at the top but swap out line 9 for the data and line 17 with the name of the comet ("C/2015 A2 (PANSTARRS)").

[brandon-rhodes]

That sounds like a quite separate issue, but while we're here, let's just re-open and continue with that new exception. I should be able to take a look in the next day or two!

[brandon-rhodes]

Drat, when I said "re-open" I meant that I was going to press the "Reopen" button right here so we could continue the conversation without your having to open a new issue.

Then?

I closed the tab and entirely forgot to accidentally press "Reopen".

Thanks for opening the new issue!

[JoshPaterson]

I’ll be able to take a look at this later tonight!