# Low Time-of-Flight Lambert Solver Trajectory Errors

I have a Lambert solver that is used to calculate the required delta-V for interplanetary transfer trajectories, but I recently ran into a problem. It works well for "long" time-of-flight trajectories (for example, a 250 day trajectory from Earth to Mars), but if I try a "short" time-of-flight trajectory (for example, a 45 day trajectory from Earth to Mars), the solver seems to break, by having its calculations not converge. I realise that a 45 day TOF is rather unrealistic, but I usually perform a mid-course correction halfway through a trajectory, so if I have a fast 90 day TOF, then I'll perform a MCC 45 days into the flight (which is when the solver breaks).

The following is the Lambert solver which uses a Universal Variables formulation (coded in the Mathematica language) and a simple iterative bisection solver to calculate the required delta-V given a require time-of-flight:

EDIT: I've since fixed the code and the Lambert solver below seems to work quite well:

G = 6.672*10^-11; (*Gravitational Constant*)
m[0] = 1.988544*10^30; (*Mass of Sun*)
TOF = (45) (86400); (*Time-of-flight in seconds*)
R[1] = {-1.1751563715176448*^11, 8.982523733108002*^10} (*Heliocentric position of Earth*)
R[2] = {-5.256112524631399*^10, -2.1604439066188406*^11} (*Heliocentric position of Mars*)
R1 = Sqrt[R[1].R[1]]
R2 = Sqrt[R[2].R[2]]
\[CapitalDelta]\[Nu] = ArcCos[R[1].R[2]/(R1 R2)] (*Change in true anomaly, in radians*)
A = Sqrt[R1 R2 (1 + Cos[\[CapitalDelta]\[Nu] ])];

iterationCount = 0;
z = 0;
zhi = 4 \[Pi]^2;
zlow = -4 \[Pi];
c[z_] := If[z > 0, (1 - Cos[Sqrt[z]])/z,
If[z < 0, (1 - Cosh[Sqrt[-z]])/z, 1/2]]
S[z_] := If[z > 0, (Sqrt[z] - Sin[Sqrt[z]])/Sqrt[z^3],
If[z < 0, (Sinh[Sqrt[-z]] - Sqrt[-z])/Sqrt[(-z)^3], 1/6]]
Y[z_] := R1 + R2 - (A (1 - S[z] z))/Sqrt[c[z]];
X[z_] := Sqrt[Y[z]/c[z]];
t[z_] := (X[z]^3 S[z] + A Sqrt[Y[z]])/Sqrt[G m[0]];

t[z] = t[z]; (*Initial value for t[z]*)

(*Iterative Bisection Solver*)
While[Norm[t[z] - \[CapitalDelta]t] > 1*10^-6 && iterationCount < 100,
c[z_] :=
If[z > 0, (1 - Cos[Sqrt[z]])/z,
If[z < 0, (1 - Cosh[Sqrt[-z]])/z, 1/2]];
S[z_] :=
If[z > 0, (Sqrt[z] - Sin[Sqrt[z]])/Sqrt[z^3],
If[z < 0, (Sinh[Sqrt[-z]] - Sqrt[-z])/Sqrt[(-z)^3], 1/6]];
Y[z_] := R1 + R2 + (A (S[z] z - 1))/Sqrt[c[z]];
(*Making sure Y>0*)
While[A > 0 && Y[z] < 0,
zlow = zlow + 0.01;
z = (zhi + zlow)/2;
c[z_] :=
If[z > 0, (1 - Cos[Sqrt[z]])/z,
If[z < 0, (1 - Cosh[Sqrt[-z]])/z, 1/2]];
S[z_] :=
If[z > 0, (Sqrt[z] - Sin[Sqrt[z]])/Sqrt[z^3],
If[z < 0, (Sinh[Sqrt[-z]] - Sqrt[-z])/Sqrt[(-z)^3], 1/6]];
Y[z_] := R1 + R2 + (A (S[z] z - 1))/Sqrt[c[z]];
];
X[z_] := Sqrt[Y[z]/c[z]];
t[z_] := (X[z]^3 S[z] + A Sqrt[Y[z]])/Sqrt[G m[0]];
If[t[z] <= \[CapitalDelta]t, zlow = z, zhi = z];
z = (zhi + zlow)/
2; (*Re-calculating z using bisection root finding method*)
Print[Norm[t[z] - \[CapitalDelta]t]];
iterationCount++;];

f = 1 - Y[z]/R1;
g = A Sqrt[Y[z]/(G m[0])];
gdot = 1 - Y[z]/R2;
vLambert[1] = (R[2] - f R[1])/g; (*Required velocity at start of trajectory*)
vLambert[2] = (gdot R[2] - R[1])/g; (*Velocity at target arrival*)


As can be seen, there are checks in place to make sure Y always stays positive, so I'm not too sure what the problem could be. Searching the internet for any academic articles or book chapters has unfortunately not given me a solution to this problem either. Any help would be greatly appreciated, thanks very much.

I don't have experience with Mathematica, but at least the line While[A > 0 && Y[z] < 0, zlow = zlow + 1]; seems to be strange. Increasing the lower bound for z without updating your criteria will result in an endless loop once it is started! A is a constant of the specific problem and Y depends on z, not on zlow.