Slight speedup in Kepler solver

[rmathar]

An approx. 20 percent speedup (measured accross a regular grid of 100x100 values in 0<M<pi and 0<e<1) in the skyfiled/keplerlib.py solver for Keplerian elliptic orbits is obtained by substituting eccentric_anomaly by eccentric_anomalyFast proposed below. The full test source looks like this:

import math
from numpy import(cos, sin,sqrt, sign)

def normpi(num):
    return (num + math.pi) % (2*math.pi)-math.pi

def eccentric_anomaly(e, M):
    """ Iterates to solve Kepler's equation to find eccentric anomaly
    Based on the algorithm in section 8.10.2 of the Explanatory Supplement
    to the Astronomical Almanac, 3rd ed.
    """
    M = normpi(M)
    E = M + e*sin(M)

    max_iters = 100
    iters = 0
    while iters < max_iters:
        dM = M - (E - e*sin(E))
        dE = dM/(1 - e*cos(E))
        E = E + dE
        if abs(dE) < 1e-14: return E
        iters += 1
    else:
        raise ValueError('Failed to converge')

def eccentric_anomalyFast(e, M):
    """ Iterates to solve Kepler's equation to find eccentric anomaly
    See arXiv:2108.03215
    """
    M = normpi(M)
    ebar = 0.25*math.pi/e-1.0
    E = 0.5*math.pi*ebar*(sign(ebar)*sqrt(1+M/(e*ebar*ebar))-1.0)

    max_iters = 10
    iters = 0
    while iters < max_iters:
        f2 = e*sin(E)
        f = E-f2-M
        f1 = 1.0-e*cos(E)
        dE = f*f1/(f1*f1-0.5*f*f2)
        E -= dE
        if abs(dE) < 1e-14:
            return E
        iters = iters+1
    else:
        raise ValueError('Failed to converge')

# numerical check: run the original and faster version and return the difference...
def comp(e,M):
    orig = eccentric_anomaly(e,M)
    rjm = eccentric_anomalyFast(e,M)
    # print( eccentric_anomaly(e,M),eccentric_anomalyFast(e,M))
    return rjm-orig

def compareVersions():
    for iM in range(1,50):
        M = math.pi*iM/50.0
        for ie in range(1,50):
            e = ie/50.0
            print(e,M,comp(e,M))

# compareVersions()

from datetime import datetime

# speed test: run original and faster code on a 100x100 grid of (M,e) values
def compareTimes():
    torig = datetime.now()
    for iM in range(1,100):
        M = math.pi*iM/100.0
        for ie in range(1,100):
            e = ie/100.0
            eccentric_anomaly(e,M)
    torig = datetime.now()-torig
    print(torig)

    tfast = datetime.now()
    for iM in range(1,100):
        M = math.pi*iM/100.0
        for ie in range(1,100):
            e = ie/100.0
            eccentric_anomalyFast(e,M)
    tfast = datetime.now()-tfast
    print( tfast)

compareTimes()

[brandon-rhodes]

Thank you for the suggestion, Dr. Mathar! I've replaced the old convergence routine with yours, and it will be used by the next version of Skyfield:

423ddb3

[Bernmeister]

I suspect I've hit an issue with the above commit. Given

#!/usr/bin/env python3


from skyfield.api import load, wgs84
from skyfield.constants import GM_SUN_Pitjeva_2005_km3_s2
from skyfield.data import mpc

import datetime


ephemeris = load( "de421.bsp" )

latitude = -33
longitude = 151
elevation = 100

timeScale = load.timescale( builtin = True )
utcNow = datetime.datetime.utcnow()
t = timeScale.utc( utcNow.year, utcNow.month, utcNow.day, utcNow.hour, utcNow.minute, utcNow.second )

sun, earth = ephemeris[ "sun" ], ephemeris[ "earth" ]
locationAtNow = ( earth + wgs84.latlon( latitude, longitude, elevation ) ).at( t )
sunAtNow = sun.at( t )

with load.open( "MPCORB-5.DAT" ) as f:
    minorPlanets = mpc.load_mpcorb_dataframe(f)

bad_orbits = minorPlanets.semimajor_axis_au.isnull()
minorPlanets = minorPlanets[ ~bad_orbits ]
print( minorPlanets.shape[ 0 ], "minor planets loaded." )

minorPlanets = minorPlanets.set_index( "designation", drop = False )

for name, row in minorPlanets.iterrows():
    minorPlanet = sun + mpc.mpcorb_orbit( row, timeScale, GM_SUN_Pitjeva_2005_km3_s2 )
    ra, dec, distanceToEarth = earth.at( t ).observe( minorPlanet ).radec()
    print( name, ra, dec )

