# p-iteration algorithm to solve Lambert's problem for interplanetary trajectories

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 = {0.473265, -0.899215, 0};
r = {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.r]
r2 = Sqrt[r.r]
deltaA = ArcCos[r.r/(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:

p[n+1]=p[n]+(deltat-t[n])(p[n]-p[n-1])/(t[n]-t[n-1])


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: • 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)$$ – InquisitiveInquirer 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? – InquisitiveInquirer Nov 28 '13 at 11:54

(* Problem Data *)
G = 6.672*10^-11;
deltat = 207*86400;(*seconds*)
M = AstronomicalData["Sun", "Mass"];
m = AstronomicalData["Earth", "Mass"];
r = {0.473265, -0.899215, 0};
r = {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.r];
r2 = Sqrt[r.r];
deltaA = ArcCos[r.r/(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 *)
Do[
AppendTo[pT, newP];
AppendTo[tT, t[pT[[-1]]]],
{10}];
tT/86400
(*
{234.446, 115., 292.58, 173.493, 201.889, 208.092, 206.971, 207., 207., 207., 207., 207.}
*)


Edit

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]

• 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? – InquisitiveInquirer 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! – InquisitiveInquirer Dec 6 '13 at 9:39
• Hmmm, this seems to fail in quite a few cases. For example if I use r = {5.9408*10^6, 3.027*10^6} and r = {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 :/ – InquisitiveInquirer Jan 12 '14 at 18:41
• I forgot to add that I also set tf = 3600 – InquisitiveInquirer Jan 12 '14 at 20:42