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

Question

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:

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:

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 *)
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]