LAMBERT’S PROBLEM

Interplanetary Mission Design

Kate Davis

INTRODUCTION

Lambert’s problem is a way to solve for the trajectory connecting two position vectors with a given time
of flight. In Figure 1, ~r0 and ~rf define the positions of the initial planet (P0 ) at the time of departure and
the final planet (Pf ) at the time of arrival. With these positions and the time of flight, solving Lambert’s
problem will define the orbital elements of the transfer orbit. Once the orbital elements of the transfer orbit
are computed, the velocities at departure and arrival may be found.

Final Planet, Pf
at time tf

!rf
Transfer Time: ∆t
∆t = tf − t0

∆ν

!r0

Initial Planet, P0
at time t0

Figure 1. Lambert’s Problem

Lambert’s Theorem

According to Lambert’s Theorem, the transfer time ∆t from P0 to Pf is independent of the orbit’s eccen-
tricity and depends only on the sum of the magnitudes of the position vectors, the semimajor axis a and the
length of the chord joining P0 and Pf .
There are many solutions to Lambert’s problem that may be found in literature. This handout will detail
the Universal Variables Formulation. The definitions of the variables used in this handout may be found in
Table 1.

1
 Table 1. Notation
a = Semimajor axis
e = Eccentricity
E = Eccentric Anomaly
i = Inclination
k = Integer number of revolutions
r = Position Vector Magnitude (r = ~r/|~r|)
∆t = Time of flight
µ = Gravitational Parameter
ν = True Anomaly

FORMULATION
Begin with the general form of Kepler’s equation:
s 
a3

tf − t0 = ∆t = 2πk + (Ef − e sin Ef ) − (E0 − e sin E0 ) , (1)
µ

where the subscripts ‘0’ and ‘f ’ refer to the initial and final states, respectively. In this formulation, we
are not concerned with multiple revolutions, so we will ignore the 2πk term. The equation is rearranged to
obtain:
√ √  
µ∆t = a3 ∆E + e (sin E0 − sin Ef )
√ √ √
µ∆t = a3 ∆E + a3 e (sin E0 − sin Ef ) . (2)

We will the use universal variable χ to rearrange Equation 2 into quantities that are already known.
Although not provided here, Bate, Mueller, & White1 and Vallado2 provide detailed explanations on how to
derive the universal variable χ. The expression for χ is
√ √
χ= a (Ef − E0 ) = a ∆E. (3)

Furthermore, define the constants c2 and c3 as

1 − cos ∆E
c2 = (4)
∆E 2
∆E − sin ∆E
c3 = . (5)
∆E 3

The quantity χ3 c3 is computed and rearranged:

√

3 ∆E − sin ∆E
χ3 c3 = a∆E 3
√ √ ∆E
3
χ c3 = a3 ∆E − a3 sin ∆E
√ √
3
a3 ∆E = χ c3 + a3 sin ∆E. (6)

Insert Equation 6 into Equation 2 to obtain:

√ √ √
µ∆t = χ3 c3 + a3 sin ∆E + a3 e (sin E0 − sin Ef ) . (7)

2
Use the trigonometric identity

sin ∆E = sin Ef cos E0 − cos Ef sin E0 , (8)

and Equation 7 becomes

√ √
µ∆t = χ3 c3 +


a3 sin Ef cos E0 − cos Ef sin E0 + e sin E0 − e sin Ef
√ √ 
µ∆t = χ3 c3 + a3 sin E0 (e − cos Ef ) − sin Ef (e − cos E0 ) .

(9)

The sine and cosine functions for true anomaly are given by:

e − cos E
cos ν = (10)
e cos E − 1
√
1 − e2 sin E
sin ν = (11)
1 − e cos E
We wish to express the terms in Equation 9 in terms of the eccentric anomaly and the true anomaly. There-
fore, we will multiply the final term in Equation 9 by β where:
√
1 − e2 (1 − e cos E0 )(1 − e cos Ef )
β=1= √ (12)
1 − e2 (1 − e cos E0 )(1 − e cos Ef )

Multiplying by β and rearranging gives:

√
√ √ 
sin E0 (e − cos Ef ) − sin Ef (e − cos E0 ) (1 − e cos E0 )(1 − e cos Ef ) 1 − e2

