File Exchange

## Robust solver for Lambert's orbital-boundary value problem

version 1.3.0.0 (36.1 KB) by Rody Oldenhuis

### Rody Oldenhuis (view profile)

Solves any Lambert-problem robustly. Can be compiled to increase efficiency.

Updated 19 Jun 2018

A Lambert-orbital boundary value problem can be stated as
"Find the orbit/trajectory of a spacecraft that flies from position [r1] to [r2], taking a time [tf] and making [m] complete orbits before arriving at [r2]. "
The solution to each Lambert-problem is NOT unique; one can travel to [r2] via the long way or the short way, and for [m > 0] there are almost always two ellipses that satisfy the boundary conditions, so that [m > 0] has four distinct solutions.
This function solves any Lambert problem *robustly*. It can handle short-way solutions (the default), long way solutions (by passing negative [tf]), or left-branch (default) or right-branch (by passing negative [m]) solutions in case [m > 0]. It uses two separate solvers; the first one it tries is a new and unpublished algorithm developed by Dr. D. Izzo from the European Space Agency [1]. This version is extremely fast, but especially for larger [m] it still fails quite frequently. In such cases, a MUCH more robust algorithm is started (the one by Lancaster & Blancard [2], with modifcations, initial values and other improvements by R.Gooding [3]), which is a lot slower partly because of its robustness.

This routine can be compiled to increase its speed by a factor of 20-50, which is certainly advisable when an application using this function requires a great number of Lambert problems to be solved. The entire routine is written in embedded MATLAB, so it can be compiled with the emlmex() function. It is described in the function's comments how to accomplish this.

Tested on WinXP/32 and Ubuntu 9.10/32. Tested on MATLAB 2008a through 2009b. As compilation is always a troublesome affair, please let me know what sort of problems you encounter so I can try to improve the code for that.

References:
[1] Izzo, D. ESA Advanced Concepts team. Code used available in MGA.M, on http://www.esa.int/gsp/ACT/inf/op/globopt.htm. Last retreived Nov, 2009.
[2] Lancaster, E.R. and Blanchard, R.C. "A unified form of Lambert's theorem." NASA technical note TN D-5368,1969.
[3] Gooding, R.H. "A procedure for the solution of Lambert's orbital boundary-value problem. Celestial Mechanics and Dynamical Astronomy, 48:145–165,1990.

If you find this work useful, please consider a donation:
https://www.paypal.com/cgi-bin/webscr?cmd=_s-xclick&hosted_button_id=6G3S5UYM7HJ3N

Vasco Grilo

### Vasco Grilo (view profile)

Thank tou very much for the solver! Just one note:

The solver is not blinded against erroneous input. The following situations were detected:
\begin{itemize}
\item For equal initial and final positions, initial or final position at the central body, null gravitational parameter, or null time of flight, the initial and final velocity components are NaN\;(not a number).
\item For negative gravitational parameter, the initial and final velocity components have an imaginary part.
\item For negative time of flight, the initial and final velocity components are real.
\end{itemize}
For all these cases, the exit flag is $1$\;, which is supposed to indicate success, instead of $-1$\;, which indicates that the given problem has no solution. Therefore the outputted results could be interpreted as correct, and undermine the following calculations.

Christopher Simpson

Hao Mei

### Hao Mei (view profile)

Great work! Thanks a lot.

BlueEyes

### BlueEyes (view profile)

Hi all.
May I provide a working example?
-----------
>> format long g
>> AU= 149597870; % astronomical unit in km
>> r1vec= [-.723554659669305e0*AU, 0.669919796762551e0*AU, -0.125124314151737e-4*AU];
>> r2vec= [0.643189121761322e0*AU, -0.338282672522157e0*AU, -0.417522705944112e-1*AU];
>> tf= 141.0; % days
>> m= 0.0; % complete nr. of orbits
>> muC= 1.327124e11; % Sun grav. parameter in km^3/s^2
>> %
>> [V1, V2, extremal_distances] = lambert(r1vec, r2vec, tf, m, muC)
------ OUTPUT ------
V1 =

17.5019752271602 20.4482145730657 -5.95060968592276

V2 =

-18.6865568988041 -31.404421393359 7.90749775832099

extremal_distances =

108895604.041308 148086264.397052
---------
Cheers

Rody Oldenhuis

### Rody Oldenhuis (view profile)

@Bielorusse: I've never seen this behavior...can you provide an example?

kenji kitamura

### kenji kitamura (view profile)

Very useful. Thank you for your great work!!

Bielorusse

### Bielorusse (view profile)

Hello,

I am trying to use exclusively the Lancaster Blanchard function, but I obtain peculiar results, with a 10 000 factor for the value of x, whereas Izzo's function's results are correct with the same input parameters.

Has anyone had similar issues? Can the "lambert_LancasterBlanchard" function be used alone without the "lambert" function?

Thanks!

Brian Holm-Hansen

Rody Oldenhuis

### Rody Oldenhuis (view profile)

@DavidSage: Corrected, thanks!

David Sage

### David Sage (view profile)

Found the problem:

Line 664 should be "elseif" instead of "else if".

David Sage

### David Sage (view profile)

Trying this on Matlab 2015a, I get the following error:

"Error: File: lambert.m Line: 783 Column: 4
The function "minmax_distances" was closed with an 'end', but at least one other function definition was not. To
avoid confusion when using nested functions, it is illegal to use both conventions in the same file."

I'm a total newbie, but can anyone suggest how I can get by this and just start to get real output?

Rody Oldenhuis

### Rody Oldenhuis (view profile)

