Propagator mean_motion hangs for some r, v vectors around Earth

[astrojuanlu]

Comes from #474, by @TimothySHamilton

The problem can be better seen by disabling Numba JIT, which shows some NumPy floating point warnings, and then turning them into errors:

$ NUMBA_DISABLE_JIT=1 ipython --no-banner

In [1]: import numpy as np

In [2]: np.seterr(all="raise")
Out[2]: {'divide': 'warn', 'over': 'warn', 'under': 'ignore', 'invalid': 'warn'}

In [3]: import numpy as np
   ...: import math
   ...: from astropy import units as u
   ...: from poliastro.bodies import Earth, Moon
   ...: from poliastro.twobody import Orbit
   ...: 
   ...: r=[8.e3, 1.e3, 0.]*u.km
   ...: v=[-0.5, -0.5, 0.]*u.km/u.s
   ...: orbit1=Orbit.from_vectors(Earth,r,v)
   ...: orbit2=orbit1.propagate(1.*u.h)
   ...: 
   ...: 
---------------------------------------------------------------------------
FloatingPointError                        Traceback (most recent call last)
<ipython-input-3-08d7b74965c9> in <module>()
      8 v=[-0.5, -0.5, 0.]*u.km/u.s
      9 orbit1=Orbit.from_vectors(Earth,r,v)
---> 10 orbit2=orbit1.propagate(1.*u.h)

~/.miniconda36/envs/poliastro37/lib/python3.7/site-packages/poliastro/twobody/orbit.py in propagate(self, value, method, rtol, **kwargs)
    403                 time_of_flight = time.TimeDelta(value)
    404 
--> 405             return propagate(self, time_of_flight, method=method, rtol=rtol, **kwargs)
    406 
    407     def sample(self, values=None, method=mean_motion):