3
µ∆t = χ c3 + a3 √
(1 − e cos E0 )(1 − e cos Ef ) 1 − e2
(13)
Now, gather the terms in Equation 13 into terms consistent with Equations 10 and 11:
"√ √ #
√ √ 1 − e 2 sin E
0 e − cos E f 1 − e 2 sin E
f e − cos E 0
3
µ∆t = χ c3 + a3 − α, (14)
(1 − e cos E0 ) (1 − e cos Ef ) (1 − e cos Ef ) (1 − e cos E0 )

where
(1 − e cos E0 )(1 − e cos Ef )
α= √ . (15)
1 − e2

Write Equation 14 in terms of sin ν and cos ν:

√ √
(1 − e cos E0 )(1 − e cos Ef )

 
33
µ∆t = χ c3 + a sin ν0 (− cos νf ) − sin νf (− cos ν0 ) √ . (16)
1 − e2
√ √ √
Use the trig identity given in Equation 8 and the fact that a4 / a = a3 to obtain:

√ a(1 − e cos E0 )a(1 − e cos Ef )

µ∆t = χ3 c3 + sin(∆ν) p . (17)
a(1 − e2 )

The magnitude of the position vector can be expressed in terms of a, e, and E as follows:

r = a(1 − e cos E). (18)

3
Therefore,
√
√ r0 rf sin ∆ν 1 − cos ∆ν
3
µ∆t = χ c3 + p √ (19)
a(1 − e2 ) 1 − cos ∆ν
√ √ √
√ 3
r0 rf sin ∆ν r0 rf 1 − cos ∆ν
µ∆t = χ c3 + √ p . (20)
1 − cos ∆ν a(1 − e2 )

Let the variables A and y be defined as:

√
r0 rf sin ∆ν
A = √ (21)
1 − cos ∆ν
r0 rf (1 − cos ∆ν)
y = . (22)
a (1 − e2 )

and Equation 20 may be written as

√ √
µ∆t = χ3 c3 + A y. (23)
In terms of the time of flight:
√
χ3 c3 + A y
∆t = √ . (24)
µ
This is the universal variable form that can be found in Bate, Mueller, & White, and Vallado. Many of the
variables can also be expressed equivalently using the variable c2 given in Equation 4:
r
y
χ = (25)
c2
q
A = DM r0 rf (1 + cos ∆ν) (26)
A ∆E 2 c3 − 1


y = r0 + rf + √ , (27)
c2

where 
+1 if ∆ν < π
DM = Direction of Motion =
-1 if ∆ν > π.
Many texts also replace the quantity ∆E 2 with ψ (ψ = ∆E 2 ).
A simple way to derive these equation in terms of ψ is to use the f and g expressions. The derivations
for f and g may be found in Bate, Mueller, & White (see Chapter 4). Then, one can iterate on ψ until
the desired time of flight is obtained. Vallado uses a secant method. The method is presented in § 7.6, pp.
463 (for the 2nd edition) and is also given in the next section. The solution to Lambert’s problem then is
an iterative process on ψ, where the variables χ, c2 , c3 , and y are computed for each new value of ψ. The
process continues until the desired time of flight is achieved. Note that the value of A is only a function of
the magnitude of the two positions and the change in the true anomaly, and is not updated.

4
ALGORITHM

Input

~r0 = Initial Position Vector

~rf = Final Position Vector
∆t0 = Desired Transfer Time
DM = Direction of Motion∗∗ (−1 or +1)

∗∗ An input of DM is OPTIONAL! You may pass it in if you want to force the direction of motion to be
(−1 or +1). Alternatively, you can use the value of ~r0 and ~rf to compute ∆ν, which gives you the DM (see
hints page).

Initialize the following values

Optional Computation: If DM is not an input

∆ν = ν2 − ν1
∗∗ See hints page on computing ∆ν

IF ∆ν < π,
DM = 1
ELSE
DM = −1
END IF

~r0 • ~rf
cos ∆ν =
|~r0 ||~rf |
p
A = DM r0 rf (1 + cos ∆ν)

IF ∆ν = 0, A = 0
Trajectory can’t be computed.
END IF

1 1
c2 = , c3 =
2 6
ψ=0
ψup = 4π 2

ψlow = −4π

