Conic transfer arcs for Kepler's problem

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.

Arnold⁷ 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

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

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 ϵ,

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)

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

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

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

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 $restricted conic$ 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 ν₁ 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/(1−e)$⁠, respectively. Necessary and sufficient conditions for intersections with the circles are $p/(1+e)≤1≤γ≤p/(1−e)$⁠. 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=γ(1−e)$⁠. For conics with $e≥1$⁠, the condition for intersection is $0≤p≤(1+e)$⁠.

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

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 γ, ν₁, 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 ν₁. By evaluating the orbital characteristics for the range of possible ν₁, 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 (c₁, c₂) of the angles $(ν1,ν2)$ are found by solving the equations

The solution defines a map

The inverse of this map has the formula

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 $f−1$ one shows that this boundary is the line segment

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 $0≤x≤2π$⁠. 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 θ₁ and θ₂ of the equation

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 $f−1[Oval(Δν)]$⁠, one uses the function g

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+γ−1−1$⁠. 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 $f−1$ 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 ν₁. 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

When the inside angle ν₁ 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.² 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

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 ν₁, 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 ν₁ 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

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 $ê$

Now replace $ê$ by e as the symbol for eccentricity. The energy of the orbit in terms of the elliptic parameters $(p, e, ω)$ is $ϵ=μm(e2−1)/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

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 $(p0 e0,ω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 R³ provided both A and B are in the x-y plane. To find all conic paths from position X to position Y in R³ 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 R³ as column vectors. Define a new orthonormal basis $(B1,B2,B3)$ for R³ 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 B₁ is orthogonal to both b₃ and B₃.

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=U′X$ and $B=U′Y$⁠.

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 R³ 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.