Lambert-Solver-Flugbahnfehler mit geringer Flugzeit

Ich habe einen Lambert-Solver, der verwendet wird, um das erforderliche Delta-V für interplanetare Transferbahnen zu berechnen, aber ich bin kürzlich auf ein Problem gestoßen. Es funktioniert gut für "lange" Flugzeiten (z. B. eine 250-Tage-Flugbahn von der Erde zum Mars), aber wenn ich eine "kurze" Flugzeit (z. B. eine 45-Tage-Flugbahn von Erde zum Mars), scheint der Solver zu brechen, da seine Berechnungen nicht konvergieren. Mir ist klar, dass ein 45-Tage-TOF ziemlich unrealistisch ist, aber ich führe normalerweise eine Kurskorrektur auf halbem Weg durch eine Flugbahn durch, wenn ich also einen schnellen 90-Tage-TOF habe, dann führe ich 45 Tage nach Beginn des Fluges ein MCC durch (was ist wenn der Solver kaputt geht).

Das Folgende ist der Lambert-Löser, der eine Formel für universelle Variablen (codiert in der Mathematica-Sprache) und einen einfachen iterativen Bisektions-Löser verwendet, um das erforderliche Delta-V bei einer erforderlichen Flugzeit zu berechnen:

EDIT: Ich habe seitdem den Code repariert und der Lambert-Solver unten scheint ganz gut zu funktionieren:

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

Wie zu sehen ist, gibt es Überprüfungen, um sicherzustellen, dass Y immer positiv bleibt, daher bin ich mir nicht sicher, was das Problem sein könnte. Auch die Suche im Internet nach wissenschaftlichen Artikeln oder Buchkapiteln hat mir leider keine Lösung für dieses Problem gebracht. Jede Hilfe wäre sehr willkommen, vielen Dank.