vallado.lambert fails for long way transfers

[paul-witsberger]

When trying to set up poliastro.iod.vallado.lambert, I receive the following error when I set short=False (long path transfer):

Traceback (most recent call last):
  File "C:/Users/pawit/.PyCharmCE2019.3/config/scratches/scratch_4.py", line 11, in <module>
    v = lambert_v(k, r0, rf, tof, short=short, numiter=numiter, rtol=rtol).__next__()
  File "C:\Users\pawit\Documents\research_env\lib\site-packages\poliastro\iod\vallado.py", line 48, in lambert
    v0, v = vallado_fast(k_, r0_, r_, tof_, short, numiter, rtol)
  File "C:\Users\pawit\Documents\research_env\lib\site-packages\poliastro\core\iod.py", line 151, in vallado
    raise RuntimeError("Maximum number of iterations reached")
RuntimeError: Maximum number of iterations reached

Here is my minimal (non)working example:

from astropy import units as u, constants as c
from poliastro.iod.vallado import lambert as lambert_v

k = c.GM_earth.to(u.km ** 3/ u.s ** 2)
r0 = [10000., 0, 0] * u.km
rf = [8000., -5000, 0] * u.km
tof = (2 * u.hour).to(u.s)
short = False
numiter = 1000
rtol = 1e-6
v = lambert_v(k, r0, rf, tof, short=short, numiter=numiter, rtol=rtol).__next__()

print('v1 = ')
print(v[0])
print('v2 = ')
print(v[1])

It has worked for all of the tests I gave it when short=True, and failed every time short=False.

[jonpsy]

I'm a newbie, but may I take this issue?.

[jorgepiloto]

Hi @paul-witsberger, I have been able to reproduce your issue. Let me try to investigate a little more 😄 By the way, @nanusai feel free to take this issue 👍

[jonpsy]

After much digging the internet and talks with @Juanlu001 and @jorgepiloto , I've drawn the following conclusions from the problem.

1. I tried to compare izzo and vallado(short = True) and apparently they both yield different answers,

from poliastro.iod.vallado import lambert as lambert_v
from poliastro.iod.izzo import lambert as lambert_i 
from astropy import units as u, constants as c


if __name__ == "__main__":
    k = (c.GM_earth.to(u.km ** 3/ u.s ** 2))
    r0 = ([10000., 0, 0] * u.km)
    rf = ([8000., -5000, 0] * u.km)
    tof = 2*u.hour
    v_vallado = lambert_v(k, r0, rf, tof,short=True) #This gets the correct anser
    v_iz = lambert_i(k,r0,rf,tof,0)
    print(next(v_vallado))
    print(next(v_iz))

Which yields

(<Quantity [ 5.85488104, -1.74457682,  0.        ] km / s>, <Quantity [-6.2545192 ,  1.72835347, -0.        ] km / s>)
(<Quantity [0.36591277, 5.8228806 , 0.        ] km / s>, <Quantity [3.99397599, 4.78236576, 0.        ] km / s>)

After some search, I found this, another algorithm which solves using Izzo and Lancaster, E.R.,Blanchad method. The solution provided by this algorithm matches with v_vallado(short=True) but fails to match with v_izzo.

2. I've noticed in the case of vallado, for a given r0 and rf, there seems to be a fixed upper and lower bound under which the user can input tof .I reckon this is the paradigm the author is facing, for example we'll feed the parameters provided by the author to lambert_v(short=False) .

from poliastro.iod.vallado import lambert as vallado
from astropy import units as u, constants as c


if __name__ == "__main__":
    k = (c.GM_earth.to(u.km ** 3/ u.s ** 2))
    r0 = [10000., 0, 0] * u.km
    rf = [8000., -5000, 0] * u.km
    tof = 2 * u.hour
    numiter = 10**3
    rtol = 1e-6
    print("\nLong\n")
    v_long = vallado(k, r0, rf, tof, short=False, numiter=numiter, rtol=rtol) 
    print(next(v_long))

This would cause "RuntimeError: Maximum number of iterations reached" as the author pointed out. But if we comment out Maximum iterations reached at iod.py, it will return velocity vectors corresponding to tof_max which is

<Quantity [-5.22675825, 1.90842479, -0. ] km / s>, <Quantity [ 5.84298844, -1.26633679, 0. ] km / s>)

The above-mentioned vectors corresponds to tof_input :1989.6321243010739 seconds which stops changing after itr 55 , so clearly tof_input(7200.0 seconds or 2*u.hour) > tof_max(1989.6321243010739) and hence the error , and hence the proposed solution.

I've tried my best to explain the scenario but I'll be glad to answer your queries ! .

[jorgepiloto]

Origin of this algorithm and possible errata ⚠️

This routine is based on "Algorithm 58" from "Fundamentals of Astrodynamics" book. After comparing poliastro's implementation of the algorithm and the one stated in the book, I noticed there is a difference in the values for psi_up and psi_low:

• poliastro: psi_up = 4.0 * pi and psi_low = -4.0 * pi
• vallado58: psi_up = 4.0 * pi ** 2 and psi_low = -4.0 * pi ** 2 (refer to the algorithm script since it seems to be an errata in the book between de values)

Proposed solution

By correcting the actual values for psiof our vallado implementation, it was possible to run the original issue code without any problem. Tests were also passing, thought an increase of iterations was required in one of them.

Regarding the pull request #845, if a numiter argument is passed, the exception for reached maximum number of iterations should be kept. Code modifications do not affect code workflow but force it to return the last value. I will let a review in the pull request.

[astrojuanlu]

Thanks @nanusai for the investigation and congratulations @jorgepiloto for finding the root cause! Let's then have a pull request that includes:

• A test that checks that both implementations give the same results, within some accuracy?
• A test that checks that long way transfers work?
• The fixes in psi_up and psi_low that make the above pass

[jonpsy]

@Juanlu001 I've increased the num iter for test_curtis53 as @jorgepiloto requested.
I've also

• Fixed psi_up and psi_low.
• Modified test_issue840 which now asserts vallado(short=False) and izzo equality with rtol =1e-6
  I've sent a pull request, check it out 😄

[jonpsy]

After much digging the internet and talks with @Juanlu001 and @jorgepiloto , I've drawn the following conclusions from the problem.

1. I tried to compare izzo and vallado(short = True) and apparently they both yield different answers,

from poliastro.iod.vallado import lambert as lambert_v
from poliastro.iod.izzo import lambert as lambert_i 
from astropy import units as u, constants as c


if __name__ == "__main__":
    k = (c.GM_earth.to(u.km ** 3/ u.s ** 2))
    r0 = ([10000., 0, 0] * u.km)
    rf = ([8000., -5000, 0] * u.km)
    tof = 2*u.hour
    v_vallado = lambert_v(k, r0, rf, tof,short=True) #This gets the correct anser
    v_iz = lambert_i(k,r0,rf,tof,0)
    print(next(v_vallado))
    print(next(v_iz))

Which yields

(<Quantity [ 5.85488104, -1.74457682,  0.        ] km / s>, <Quantity [-6.2545192 ,  1.72835347, -0.        ] km / s>)
(<Quantity [0.36591277, 5.8228806 , 0.        ] km / s>, <Quantity [3.99397599, 4.78236576, 0.        ] km / s>)

After some search, I found this, another algorithm which solves using Izzo and Lancaster, E.R.,Blanchad method. The solution provided by this algorithm matches with v_vallado(short=True) but fails to match with v_izzo.

1. I've noticed in the case of vallado, for a given r0 and rf, there seems to be a fixed upper and lower bound under which the user can input tof .I reckon this is the paradigm the author is facing, for example we'll feed the parameters provided by the author to lambert_v(short=False) .

from poliastro.iod.vallado import lambert as vallado
from astropy import units as u, constants as c


if __name__ == "__main__":
    k = (c.GM_earth.to(u.km ** 3/ u.s ** 2))
    r0 = [10000., 0, 0] * u.km
    rf = [8000., -5000, 0] * u.km
    tof = 2 * u.hour
    numiter = 10**3
    rtol = 1e-6
    print("\nLong\n")
    v_long = vallado(k, r0, rf, tof, short=False, numiter=numiter, rtol=rtol) 
    print(next(v_long))

This would cause "RuntimeError: Maximum number of iterations reached" as the author pointed out. But if we comment out Maximum iterations reached at iod.py, it will return velocity vectors corresponding to tof_max which is

<Quantity [-5.22675825, 1.90842479, -0. ] km / s>, <Quantity [ 5.84298844, -1.26633679, 0. ] km / s>)

The above-mentioned vectors corresponds to tof_input :1989.6321243010739 seconds which stops changing after itr 55 , so clearly tof_input(7200.0 seconds or 2*u.hour) > tof_max(1989.6321243010739) and hence the error , and hence the proposed solution.

I've tried my best to explain the scenario but I'll be glad to answer your queries ! .

Before this thread is closed I wanted to clear the misconception I might've created with the above post.

1. The izzo algorithm is indeed working perfectly fine, it's value is matching with vallado(short=False).
2. This hypothesis is flawed and doesn't hold its ground when the algorithm is fixed. Although putting in unrealistically high/low values as tof would cause "Maximum iterations reached" as is what's expected of it.

The hypothesis was just a logical explanation of all the things I observed which proved to be incorrect. Again kudos to @jorgepiloto for actually solving the problem and @Juanlu001 for guiding me towards my first merged PR.

This thread maybe closed now.

[astrojuanlu]

Fixed in #845! @paul-witsberger if you want to test with master you're welcome to do so :) And if the issue was still not fixed, please comment here.