A fundamental problem in spacecraft mission design is to find free-flight paths from one place to another that satisfy various design criteria. We explore the geometry of free-flight paths between departure and arrival points for Kepler's problem. Newton showed that these paths are conic arcs. We find the parameters for all conic paths between a departure and an arrival point as a function of one key variable called the inside angle. Once the paths are written in terms of this single parameter, then it is straightforward to find the path that takes a specified travel time (the Lambert problem) or to perform other optimizations such as minimizing the fuel costs.

We aim to solve boundary value problems for Newton's inverse square central force gravitational model. We first find the conic parameters (p,e,ω) for all conic arcs between departure and arrival points. Finding these parameters is complicated.1,2 In Sec. II, we give a new geometric solution to this problem. Our solution is direct and uses the symmetry of the central force model. We avoid using the eccentric anomaly and Kepler's equation.

As an application, we solve Lambert's problem. This problem has a long history.1–6 In the standard solution of Lambert's problem, the travel time is given, and the parameters of the conic path that solves the problem are unknown. Instead of starting with the travel time, we find the conic parameters (p,e,ω) for all conic arcs between these points as a function of a key variable called the inside angle and then find the parameters for a specific travel time. We also find conic transfer arcs that optimize other quantities of interest such as the energy needed for the transfer.

Arnold7 provided an elegant derivation of Kepler's laws, and we present that derivation here, both as a reminder of the general orbit problem and as a way to introduce the parameters of the elliptical orbits. The derivation starts with Newton's two body gravitational equation

r¨=μ|r|3r,
(1)

where r is the vector from the central mass to the secondary body, μ=MG, G is the gravitational constant, and M is the central mass, which we assume is much greater than the orbiter's mass and is, therefore, fixed at the origin.

Angular momentum is defined as the vector mr×ṙ, where m is the mass of the orbiting body, and conservation of angular momentum is established by using Eq. (1) to show that d/dt(r×ṙ)=0. Solutions with non-zero angular momentum travel in the plane orthogonal to the angular momentum vector. Let L denote the magnitude of the angular momentum vector and choose polar coordinates (r,ν) in this plane. (The symbol ν is commonly used in celestial mechanics to represent the polar angle, and it is referred to as the “true anomaly.”) Then Newton's equation in polar coordinates is the system of equations

mr2ν̇=L,r¨=L2m2r3μr2.
(2)

For these equations, one can define an “effective potential energy” V(r)=L2/(2mr2)μm/r and define energy as ϵ(r,ṙ)=mṙ2/2+V(r). Using Eq. (2), one shows that dϵ/dt=0 so that energy is conserved. For a fixed energy ϵ,

ṙ=2[ϵV(r)]/m.
(3)

Attempting to solve this differential equation by separation of variables leads to an elliptic integral that does not have an antiderivative given in terms of elementary functions. However, a satisfying alternative is found by finding the angle as a function of the radius using Eq. (2)

ν̇=dνdt=dνdrdrdt=Lmr2,
(4)
dνdr=Lr22m[ϵV(r)].
(5)

Equation (5) can be solved by separation of variables, and the equation for the angle as a function of the radius can be inverted with the result that

r(ν)=p1+ecos(ν),p=L2m2μ,e=1+2ϵL2μ2m3.
(6)

This gives Kepler's first law that orbits trace conic arcs. Using Eq. (2), the equation for the angle as a function of time is

dνdt=Lmp2(1+ecos(ν))2.
(7)

To solve Eq. (7) by separation of variables leads to an elliptic integral. The integral caused a big problem before the invention of computers, because the integral did not have an antiderivative that could be expressed in terms of elementary functions. This led to the introduction of eccentric anomalies and Kepler's equation. We will instead use the true anomaly as discussed below.

Our approach begins by first finding all possible conic arcs between two points. We start with the special case where the first point (1,ν) is on the unit circle and the second point (γ,ν+Δν) is on the circle of radius γ>1. The conic parameters (p,e,0) of all possible conic arcs with ω = 0 between these points are found as a function of the variables (γ,Δν,ν) as discussed in Sec. II.

