I've been reading up on using the p-iteration method to solve Lambert's problem for choosing the correct interplanetary trajectory between two planets given the initial position of the source planet, final position of the target planet as well as a time of flight of the spacecraft. I can solve the equations manually, but what I wanted was to make a while loop in Mathematica that would automatically find the required time of flight. As an example, here is what I have so far (which I got from http://www.braeunig.us/space/interpl.htm):

G = 6.672*10^-11;
M = AstronomicalData["Sun", "Mass"];
m = AstronomicalData["Earth", "Mass"];
deltat = 207*86400; (*seconds*)
r[1] = {0.473265, -0.899215, 0};
r[2] = {0.066842, 1.561256, 0.030948 };

Subscript[μ, s] = G M/(149.597870*10^9)^3 (*In AU^3/s^2*)

(*p-iteration constants*)
r1 = Sqrt[r[1].r[1]]
r2 = Sqrt[r[2].r[2]]
deltaA = ArcCos[r[1].r[2]/(r1 r2)]/Degree
k = r1 r2 (1 - Cos[deltaA Degree])
l = r1 + r2
m = r1 r2 (1 + Cos[deltaA Degree])

(*p-iteration calculations*)
a = (m k p)/((2 m - l^2) p^2 + 2 k l p - k^2)
f = 1 - r2/p (1 - Cos[deltaA Degree])
g = (r1 r2 Sin[deltaA Degree])/Sqrt[Subscript[μ, s] p]
fdot = Sqrt[Subscript[μ, s]/p]
Tan[(deltaA Degree)/2] ((1 - Cos[deltaA Degree])/p - 1/r1 - 1/r2)
gdot = 1 - r1/p (1 - Cos[deltaA Degree])
deltaE = ArcCos[1 - r1/a(1 - f)] (*In radians*) (*used if a is positive*)
deltaF = ArcCosh[1 - r1/a(1 - f)] (*In radians*) (*used if a is negative*)
t = g + Sqrt[a^3/Subscript[μ, s]] (deltaE - Sin[deltaE ]) (*used if a is positive*)
t = g + Sqrt[(−a)3/Subscript[μ,s]](Sinh[deltaF]−deltaF) (*used if a is negative*)

I've been trying to make a while loop to perform iterations on the above calculations that will choose a value for p in order to match t as close as possible to deltat but have so far been unable to get a working solution. The way to choose a new value for p for each new iteration is done using the following equation:


where t[n-1] and t[n] are the first and second iterations of the above algorithm respectively.

Does anyone know how I could do something like this in a while loop? Any help would be great appreciated.

EDIT: With regards to using the FindRoot function, this is what I tried (I essentially replaced g, a and deltaE in

 t = g + Sqrt[a^3/Subscript[μ, s]] (deltaE - Sin[deltaE ]) 

with the basic equations I started with) which unfortunately brought up a few errors:

enter image description here

  • $\begingroup$ Apologies for the late reply, work has been a bit hectic. Yup, a can be negative in which case I would use $$t = g + Sqrt[(-a)^3/Subscript[\[Mu], s]](Sinh[deltaF] - deltaF)$$ where $$Cosh[deltaF]=1-(r1/a)(1-f)$$ $\endgroup$ Nov 26, 2013 at 15:17
  • $\begingroup$ It looks like the secant method of root-finding, which FindRoot can do automatically. You can google it and find many, many different lecture notes on implementing it. $\endgroup$
    – Michael E2
    Nov 27, 2013 at 16:21
  • $\begingroup$ I've been trying to use FindRoot to no avail and unfortunately the documentation doesn't have examples that help me too much. How would I used FindRoot in the context of the above equations? $\endgroup$ Nov 28, 2013 at 11:54

1 Answer 1


Solving your linked exercise in the below comment's link, assuming a >0 :

(* Problem Data *)
G = 6.672*10^-11;
deltat = 207*86400;(*seconds*)
M = AstronomicalData["Sun", "Mass"];
m = AstronomicalData["Earth", "Mass"];
r[1] = {0.473265, -0.899215, 0};
r[2] = {0.066842, 1.561256, 0.030948};
tf = 207 86400; (*207 days*)
muS = G M/(149.597870*10^9)^3 ;(*In AU^3/s^2*)

(*p-iteration constants*)
r1 = Sqrt[r[1].r[1]];
r2 = Sqrt[r[2].r[2]];
deltaA = ArcCos[r[1].r[2]/(r1 r2)];
k = r1 r2 (1 - Cos[deltaA ]);
l = r1 + r2;
m = r1 r2 (1 + Cos[deltaA ]);

(*p-iteration calculations*)
a[p_] := (m k p)/((2 m - l^2) p^2 + 2 k l p - k^2)
f[p_] := 1 - r2/p (1 - Cos[deltaA ])
g[p_] := (r1 r2 Sin[deltaA ])/Sqrt[muS p]
deltaE[p_] := If[a[p] >= 0, ArcCos[1 - r1/a[p] (1 - f[p])], Throw["a <0"]] 
t[p_] := g[p] + Sqrt[a[p]^3/muS] (deltaE[p] - Sin[deltaE[p]]) 

(* initial guesses*)
pT = {k/(l + Sqrt[2 m]) 1.1, k/(l - Sqrt[2 m] 0.95)}; 
tT = t /@ pT;
(*Recursive calc*)
newP := pT[[-1]] + (tf - tT[[-1]]) (pT[[-1]] - pT[[-2]])/(tT[[-1]] - tT[[-2]])

(* calc it *)
  AppendTo[pT, newP];
  AppendTo[tT, t[pT[[-1]]]],
{234.446, 115., 292.58, 173.493, 201.889, 208.092, 206.971, 207., 207., 207., 207., 207.}


Instead of the ugly Do[ ] you can do:

iter[{pT_, tT_}] := {Join[pT, {#}], Join[tT, {t[#]}]} &@
                         (pT[[-1]] + (tf - tT[[-1]]) (pT[[-1]] - pT[[-2]])/(tT[[-1]] - tT[[-2]]))
Nest[iter, {pT, tT}, 10]
  • $\begingroup$ I'm a little confused by the results as unfortunately they do not match the answer given in the example question on the linked website (braeunig.us/space/problem.htm#5.3, where the value of p that gives the correct answer for t is p = 1.259067). Is there a way to make FindRoot accept only an initial value for the time of flight (deltat) and p (i.e. no guesses for a and deltaE) and come out with a result for t as the algorithm in the paper implied? $\endgroup$ Dec 5, 2013 at 15:50
  • 1
    $\begingroup$ @user7388 See the new (iterative) version $\endgroup$ Dec 5, 2013 at 18:20
  • $\begingroup$ Thank you belisarius, bounty points well deserved! Now I can learn what each part of your algorithm does as my knowledge of Mathematica's functions and syntax are quite weak to say the least. Cheers! $\endgroup$ Dec 6, 2013 at 9:39
  • $\begingroup$ Hmmm, this seems to fail in quite a few cases. For example if I use r[1] = {5.9408*10^6, 3.027*10^6} and r[2] = {0, 6.9675*10^6} I get an error saying GreaterEqual::nord: "Invalid comparison with -1.06784*10^6-520653.\ i attempted.". I'm guessing this has to do with calculating deltaE where a[p] is an imaginary number (and imaginary numbers can only have their magnitudes compared). I've looked at the original source and it mentions nothing about what to do if imaginary numbers crop up though, so I'm quite lost at the moment :/ $\endgroup$ Jan 12, 2014 at 18:41
  • $\begingroup$ I forgot to add that I also set tf = 3600 $\endgroup$ Jan 12, 2014 at 20:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.