5
WHILE |∆t − ∆t0 | > 1 × 10−6
A (ψc3 − 1)
y = r0 + rf + √
c2
IF (A > 0.0 & y < 0.0)
Readjust ψlow until y > 0.0 ∗∗ See Hints Section on ways to adjust ψ

END IF
y
r
χ=
c2
√
χ3 c3 + A y
∆t = √
µ
IF (∆t ≤ ∆t0 )
ψlow = ψ
ELSE
ψup = ψ
END IF
ψup + ψlow
ψ=
2
IF (ψ > 1 × 10−6 )
√ √ √
1.0 − cos ψ ψ − sin ψ
c2 = , c3 = p
ψ ψ3
ELSE IF (ψ < −1 × 10−6 )
√ √ √
1.0 − cosh −ψ sinh −ψ − −ψ
c2 = , c3 = p
ψ (−ψ)3
ELSE
1 1
c2 = , c3 =
2 6
END IF
END WHILE

Compute
y y y
r
f =1− , ġ = 1 − , g=A
r0 rf µ

Output
~rf − f~r0 ġ r~f − ~r0
~v0 = , ~vf =
g g

6
Helpful Hints

• You do not necessarily have to pass the value for DM into the code. You may want to pass it in,
depending on mission requirements or you may want to compute it with the code. If you are going
to be running the code many times within loops, such as for making Pork Chop Plots, it is probably
easiest to compute DM inside the code itself. Although the quantity cos∆ν is computed within the
algorithm, do NOT take the arc cosine to get the value of ∆ν; this will have a sign ambiguity.
There is a good way to approximate the values for ν1 and ν2 which will account for the signs. If we
assume that the positions of the planets are in the ecliptic, then we can use the following to compute
ν1 and ν2 :
ry
tan ν = (28)
rx
In MATLAB, there is a function called atan2 which will guarantee the correct quadrant when you
take the arctangent. You can use that function to compute ν1 and ν2 , and then compute ∆ν and DM.
Remember to do a check to make sure 0 < ∆ν < 2π.
NOTE: Because we assumed the planets positions are in the ecliptic, if you take the cosine of the
value of ∆ν that is computed using Equation 28, the value will be slightly different than the quantity
cos∆ν that is computed within the algorithm using the dot-product, since it uses all three components
of the planets positions. However, this assumption is generally valid for most of the applications in
this class.

• Sometimes, the solver won’t converge very quickly, or it won’t converge at all for very long or very
short transfers, e.g., you might not be able to find a solution for an Earth-Saturn transfer that takes 12
hours. So it’s probably a good idea to put some internal checks into your Lambert’s solver to prevent
it from even attempting a ridiculous transfer. Don’t be too concerned if it doesn’t work sometimes,
but it should always work for reasonable Type I and Type II transfers.

• In the algorithm shown above, you’ll notice the line

Readjust ψlow until y > 0.0

This is taken directly from Vallado’s algorithm. The algorithm in Bate, Mueller, and White, simply
says, “Adjust the trial value of ψ ...” Neither algorithm give a method for exactly how to adjust ψ. In
the class, people have used different methods to adjust ψ. The reason for adjusting ψ is to make sure
that y is positive because we take the square root of it later.
We have found that simply incrementing ψ by 0.1 seems to work:
while y < 0
ψ = ψ + 0.1
Resolve for y
end
Recall that y is a function of c2 and c3 , which change as ψ changes. So, within the while loop, you’ll
need to recompute values for c2 and c3 .
There is another method that works more quickly: Since we know that we want to drive y > 0, set y
equal to 0 in the equation:
A (ψc3 − 1)
y = r0 + rf + √
c2

7
 and solve for ψ. Then, you can set multiply that value of ψ by some factor N to ensure that y > 0. I
use a factor N = 0.8. You’ll still need to use the while loop and recompute values for y, c2 , and c3 .
 √ 
c
ψ = N c13 1 − A2 (r0 + rf )

As a side note, when you execute the code, there are very few times when you actually have to adjust
ψ, so don’t stress too much over how you increment ψ.
Also, if you think of a good way to increment ψ, let us know!

REFERENCES
[1] R. R. Bate, D. D. Mueller, and J. E. White, Fundamentals of Astrodynamics. Dover Publications, Inc., New York,
1971.
[2] D. A. Vallado, Fundamentals of Astrodynamics and Applications. McGraw-Hill Companies, Inc., 1997.