To produce all the conic arcs from a point (R,α) to a point (γR,α+Δν), we rescale and rotate the conic arcs for the special case. As discussed in Sec. III, this general solution is used to choose conic arcs that minimize characteristics of interest such as the eccentricity or the energy of a transfer orbit in addition to the travel time.

In Sec. IV, we use a recent NASA mission to Mars to pose and solve a boundary value problem and apply our methods. In Sec. V, we discuss how to patch together conic arcs between several departure and arrival positions.

For two points with polar coordinates (R,α) and (Rγ,α+Δν) with γ>1, the first step of the problem is to find the parameters of all conic arcs between them. The equation of a general conic has the form

r(θ)=p1+ecos(θω).
(8)

The angle ν=θω is called the true anomaly. A general conic is defined by its parameters (p,e,ω). The scale parameter p and the elliptic parameter, or eccentricity, e are non-negative. A restrictedconic as we define it here has the argument of periapse ω = 0. Coordinate rotation by ω takes the restricted conic with parameters (p,e,0) to the conic with parameters (p,e,ω). Suppose a restricted conic with parameters (p,e,0) contains the points with polar coordinates (1,ν1) and (γ,ν1+Δν) as shown in Fig. 1. Here, the initial true anomaly ν1 is also referred to as the inside angle, and the change in the true anomaly Δν is called the transfer angle.

Then the rotated and rescaled conic with parameters (Rp,e,αν1) contains the points (R,α) and (Rγ,α+Δν). We will rotate and rescale to produce transfer conics from restricted conics. This will be done in the final two paragraphs of this section. Until that point, all conics will be restricted conics with ω = 0.

Given circles of radius 1 and radius γ>1 centered at the origin of a plane, we first find all ellipses that intersect both circles. The minimum and maximum radii of an ellipse with parameters (p,e,0) are equal to p/(1+e) and p/(1e), respectively. Necessary and sufficient conditions for intersections with the circles are p/(1+e)1γp/(1e). The set of points in the (p, e) plane which satisfy these conditions forms a triangle T(γ) as shown in Fig. 2 for γ=1.524. The value 1.524 was used because this is the roughly the distance of the orbit of Mars from the Sun in astronomical units. The triangle is bounded by the lines e=1,p=1+e, and p=γ(1e). For conics with e1, the condition for intersection is 0p(1+e).

A conic with parameters (p,e,0) generates two equations for its intersections with the two circles

p1+ecos(ν1)=1,
(9)
p1+ecos(ν2)=γ.
(10)

The condition for the intersection points to span the angle Δν=ν2ν1 further restricts the possible values of e and p. Solving Eqs. (9) and (10) for e and p as functions of the parameters γ, ν1, and Δν=ν2ν1 yields Eqs. (21) and (22). These equations allow all orbits with the required transfer angle Δν between the two orbital radii to be parameterized in terms of the single parameter ν1. By evaluating the orbital characteristics for the range of possible ν1, one can find the ellipse that satisfies a particular condition such as the transit time (the Lambert problem) or the minimum fuel transfer orbit.

However, before examining the results of those equations, we will first derive them through the following transformations, which can provide additional insight into the allowed values.

The cosines (c1, c2) of the angles (ν1,ν2) are found by solving the equations

p=1+ec1,
(11)
pγ=1+ec2.
(12)

The solution defines a map

f:T(γ)[1,1]×[1,1],
(13)
f(p,e)=(c1,c2)=(p1e,pγeγ).
(14)

The inverse of this map has the formula

f1(c1,c2)=(1+c1(γ1)c1γc2,γ1c1γc2).
(15)

The set T*(γ)=f(T(γ)) is also a triangle contained in the square [1,1]×[1,1]. The top boundary of T*(γ) is the image of the boundary of T(γ) with e =1. By using the formula for f1 one shows that this boundary is the line segment

c2=c1+1γ1,1c11.
(16)

The triangle T*(1.524) is shown in the unshaded region in Fig. 2. The cosine oval Oval(Δν) shown in Fig. 3 is the parameterized curve Oval(x)=(cos(x),cos(x+Δν)) for 0x2π. The arc of the oval that is contained in the interior of the triangle will be used to solve for the parameters of all ellipses that intersect the circles at points with coordinates of the form (1,ν1) and (γ,ν1+Δν). The complement of this arc may be used to solve for the parameters of all parabolic and hyperbolic transfer conics.