@Simon: you are absolutely correct; I've patched it with your proposed fix. Thanks!

Simon Tardivel

### Simon Tardivel (view profile)

@Rody

I think there's an unexpected behavior that prevents the function from computing some solvable situations.

As you tested the function in previous Matlab functions, it may have been only present in recent matlab versions. The culprit is in the robust solver, around line 280, in the calculation of d. Right now, the function outputs a NaN if E<0 and f+g<0. Indeed, it computes d as "something + 0*-Inf" which yields NaN.

I'm not sure this is desired behavior. I propose the following fix (but other syntaxes are possible):
if E<0
d = atan2(f, g) + pi*m;
else
if E == 0
d = 0;
else
d = log(max(0, f + g));
end
end

Rody Oldenhuis

### Rody Oldenhuis (view profile)

@Simon I updated the code with your suggestions. Thanks!

Simon Tardivel

### Simon Tardivel (view profile)

@Rody Alright!

So indeed the new function is codegen. I haven't looked too much into it, but it seems to be the central platform of translating matlab code into "other" stuff and notably their sped-up version, the mex file.

The syntax is extremely similar, as you pointed out. Using your example_input, simply type:
codegen lambert -args example_input

This outputs a C compiled, optimized lambert_mex.mex***. The *** depends on your platform, it gives me w64 as i'm on Windows 64. Of course, this function can be renamed lambert.mexw64, placed in the same folder as the other lambert.m it will be called over the .m lambert.

Now, on the specific example I run, it speeds up the computation of lambert by a factor 9, not 20-50. On my computer, running matlab's profiler, the mex gives me 62 000 lamberts per second with mex, compared to 6700 lamberts per second with the .m. The speed may also be problem dependent, I have just have my own example.

That may be a sign that matlab is actually better now than before at running its own code haha! Still, 12h or 110h makes very big difference! Maybe one should play with the options of codegen (e.g. -O for the optimization options) to try and make that faster, but right now it is very sufficient and so easy!

Thanks for the help, and again thanks for the code!

Rody Oldenhuis

### Rody Oldenhuis (view profile)

@Simon: Indeed, I wrote this on a much older MATLAB version. In the meantime the MathWorks switched to a different compilation tool. I know they have issued warnings about using emlmex being deprecated, but apparently, they finally dropped support for it in R2014b.

I'm also a bit lost (I don't have anything newer than R2010a), but I know the google keyword is 'codegen' or 'coder'. Supposedly, from what I found, there is a one-to-one correspondence to emlmex syntax and codegen syntax (with codegen supporting more), but I have never been able to play with this in order to verify this claim.

If you do figure it out, please let me know the precise steps you took, then I can update the documentation here :)

Simon Tardivel

### Simon Tardivel (view profile)

First of all, thanks for this great piece of code! It is really useful.

I have a question regarding the compilation. I've changed computers (mac to windows) and Matlab distributions (i'm now on matlab r2014b). Previously I had compiled the function and it was very fast. But when I try compiling it, now it can't find the function emlmex.

I'm rather new to this compilation of matlab code. Most of the research i've done on internet seems to refer to obsolete functions/libraries -- or I may be using them wrong. The only actual compiling library i have found is the fixed-point toolbox. And the toolbox complains about the trigonometric functions, NAN and Inf, so I cannot use it.

Do you know what I am doing wrong here? I would really like to use the compiled version, as it is much faster!

Thanks!

Abel

### Abel (view profile)

No question ... the best lamber implementation. Very well done.

Rody Oldenhuis

### Rody Oldenhuis (view profile)

@Martin: do you happen to have some reliable comparisons between Battin's and Izzo's methods regarding their execution speed, both for 0-rev and m-rev cases?

I've happen to have been in contact recently with Dr. Izzo regarding an update he's done:

http://sourceforge.net/p/keptoolbox/code/ci/development/tree/src/lambert_problem.cpp

He claims he has further improved the routine's speed, and measured a 1.8× improvement. I'm real short on time these days, so if you can find the time, can you verify this statement? And/or compare this new version to Battin's method, both regarding speed and reliability?

Rody Oldenhuis

### Rody Oldenhuis (view profile)

@Martin: Thank you for the kind words :) I have just uploaded an update with the correction. I haven't cross-checked it with Gooding's paper, but I trust your statement ^_^

Martin

### Martin (view profile)

I have tested this lambert targeter against MANY others (Universal variables, Battin's method, Izzo's alone, Lancaster & Blancard alone, etc.) and it is the most robust and almost as fast as Battin's method in zer0-rev cases. I used this code with great success in my Masters and PhD work. Although in my testing of ALL cases (shortway/longway, leftbranch, rightbranch, etc.) I found an typo/error on line 515:

W = x00 + 1.7*sqrt(2*(1 - phr));

which should actually be:

W = x00 + 1.7*sqrt(2 - phr/pi);

as agrees with this interpretation of the Gooding paper.

This is the only error I have found with EXTENSIVE testing and it only pertains to choosing rightbranch for multi-rev solutions.

Thanks for the great tool!

Nathan

### Nathan (view profile)

I have used this Lambert's solver extensively for my undergraduate spacecraft design class. Very fast for MATLAB, and it is quite robust, too.

Hermes

### Hermes (view profile)

I have a problem using this program. First, I created an example_input and ran the command
emlmex -eg example_input lambert.m

Now I have an input for my problem which is:

input = {[22592.145603,-1599.915239,-19783.950506],[1922.067697,4054.157051,-8925.727465],0.417,0,324859};

My question is how I should run this input and how I can access the output of the m-file.

Thank you,