5
$\begingroup$

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.

$\endgroup$
4
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Looking closer I see that that was indeed where my problem lay and it has since been fixed. Thanks very much! $\endgroup$ – indigoblue Jul 10 '15 at 15:45
  • $\begingroup$ indigoblue, would you kindly post your newly corrected code here? $\endgroup$ – tckosvic Jul 10 '15 at 16:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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