For a given transfer angle, draw the oval and intersect this oval with the image T*(γ) of the map f. The intersection points of the cosine oval Oval(Δν) with the upper boundary of the triangle T*(γ) are found by finding the solutions θ1 and θ2 of the equation

cos(θ+Δν)=cos(θ)+1γ1.
(17)

The intersection points determine an interval of inside angles θ1<ν1<θ2. Angles in this interval parametrize the possible elliptic transfer orbits.

A typical intersection is shown in Fig. 3 for Δν=2π/3. The preimage of the intersection Oval(Δν)T*(γ) as viewed in the p-e plane is shown in Fig. 4. This shows the elliptic parameters of the transfer ellipses with transfer angle 2π/3.

To compute the set of conic parameters f1[Oval(Δν)], one uses the function g

g:[0,2π]R2,
(18)
g(ν1)=f1(cos(ν1),cos(ν1+Δν)),
(19)
g(ν1)=(p(γ,Δν,ν1),e(γ,Δν,ν1)),
(20)
e(γ,Δν,ν1)=γ1cos(ν1)γcos(ν1+Δν),
(21)
p(γ,Δν,ν1)=e(γ,Δν,ν1)cos(ν1)+1.
(22)

We restrict the conic parameters to non-negative values. This requires that [cos(ν1)γcos(ν1+Δν)]>0. Viewed in the cosine-cosine square with c1=cos(ν1) and c2=cos(ν1+Δν), this shows a further restriction on “allowed” inside angles that generate non-negative conic parameters. The allowed angles for the non-negative conic parameters (p,e) of parabolic and hyperbolic transfers are those where the oval is below the line c2=γ1c1. The eccentricity is infinite on this line.

The allowed angles for elliptic parameters are those where the “Oval” is below the line c2=γ1c1+γ11. When the oval intersects this line, the eccentricity is equal to 1.

The function g restricted to the interval (θ1,θ2) determines the set of allowed elliptic parameters of transfer ellipses as a function of the outer radius γ and the transfer and inner angles. The function f1 applied to the set Oval(Δν) also determines the set of all conic transfer parameters.

To create elliptic arcs joining the points (R1,α) and (R2,α+Δν), choose γ=R2/R1. Using Eqs. (21) and (22), choose elliptic parameters (p(γ,Δν,ν1),e(γ,Δν,ν1)). These are the parameters of an ellipse that intersects the circle of radius 1 and the circle of radius γ with the transfer angle Δν and inside angle ν1. Set ω=αν1. The angle is the argument of the periapse of a rotated restricted ellipse. Then the transfer ellipse between the points (R1,α) and (R2,α+Δν) is the general conic

r(θ)=R1p(γ,Δν,ν1)1+e(γ,Δν,ν1)cos(θω).
(23)

When the inside angle ν1 is unconstrained, formula (23) produces ellipses, parabolas, and hyperbolas. Using Eq. (21), the eccentricity is zero when tan(ν1)=[cos(Δν)1/γ]/sin(Δν) and is infinite when [cos(ν1)γcos(ν1+Δν)]=0. We chose to focus on elliptic transfer paths, because these have the lowest energy needed for the transfer.

For a given gravitational parameter μ, travel times from an inner radius R to an outer radius γR depend on the values of (R,μ,γ,Δν,ν1). Consider travel times for restricted conics with ω = 0 as a function of these variables. Use Eqs. (21) and (22) to produce parameters e(γ,Δν,ν1),p(γ,Δν,ν1) for a conic with R =1 and the chosen transfer and inner angles. Then compute travel times using Eq. (7) and numerical integration. Alternately, use the standard approach using Kepler's equation.2 Then rescale and rotate.

Rescaling changes the travel time from τ/R3/2 to τ. This follows from Eq. (7). The time to travel from (R,α) to (γR,α+Δν) on the transfer conic

r(θ)=Rp(γ,Δν,ν1)1+e(γ,Δν,ν1)cos(θω),
(24)

with ω=αν1 is the same as the time to travel on the corresponding conic with ω = 0. Thus, the travel time does not depend on ω.

