Poliastro: Examine solutions in Lambert notebook

[astrojuanlu]

🐞 Problem

With the current master version, I'm getting solutions in strange locations of the plot:

I am not sure this was always the case, so it's worth it to rerun the notebook and assess whether these results are correct.

📋 Steps to solve the problem

• Comment below about what you've started working on.
• Add, commit, push your changes
• Submit a pull request and add this in comments - Addresses #<put issue number here>
• Ask for a review in comments section of pull request
• Celebrate your contribution to this project 🎉

[astrojuanlu]

This has been like that at least for two years https://docs.poliastro.space/en/v0.7.0/examples/Revisiting%20Lambert's%20problem%20in%20Python.html perhaps it was never quite right?

[astrojuanlu]

The plot is taken from the Izzo paper:

https://arxiv.org/abs/1403.2705

[jorgepiloto]

I also saw that the labels Elliptic and Hyperbolic are switched, I was uploading this image but you were faster than me :smile:

[jorgepiloto]

The current state of this:

1. After comparing the poliastro Izzo's algorithm against pykep it got exactly the same results, so algorithm implementation seems to be right.

2. I focused my attention on the _compute_T_min private function and realized that wrong results were obtained only for those values of ll = -1. If we increase the y-axis limit to also see the solutions M=3 we can notice that red crosses corresponding to ll = -1 are the problematic points:

1. This value of lambda is a problematic one as this comment states. Therefore, I will try to further inspect this code section :+1:

[jorgepiloto]

Alright, so it seems that the initial condition x_i when ll = -1 tremendously affects the output of the function _compute_T_min. By forcing each initial condition x_i to each M, those red crosses matched the lines. See:

However, this is not a solution, of course. A deep study on _tof_equation and _halley method is required to better see why the initial condition makes the solution to diverge if ll = -1 :confused:

[astrojuanlu]

From the paper:

Values of λ^2 close to unity indicate a chord of zero length, a case which is indeed extremely interesting in interplanetary trajectory design as it is linked to the design of resonant transfers.

[astrojuanlu]

Encouraged by @jorgepiloto progress on this issue, I couldn't help but look at it... And I found the problem!

@@ -370,7 +378,8 @@ def _compute_T_min(ll, M, numiter, rtol):
             y = _compute_y(x_i, ll)
             T_i = _tof_equation(x_i, y, 0.0, ll, M)
             x_T_min = _halley(0.1, T_i, ll, rtol, numiter)
-            T_min = _tof_equation(x_T_min, y, 0.0, ll, M)
+            new_y = _compute_y(x_T_min, ll)
+            T_min = _tof_equation(x_T_min, new_y, 0.0, ll, M)

     return [x_T_min, T_min]

The value of y is being reused with an old value of x for the case ll != 1, M != 0. For some reason, this seems to affect most the ll == -1 case.

The reason why the value y is computed outside of the time of flight equation is to avoid a repetitive calculation when evaluating the Householder iteration scheme:

https://github.com/poliastro/poliastro/blob/c9afdc9b5fcbd666f8da3d0c231842a13f72d3bb/src/poliastro/core/iod.py#L448-L453

Which was introduced in https://github.com/poliastro/poliastro/commit/b2d1fce07af5ebd2f6162290f717436949627e11 in a buggy way, therefore breaking the algorithm.

perhaps it was never quite right?

I knew it was right _at some point_!

https://github.com/poliastro/poliastro/commit/986b85d151bcd22a09e19a4fd74cab610375e7f6

https://github.com/poliastro/poliastro/blob/986b85d/examples/Revisiting%20Lambert's%20problem%20in%20Python.ipynb