Exercises — Gauss's method for Lambert's problem

Inputs: two position vectors ${\mathbf{r}}_1,{\mathbf{r}}_2$ (measured from the same focus = central body) and the time-of-flight $\Delta t$ between them. Outputs: the connecting Kepler orbit, delivered concretely as the two velocities ${\mathbf{v}}_1$ (at point 1) and ${\mathbf{v}}_2$ (at point 2). From ${\mathbf{r}}_1,{\mathbf{v}}_1$ you can compute every classical orbital element, so the orbit is fully determined.

$y=\frac{\text{area of the elliptic sector swept}}{\text{area of the triangle}}$ formed by the focus and the two position points. $y\approx 1$ means the swept "pie slice" is nearly the same as the flat triangle between the points — the arc barely bulges out, which happens for small transfer angles $\theta$. Look at the figure: the yellow sector hugs the blue triangle.

The time-of-flight must enter squared: $m=\frac{\mu \Delta t^2}{(2\sqrt{r_1{r}_2}\cos (\theta /2){)}^3}$. Why squared? This is Kepler's-third-law scaling — the parent note packs the time information as $\propto \Delta t^2$. Using $\Delta t$ linearly is dimensionally and physically wrong.

What we do: apply the law of cosines to the focus–point–point triangle. Why: $c$ is the straight-line distance; everything on the right is already known. $c^2=r_1^2+r_2^2-2r_1{r}_2\cos \theta =4+9-2(6)(0.76604)=13-9.19252=3.80748,$ $c=\sqrt{3.80748}=\mathrm{1.9513.}$

$\cos (\theta /2)=\cos 20^{\circ }=0.93969$, $\sqrt{r_1{r}_2}=\sqrt{6}=2.44949$. $\ell =\frac{r_1+r_2}{4\sqrt{r_1{r}_2}\cos (\theta /2)}-\frac{1}{2}=\frac{5}{4(2.44949)(0.93969)}-\mathrm{0.5.}$ Denominator $=4(2.44949)(0.93969)=9.2065$, so $\ell =\frac{5}{9.2065}-0.5=0.54309-0.5=0.04309$.

The cube's base is $2\sqrt{r_1{r}_2}\cos (\theta /2)=2(2.44949)(0.93969)=4.60326$. $m=\frac{\mu \Delta t^2}{(4.60326{)}^3}=\frac{(1)(1.8{)}^2}{97.552}=\frac{3.24}{97.552}=\mathrm{0.033213.}$

Step A — invert Gauss's first equation to get $x$. The first equation is $y^2=\frac{m}{\ell +x}$; solving it for $x$ gives $\ell +x=\frac{m}{y^2}$, i.e. $x=\frac{m}{y^2}-\ell =\frac{0.033213}{1^2}-0.04309=-\mathrm{0.009877.}$ Step B — series $X(x)$. Keep two terms: $X(x)=\frac{4}{3}\left(1+\frac{6}{5}x+\cdots \right)=\frac{4}{3}(1+1.2(-0.009877))=\frac{4}{3}(0.988148)=\mathrm{1.31753.}$ Step C — update $y$ via $y=1+(\ell +x)(x-\frac{1}{2}+X)$ (this is Gauss's second equation divided by $y^2$, using $m/y^2=\ell +x$):

• $\ell +x=0.04309-0.009877=0.033213$
• $x-\frac{1}{2}+X=-0.009877-0.5+1.31753=0.807653$
• $y=1+(0.033213)(0.807653)=1+0.026823=\mathrm{1.02682.}$

What it means: the guess $y=1$ was slightly too small; the sector is $\sim 2.7\%$ bigger than the triangle. One more iteration barely moves it.

New $x$: $x=\frac{0.033213}{1.02682^2}-0.04309=\frac{0.033213}{1.054360}-0.04309=0.031501-0.04309=-0.011589$. New $X$: $X=\frac{4}{3}(1+1.2(-0.011589))=\frac{4}{3}(0.986093)=1.314791$. New $y$:

• $\ell +x=0.04309-0.011589=0.031501$
• $x-\frac{1}{2}+X=-0.011589-0.5+1.314791=0.803202$
• $y=1+(0.031501)(0.803202)=1+0.025302=\mathrm{1.025302.}$

$y$ moved from $1.02682$ to $1.02530$ — a change of only $0.0015$. The map is a contraction near $y\approx 1$; another step settles $y\approx 1.0253$. This is the "3–5 iterations" claim made concrete — trace it on the figure: the blue markers step from $1$ toward the green dashed fixed line, each jump smaller than the last.

Semi-latus rectum. Using the parent's extraction formula $p=\frac{r_1{r}_2(1-\cos \theta )}{r_1+r_2-2\sqrt{r_1{r}_2}\cos (\theta /2)(1-2x)}.$ Numerator: $r_1{r}_2(1-\cos \theta )=6(1-0.76604)=6(0.23396)=1.40376$. Denominator: $2\sqrt{r_1{r}_2}\cos (\theta /2)=4.60326$; $(1-2x)=1-2(-0.011589)=1.023178$; product $=4.71028$; so denominator $=5-4.71028=0.28972$. $p=\frac{1.40376}{0.28972}=\mathrm{4.8452.}$ Lagrange $f,g,\dot{g}$ (parent formulas):

