`farnocchia` propagator returns `NaN` for some orbits

[Yash-10]

I tried to create ephemerides for an orbit that had a very small period (9.3605811e-8s).

import numpy as np; from poliastro.twobody.orbit import Orbit; from poliastro.twobody.propagation import cowell; from poliastro.ephem import Ephem

import astropy.units as u; from poliastro.bodies import Earth; from astropy.time import Time

r = np.array([-500, 1500, 4012.09]) << u.km
v = np.array([5021.38, -2900.7, 1000.354]) << u.km / u.s
orbit = Orbit.from_vectors(Earth, r, v, epoch=Time("2020-01-01"))

tofs = [100, 125, 150] << u.s
# Propagate the secondary body to generate its position coordinates.
rr, vv = cowell(
    Earth.k,
    orbit.r,
    orbit.v,
    tofs,
)

ephem = Ephem.from_orbit(orbit, orbit.epoch+tofs)
print(ephem.rv())

Output

(<Quantity [[nan, nan, nan],
           [nan, nan, nan],
           [nan, nan, nan]] km>, <Quantity [[nan, nan, nan],
           [nan, nan, nan],
           [nan, nan, nan]] km / s>)

For any t < 100 * u.s , the values were non-nan, whereas they start to become nan from the 100s mark.

I initially thought this might be because the period of the orbit is small whereas tof is high. I tried with another orbit, again with a small period of (8.3325481e-5s) but I got a non-nan result:

r = np.array([-5000, 1500, 4012.09]) << u.km
v = np.array([502.38, -2900.7, 1000.354]) << u.km / u.s
orbit = Orbit.from_vectors(Earth, r, v, epoch=Time("2020-01-01"))

tofs = [100, 125, 150] << u.s
rr, vv = cowell(
    Earth.k,
    orbit.r,
    orbit.v,
    tofs,
)

ephem = Ephem.from_orbit(orbit, orbit.epoch+tofs)
print(ephem.rv())

gives

(<Quantity [[  45239.34869765, -288568.58094236,  104045.47706077],
           [  57799.18822177, -361085.69078418,  129053.80351732],
           [  70359.02749432, -433602.79905165,  154062.12941078]] km>, <Quantity [[  502.39358751, -2900.68443511,  1000.33307316],
           [  502.39357529, -2900.68435797,  1000.33304546],
           [  502.393567  , -2900.68430656,  1000.33302714]] km / s>)

Also, for an orbit with a period of 1.0245247e-7s, the output values are nan.

Any suggestions on this would be helpful. Thanks!

[applysomesky]

Hi Yash-10,

I noticed a few things:

1. I quickly looked at your inputs for your first run, and it seems the total distance of r is smaller (at 4312 km) than the radius of the Earth (which is about 6378km at the equator), which for your second orbit is not the case. Not sure if that is the cause for the error, I am personally not sure what happens if poliastro detects a radius lower than the main attractor's radius. To rule out this maybe you can use the moon as the main attractor for the first case.
2. I also noticed your orbital periods, after checking your semi major axis, it is negative, which means its a hyperbolic orbit, in other words it does not have an orbital period since it only passes the main object once. I am not sure how you calculated the orbital periods. (I am also not sure what kind of thing you are trying to simulate, but 6000 km/s is very very fast).

Note: I haven't used poliastro in a while, so I am not an expert on the package.

I hope this helps a bit.

[astrojuanlu]

Thanks @WolfsSky for chiming in! Indeed, both orbits are hyperbolic.

In any case, the problem seems to be that, while cowell returns an answer, the default propagator fails beyond a certain point:

In [25]: orbit
Out[25]: 4285 x -4285 km x 104.0 deg (GCRS) orbit around Earth (♁) at epoch 2020-01-01 00:00:00.000 (UTC)

In [26]: orbit.rv()
Out[26]: 
(<Quantity [-500.  , 1500.  , 4012.09] km>,
 <Quantity [ 5021.38 , -2900.7  ,  1000.354] km / s>)

In [27]: propagate(orbit, [73] << u.s)
Out[27]: 
<CartesianRepresentation (x, y, z) in km
    [(366059.84053624, -210250.93953118, 77036.53364381)]
 (has differentials w.r.t.: 's')>

In [28]: propagate(orbit, [74] << u.s)
Out[28]: 
<CartesianRepresentation (x, y, z) in km
    [(nan, nan, nan)]
 (has differentials w.r.t.: 's')>

In [29]: propagate(orbit, [74] << u.s, method=cowell)
Out[29]: 
<CartesianRepresentation (x, y, z) in km
    [(371081.20762409, -213151.63704118, 78036.8682027)]
 (has differentials w.r.t.: 's')>

Therefore, let's rule out Ephem from this issue.

[Yash-10]

Had some time so did a bit of debugging. All outputs are nan because of this line:

poliastro/src/poliastro/core/propagation/farnocchia.py

Line 284 in 069167a

F = M_to_F(M, ecc)

.

_newton_hyperbolic does not converge for propagate(orbit, [74] << u.s) (from the above example) due to which it returns nan:

poliastro/src/poliastro/core/angles.py

Line 38 in 069167a

return np.nan

.
Maybe the solution could be to raise a non-convergence warning. However, it might be also good to investigate the cause of this non-convergence.

[astrojuanlu]

Thanks a lot for the thorough debugging @Yash-10 ! Do you know if it does converge increasing the number of iterations? Or does it just diverges?

[Yash-10]

Interestingly, it does seem to converge (i.e. abs(p-p0) < tol) at the 101th iteration.

The default no. of iterations for _newton_hyperbolic in M_to_F is 100, and so the failure is because it takes just one iteration more!

poliastro/src/poliastro/core/angles.py

Line 327 in 069167a

F = _newton_hyperbolic(F0, args=(M, ecc), maxiter=100)

So maybe increasing maxiter might be a solution...

[astrojuanlu]

Interesting!

Well, the problem here I think is that there will always be an orbit that will not converge. Unless we set maxiter to float("inf") or something like that :) Otherwise, just increasing the number of iterations to 101 seems a bit silly, because we'd be able to find an orbit that would need 102, and so forth.

I see three options:

• Increase maxiter to some other finite number. Not a super solid solution, but it doesn't harm.
• Set maxiter to float("inf"). A bit dangerous, are we sure that the convergence is guaranteed?
• Allow the user to more easily set the maxiter in the farnocchia propagator.

Thoughts?

[Yash-10]

Well, I recall #1247. It turns out it's a mistake from my side. I reverted to an older commit before that PR, and now I do not get any nan values.

That approach was taken from Vallado (discussion in #1240). Maybe we could work on understanding it better?

[astrojuanlu]

Oh, I see. Well, I guess slightly different initial guesses will have slightly different numerical performance. If reverting #1247 fixes this issue, then let's do that and add a regression test, so we don't forget again.

Maybe Vallado didn't test his initial guess with super high eccentricity orbits like this one? And even if they seem rare, I wouldn't be surprised if orbits like this emerge in impact trajectories or things like that.