I took the first 5 entries from MPCORB.DAT and saved to the file MPCORB-5.DAT:

00001    3.31  0.15 K221L 291.37562   73.63703   80.26858   10.58769  0.0785010  0.21424745   2.7660431  0 MPO681823  7301 122 1801-2022 0.65 M-v 30k Pan        0000      (1) Ceres              20220126
00002    4.11  0.15 K221L 272.47992  310.69724  172.91658   34.92531  0.2299930  0.21366046   2.7711069  0 MPO681823  8875 119 1804-2022 0.58 M-c 28k Pan        0000      (2) Pallas             20220105
00003    5.28  0.15 K221L 261.29865  247.94175  169.84777   12.99186  0.2568870  0.22606442   2.6687911  0 MPO660785  7440 112 1804-2021 0.61 M-v 3Ek Pan        0000      (3) Juno               20211107
00004    3.40  0.15 K221L   7.03155  151.09088  103.80396    7.14172  0.0882342  0.27163553   2.3612665  0 MPO660785  7532 107 1821-2021 0.62 M-p 18k Pan        0000      (4) Vesta              20210903
00005    6.98  0.15 K221L 160.98208  358.74039  141.57036    5.36762  0.1900994  0.23850283   2.5751766  0 MPO660785  3175  84 1845-2021 0.70 M-v 3Ek Pan        0000      (5) Astraea            20211105

When I run with version 1.43.1, I get the output

5 minor planets loaded.
(1) Ceres 09h 51m 12.91s +19deg 43' 36.4"
(2) Pallas 06h 11m 59.08s -10deg 41' 44.5"
(3) Juno 22h 55m 21.76s -05deg 29' 30.3"
(4) Vesta 22h 00m 17.02s -22deg 19' 43.6"
(5) Astraea 23h 06m 42.75s -09deg 24' 35.7"

But when I run with version 1.44, I get

5 minor planets loaded.
(1) Ceres 09h 51m 13.14s +19deg 43' 35.5"
/usr/local/lib/python3.8/dist-packages/skyfield/keplerlib.py:281: RuntimeWarning: invalid value encountered in sqrt
  E = 0.5 * pi * ebar * (sign(ebar) * sqrt(1 + M/(e*ebar*ebar)) - 1.0)
Traceback (most recent call last):
  File "test.py", line 35, in <module>
    minorPlanet = sun + mpc.mpcorb_orbit( row, timeScale, GM_SUN_Pitjeva_2005_km3_s2 )
  File "/usr/local/lib/python3.8/dist-packages/skyfield/data/mpc.py", line 97, in mpcorb_orbit
    minor_planet = _KeplerOrbit._from_mean_anomaly(
  File "/usr/local/lib/python3.8/dist-packages/skyfield/keplerlib.py", line 194, in _from_mean_anomaly
    E = eccentric_anomaly(eccentricity, M)
  File "/usr/local/lib/python3.8/dist-packages/skyfield/keplerlib.py", line 292, in eccentric_anomaly
    raise ValueError('eccentric anomaly failed to converge')
ValueError: eccentric anomaly failed to converge

I switched Skyfield versions back and forth by first doing an uninstall and then installing a specific version using sudo pip3 install -Iv skyfield==1.44 and similarly for 1.43.1. I also tested on other files from MPC and the main database from Lowell and got a similar result.

I searched for eccentric in case I've missed something (screen out dodgy data) and all I could find was one related to comets and another related to minor planets mentioning eccentricity.

Getting late, so apologies if I've blundered somewhere...

[brandon-rhodes]

Indeed, it looks like you uncovered a problem with this newer and faster routine! It only works when the mean anomaly is positive. When, instead, $-\pi ≤ M &lt; 0$, it tosses a negative number into a square root and, depending on its mood, either throws an immediate exception or else returns nan.

I was tipped off to the problem by a quick read of @rmathar’s paper at

https://arxiv.org/pdf/2108.03215.pdf

which says “To simplify the notation, we discuss only the cases where M ≥ 0” and suggests flipping the sign of both $M$ and the output $E$ if $M$ is negative.

Here's a fix that I've just committed:

c543eac

Can you try out the development version of Skyfield and seeing if the fix works for your case? Thanks!

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

[Bernmeister]

Fixed; thank you!

[brandon-rhodes]

Version 1.45 with this fix has just been released. Thanks for double-checking me, @Bernmeister!