The algorithm solves Lambert's problem by first finding the conic parameters of all transfer conics and then finding the parameters for a specific travel time. Rescaling and rotating conic arcs also simplify other boundary value problems.

As an example, we consider the Mars 2020 mission launched by NASA on July 30, 2020 and landed on February 18, 2021. The transfer angle between departure and arrival points is about 143.2°, and the travel time is about 203 days. The allowable range of inside angles (θ1,θ2) for elliptic transfer orbits may be computed by first finding the intersections of the oval with the hypotenuse of the triangle as shown in Fig. 5.

By plotting the transfer time function of T(ν1), the different possible solutions between the two desired points may now be clearly seen as shown in Fig. 6. By using different initial guesses for ν1, a root finding algorithm can converge to find a zero of the function f(ν1)=T(ν1)TransferTime.

By limiting the search space to the allowable ν1 values and using the algorithm, our solution converged to ν1=17.3° and (e,p)=(0.219,1.209). The geometry of this ellipse is shown in Fig. 7. While these results were computed for the dimensionless orbits in the plane, the actual transfer orbits may be computed in three dimensions by simply scaling and rotating the resulting orbits in combination with using states from the ephemeris.

A similar algorithm can be used to minimize fuel costs, which are an important factor in mission design. We are biased towards using elliptic transfer orbits because these have the lowest energy cost. Equation 33 below shows the role played by the conic parameters in the speed at a point on a transfer orbit. Having identified all possible trajectories, a simple calculation finds the parameters of the least energy transfer. To minimize energy as a function of the inside angle, first find the scale and eccentricity parameters (p, e) as functions of the inside angle, then express the energy as a function of these parameters, and then minimize the energy as a function of the inside angle.

Equation 21 gave the eccentricity as a function of the inside angle for a fixed transfer angle. The minimum eccentricity is found at the inside angle implicitly defined by the equation

tan(ν1)=γsin(Δν)/(1γcos(Δν)).
(25)

In this section, the formulas for conic arcs and transfer times are applied to design patched conic orbits. Patched conic orbits are used for preliminary mission design where the velocity change of a spacecraft at the patch point is accomplished by applying a thrust impulse from a rocket engine. To transfer from one arc to the next on a patched conic path, one must know the angle between the tangent lines to the arcs at the patch point. The angle Φ of a tangent vector to the conic equation (31) at a point (r(ν),ν) on the conic is required for this calculation. To find this angle, one can calculate using complex variables. For this calculation the notation for eccentricity conflicts with exponential notation. Just for this calculation, the symbol used for eccentricity will be ê

r(ν)=p1+êcos(ν),
(26)
q(ν)=r(ν)eiν,q(ν)=r(ν)eiν+r(ν)ieiν=s(ν)eiΦ,
(27)
r(ν)+r(ν)i=s(ν)ei(Φν),
(28)
tan(Φν)=r(ν)r(ν),
(29)
Φ=ν+arctan(r(ν)r(ν)),
(30)
r(ν)=êpr2(ν)sin(ν),
(31)
r(ν)r(ν)=1+êcos(ν)êsin(ν).
(32)

Now replace ê by e as the symbol for eccentricity. The energy of the orbit in terms of the elliptic parameters (p,e,ω) is ϵ=μm(e21)/2p. For an elliptic orbit of the Kepler problem r¨=μr|r|3, the Cartesian velocity v(r,ν)=ṙ is the velocity at the point (r,ν) on the conic and the speed is

s(r,ν)=|ṙ|=2(ϵ+μmr)=μm(e21p+2r).
(33)

The velocity in polar coordinates is (s(r,θ,μ,e,p),Φ(r,θ,e,p,ω)). To transfer from an incoming elliptic orbit to an outgoing elliptic orbit at a common patch point (r,θ), the velocity change needed is the difference between the outgoing and incoming velocities at this point. These velocities are determined for polar coordinates by the above formulas.