~/.miniconda36/envs/poliastro37/lib/python3.7/site-packages/poliastro/twobody/propagation.py in propagate(orbit, time_of_flight, method, rtol, **kwargs)
    177 
    178     """
--> 179     r, v = method(orbit, time_of_flight.to(u.s).value, rtol=rtol, **kwargs)
    180     return orbit.from_vectors(orbit.attractor, r * u.km, v * u.km / u.s, orbit.epoch + time_of_flight, orbit.plane)

~/.miniconda36/envs/poliastro37/lib/python3.7/site-packages/poliastro/twobody/propagation.py in mean_motion(orbit, tofs, **kwargs)
    118 
    119     if not hasattr(tofs, '__len__'):
--> 120         return mean_motion_fast(k, r0, v0, tofs)
    121 
    122     results = [mean_motion_fast(k, r0, v0, tof) for tof in tofs]

~/.miniconda36/envs/poliastro37/lib/python3.7/site-packages/poliastro/core/propagation.py in mean_motion(k, r0, v0, tof)
     33 
     34     # get the initial mean anomaly
---> 35     M0 = nu_to_M(nu0, ecc)
     36     # strong elliptic or strong hyperbolic orbits
     37     if np.abs(ecc - 1.0) > 1e-2:

~/.miniconda36/envs/poliastro37/lib/python3.7/site-packages/poliastro/core/angles.py in nu_to_M(nu, ecc, delta)
    183     else:
    184         D = nu_to_D(nu)
--> 185         M = D_to_M(D, ecc)
    186     return M
    187 

~/.miniconda36/envs/poliastro37/lib/python3.7/site-packages/poliastro/core/angles.py in D_to_M(D, ecc)
    155 @jit
    156 def D_to_M(D, ecc):
--> 157     M = _kepler_equation_parabolic(D, 0.0, ecc)
    158     return M
    159 

~/.miniconda36/envs/poliastro37/lib/python3.7/site-packages/poliastro/core/angles.py in _kepler_equation_parabolic(D, M, ecc)
     26 @jit
     27 def _kepler_equation_parabolic(D, M, ecc):
---> 28     return M_parabolic(ecc, D) - M
     29 
     30 

~/.miniconda36/envs/poliastro37/lib/python3.7/site-packages/poliastro/core/angles.py in M_parabolic(ecc, D, tolerance)
     41     k = 0
     42     while not small_term:
---> 43         term = (ecc - 1.0 / (2.0 * k + 3.0)) * (x ** k)
     44         small_term = np.abs(term) < tolerance
     45         S += term

FloatingPointError: overflow encountered in double_scalars

And then, mean_motion fails for true anomalies very close to pi in the near parabolic case:

In [14]: M_parabolic(ecc, nu_to_D(np.pi + 0.01))  # overflow

In [15]: M_parabolic(ecc, nu_to_D(np.pi - 0.01))  # overflow

This deserves closer inspection. @nikita-astronaut was the original author that read the corresponding paper, but I don't think he will have the time to answer. Perhaps one solution would be to use the elliptic solution when the true anomaly is close to pi even if the eccentricity is close to one:

poliastro/src/poliastro/core/angles.py

Lines 176 to 186 in 1f12764

	def nu_to_M(nu, ecc, delta=1e-2):
	if ecc > 1 + delta:
	F = nu_to_F(nu, ecc)
	M = F_to_M(F, ecc)
	elif ecc < 1 - delta:
	E = nu_to_E(nu, ecc)
	M = E_to_M(E, ecc)
	else:
	D = nu_to_D(nu)
	M = D_to_M(D, ecc)
	return M

We should also find out if we could prevent at least the hangs by compiling some functions with error_model="python", as advised in numba/numba#1256 (comment).

[astrojuanlu]

While trying to answer #494, I discovered that most orbits coming from the Lambert problem in my particular case were difficult to propagate, because they were almost parabolic. Adding a or nu < 3 condition was too naïve and led to clearly wrong results, also because there's a leap:

>>> nu_to_M(t1.nu, t1.ecc, delta=1e-2)  # Strange!
-667.3977108062284 rad
>>> nu_to_M(t1.nu, t1.ecc, delta=1e-2) % (2 * np.pi * u.rad)
4.903117061987366 rad
>>> nu_to_M(t1.nu, t1.ecc, delta=1e-3) % (2 * np.pi * u.rad)
5.979356985157597 rad

So there's some sort of discontinuity or instability that should be investigated.

Moreover, the nu_to_M and M_to_nu exhibit some weird behavior very localized in the vicinity of nu = pi, 0.985 < ecc < 1.0 (for delta=1e-2), as can be seen in this plot:

Here I'm plotting the difference between nu and the roundtrip conversion, M_to_nu_fast(nu_to_M_fast(nu, ecc, delta), ecc, delta). The black region is where the algorithm did not even converge.

Code to generate:

This code is quite slow to run, because the JIT has to be disabled (otherwise the propagator just hangs) and it turns out that catching the exception is slow too. With more time, a plot with better resolution could be generated. Also, it would be interesting to study the influence of `delta`.

import numpy as np
from poliastro.core.angles import M_to_nu as M_to_nu_fast, nu_to_M as nu_to_M_fast

np.seterr(all="raise")

N = 501

eccs = np.linspace(0.985, 1, num=N)
nus = np.linspace(0, 6, num=N)
delta = 1e-2

res = np.zeros((N, N))

for ii, ecc in enumerate(eccs):
    for jj, nu in enumerate(nus):
        try:
            nu_final = M_to_nu_fast(nu_to_M_fast(nu, ecc, delta), ecc, delta)
        except FloatingPointError:
            nu_final = 20  # This signals that something wrong happened

        res[jj, ii] = abs(nu - nu_final)

res[res == 0.0] = np.nan  # To have a white pixel in the plot
res[res == 2 * np.pi] = np.nan  # To have a white pixel in the plot

[jorgepiloto]

I am going to start working on this. Those extreme cases where numerical methods struggle to obtain the most accurate result are really interesting. Hope #577 be ready for today.

[jorgepiloto]

Some days ago I found in Battin 1999 an algorithm developed by Gauss that is said to work really well for this extreme eccentricity conditions. After found a paper where the algorithm is explained, coefficients are detailed and it even includes some test at the end.

I will study more in depth this: https://iopscience.iop.org/article/10.1086/128750/pdf

[jorgepiloto]

After successfully implementing on my local machine Gauss algorithm for near parabolic orbits, problem with extreme eccentricities was solved. Algorithm converged really fast, even for orbits with ecc=0.9999 and no hangout of Jupyter kernel was received.

A test using Battin 1999 example was done in order to check that the algorithm was working properly. However, when applying it to #474, I obtain an error related with a negative quantity inside an sqrt in the algorithm. Still working on this... 🚀

[astrojuanlu]

Our hero! 😍 Looking forward to that PR!

…

On Sun, 24 Mar 2019, 11:51 Jorge Martinez, ***@***.***> wrote: After successfully implementing on my local machine Gauss algorithm for near parabolic orbits, problem with extreme eccentricities was solved. Algorithm converged really fast, even for orbits with ecc=0.9999 and no hangout of Jupyter kernel was received. A test using Battin 1999 example was done in order to check that the algorithm was working properly. However, when applying it to #474 <#474>, I obtain an error related with a negative quantity inside an sqrt in the algorithm. Still working on this... 🚀 — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#475 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AATUZUjVRpgGcGvPrs8jyGFr2h4t3RkBks5vZ1ifgaJpZM4XzZJq> .