Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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

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

share|improve this question
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)$$ – AlexR Nov 26 '13 at 15:17
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. – Michael E2 Nov 27 '13 at 16:21
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? – AlexR Nov 28 '13 at 11:54
up vote 3 down vote accepted

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]
share|improve this answer
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 (, 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? – AlexR Dec 5 '13 at 15:50
@user7388 See the new (iterative) version – Dr. belisarius Dec 5 '13 at 18:20
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! – AlexR Dec 6 '13 at 9:39
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 :/ – AlexR Jan 12 '14 at 18:41
I forgot to add that I also set tf = 3600 – AlexR Jan 12 '14 at 20:42

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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