As an interesting exercise, one may patch elliptic arcs together in an outward spiral. Choose conic arcs between circles of radii 1,γ,γ2, with a fixed transfer angle Δν between patch points and choose allowed elliptic parameters (p0e0,ω0). Then choose an initial scale parameter and generate a sequence of elliptic parameters p1=γp0,e1=e0,ω1=ω0+Δν, Elliptic arcs, travel times, and velocity changes can now be computed from these data. A plot for the choice γ=1.524,Δν=2π/3,p0=1.202,e0=0.2385 is shown in Fig. 8.

The problem of finding all conic paths from one position to another in three-dimensions with focus at the origin is an exercise in linear algebra. We have already shown how to find all conic paths from A to B in R3 provided both A and B are in the x-y plane. To find all conic paths from position X to position Y in R3 with focus at the origin first rotate the X-Y plane to the x-y plane around the line where these planes intersect. This line is called the line of nodes. Find the transformed points A and B. Find all conic paths from A to B and rotate the paths back to the X-Y plane. The details of this process will use the properties of vector cross products.

Let (b1,b2,b3) denote the standard basis of R3 as column vectors. Define a new orthonormal basis (B1,B2,B3) for R3 starting with B3=(X×Y)/norm(X×Y). Both X and Y are orthogonal to B3.

Define B1=b3×B3. This is the carefully chosen new first-axis in the plane spanned by X and Y. The basis vector B1 is orthogonal to both b3 and B3.

Define B2=B3×B1. This is a new second-axis in the plane spanned by X and Y. The new basis (B1,B2,B3) preserves orientation. The linear transformation from the standard basis to the new basis is given by the orthonormal matrix U=(B1,B2,B3). The inverse of the matrix U is U (the transpose of U). Transform the position vectors X and Y to the x-y plane by defining A=UX and B=UY.

These are vectors in the x-y plane. Now find all conic arcs from A to B. Then transform the arcs to arcs from position X to position Y in R3 using the matrix U. The elliptic parameters of these arcs in the x-y plane are the same as the elliptic parameters of the transformed arcs in the plane using the basis (B1,B2,B3).

To find the departure and arrival velocities in 3D for a given travel time from position X to position Y, rotate the positions into the x-y plane, and solve the 2 D problem. Then calculate the velocities at the departure and arrival positions and rotate these velocities back into the X-Y plane.

A new geometric approach to computing all conic transfer arcs between two points is derived as a function of the polar coordinates of the points. The travel time between points on a given conic transfer arc is calculated as a function of the inside angle. This method was successfully used to compute conic transfers across a range of possible transfer times for the direct case using elliptic arcs. By implementing a straightforward zero finding method, it is possible to compute specific orbits with desired transfer times. This approach solves Lambert's problem once subsequent standard rotations are performed to orient the orbit correctly in a three-dimensional space. The present approach also allows one to optimize other possible astrodynamic parameters for the transfer. This may reduce the search space and speed up the computation of the desired solution.

Part of the research presented in this paper has been carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (No. 80 NM0018 D0004). The authors would like to thank Jonathan Levine, James Meiss, and David Sterling for valuable comments and suggestions.

The authors have no conflicts to disclose.

1.
R. R.
Bate
,
D. D.
Mueller
, and
J. E.
White
,
Fundamentals of Astrodynamics
(
Dover Publications, Inc.
,
New York
,
1971
).
2.
R. H.
Battin
,
An Introduction to the Mathematics and Methods of Astrodynamics
, revised ed., AIAA Education Series (
AIAA
,
Sterling
,
1999
).
3.
R. H.
Gooding
, “
A procedure for the solution of lambert's orbital boundary-value problem
,”
Celestial Mech. Dynamical Astron.
48
(
2
),
145
165
(
1990
).
4.
J. F.
Jordan
, “
The application of Lambert's theorem to the solution of interplanetary transfer problems
,” Technical Report,
Jet Propulsion Laboratory, California Institute of Technology, California
.
5.
R. P.
Russell
, “
On the solution to every Lambert problem
,”
Celestial Mech. Dyn. Astron.
131
,
50
(
2019
).
6.
D. A.
Vallado
, “
Fundamentals of astrodynamics and applications
,” 4th ed.
Space Technology Library
(
Microcosm Press
,
Oregon
,
2013
).
7.
V. I.
Arnold
,
Mathematical Methods of Classical Mechanics
(
Springer-Verlag
,
Berlin
,
1978
).