• $f=1-\frac{r_2}{p}(1-\cos \theta )=1-\frac{3}{4.8452}(0.23396)=1-0.14487=\mathrm{0.85513.}$
• $g=\frac{r_1{r}_2\sin \theta }{\sqrt{\mu p}}=\frac{6\sin 40^{\circ }}{\sqrt{4.8452}}=\frac{6(0.64279)}{2.20118}=\frac{3.85674}{2.20118}=\mathrm{1.75212.}$
• $\dot{g}=1-\frac{r_1}{p}(1-\cos \theta )=1-\frac{2}{4.8452}(0.23396)=1-0.09658=\mathrm{0.90342.}$

Why these: $f,g$ are the closed-form propagator coefficients (Lagrange f and g functions), so ${\mathbf{r}}_2=f{\mathbf{r}}_1+g{\mathbf{v}}_1$ — inverting gives the velocity without re-solving Kepler's equation.

${\mathbf{v}}_1=\frac{{\mathbf{r}}_2-f{\mathbf{r}}_1}{g}.$ $f{\mathbf{r}}_1=(0.85513\cdot 2,0)=(1.71026,0)$. ${\mathbf{r}}_2-f{\mathbf{r}}_1=(2.29813-1.71026,1.92836-0)=(0.58787,1.92836)$. ${\mathbf{v}}_1=\frac{1}{1.75212}(0.58787,1.92836)=(0.33552,1.10059).$ ${\mathbf{v}}_2=\frac{\dot{g}{\mathbf{r}}_2-{\mathbf{r}}_1}{g}.$ $\dot{g}{\mathbf{r}}_2=(0.90342\cdot 2.29813,0.90342\cdot 1.92836)=(2.07623,1.74212)$. Subtract ${\mathbf{r}}_1=(2,0)$: $(0.07623,1.74212)$. ${\mathbf{v}}_2=\frac{1}{1.75212}(0.07623,1.74212)=(0.04351,0.99429).$ These are the connecting orbit's velocities — the output of Lambert's problem.

As $\theta \to 180^{\circ }$, $\theta /2\to 90^{\circ }$, so $\cos (\theta /2)\to 0$.

• $\ell =\frac{r_1+r_2}{4\sqrt{r_1{r}_2}\cos (\theta /2)}-\frac{1}{2}\to \frac{2}{4\cdot 1\cdot 0^{+}}-\frac{1}{2}\to +\infty .$
• $m=\frac{\mu \Delta t^2}{(2\sqrt{r_1{r}_2}\cos (\theta /2){)}^3}\to \frac{\mu \Delta t^2}{(0^{+}{)}^3}\to +\infty .$

Physical reason: the two points are becoming antipodal (opposite sides of the focus). The transfer plane is no longer uniquely defined by ${\mathbf{r}}_1,{\mathbf{r}}_2$ — infinitely many orbital planes contain both points — so the geometry degenerates. Gauss's half-angle parameterization sits exactly on this singularity. Fix: switch to a solver that has no half-angle in the denominator — the Universal variable formulation Lambert solvers (Battin's method, Izzo Lambert solver) stay well-behaved through $180^{\circ }$. Rule of thumb: Gauss is best for $\theta \lesssim 90^{\circ }$–$120^{\circ }$.

• $\cos 320^{\circ }=\cos 40^{\circ }=0.76604$ — unchanged (cosine is even about $360^{\circ }$). So the chord $c$ from the law of cosines is identical! That is why the dot product alone cannot distinguish the branches.
• $\theta /2=160^{\circ }$, $\cos 160^{\circ }=-0.93969$ — now negative.

Then $2\sqrt{r_1{r}_2}\cos (\theta /2)=-4.60326$ (negative). Consequences:

• $\ell =\frac{5}{4(2.44949)(-0.93969)}-0.5=\frac{5}{-9.2065}-0.5=-0.54309-0.5=-\mathrm{1.04309.}$
• $m=\frac{3.24}{(-4.60326{)}^3}=\frac{3.24}{-97.552}=-0.033213$ (negative).

What it means: the sign flip encodes the long-way geometry. The iteration and extraction formulas must be run consistently with these signed constants — you cannot mix a short-way $\cos (\theta /2)$ with a long-way $\cos \theta$. This is the mathematical face of the L1 "short vs long" trap.

A porkchop plot fixes a grid of (departure date, arrival date) pairs. Each pair fixes ${\mathbf{r}}_1$ (planet A at departure), ${\mathbf{r}}_2$ (planet B at arrival), and $\Delta t$ (their date difference) — exactly the three Lambert inputs. Solving Lambert at every grid cell yields ${\mathbf{v}}_1,{\mathbf{v}}_2$; subtracting the planets' own velocities gives the hyperbolic excess speeds, and $C_3=v_{\infty }^2$ (departure energy) is what the contours color. So the plot is literally a millions-fold sweep of Lambert's problem, and Gauss's method (or its robust cousins near $180^{\circ }$) is the engine underneath.

• Inputs/outputs? ::: In: ${\mathbf{r}}_1,{\mathbf{r}}_2,\Delta t$. Out: ${\mathbf{v}}_1,{\mathbf{v}}_2$.
• What is $y$? ::: sector-area ÷ triangle-area, $\approx 1$ for small $\theta$.
• Why square $\Delta t$? ::: Kepler-third-law scaling in $m$.
• Where does Gauss die? ::: $\theta \to 180^{\circ }$, $\cos (\theta /2)\to 0$, $\ell ,m\to \infty$.
• How many iterations? ::: 3–5 iterations to reach $10^{-8}$ accuracy, because near $y\approx 1$ the update map is a contraction.