The Oblate Lambert Problem: Geometric Formulation and Solution of an Unperturbed, Generalized Lambert Problem Governed by Vinti’s Potential

Numerous methods exist for solving the Lambert problem, the two-point boundary value problem (BVP) governed by two-body dynamics. Many applications would benefit from a solution to a perturbed Lambert problem; a few studies have attempted to solve one. Establishing a larger pool of alternative solution methods gives practitioners greater latitude in choosing the solution that best suits their needs. To that end, a novel Lambert-type BVP is constructed in this work that includes oblateness by way of Vinti’s potential, rendering the problem mathematically unperturbed. This BVP is first defined and then converted to a system of equations that is amenable to an iterative solution. The formulation, which is valid for both the zero- and multiple-revolution problems, couples oblate spheroidal (OS) universal variables and OS equinoctial orbital elements together to sow robustness across all orbital regimes, only excepting orbits that are sufficiently rectilinear. For the first time, the solution space is broadly explored, exposing multiple new insights of significant practical use. Initial guess and root-solve techniques are offered to solve the system of equations. When assessed at Earth for robustness, accuracy, and computational efficiency, the zero-revolution algorithm excels across all three performance metrics, with runtimes averaging only about 15 times slower than a typical two-body Lambert solver. The multiple-revolution algorithm, while not yet evaluated as extensively, also exhibits high levels of performance, the formulation generally characterizing the existence of solutions around oblate bodies more accurately than its Keplerian counterpart.


Introduction
After more than 200 years of study, the Lambert problem, a gravitational twobody boundary value problem (BVP), is considered one of the most fundamental and well understood problems of astrodynamics [1,2].Each solution draws a ballistic arc between two given terminal positions for a specified time of flight (TOF), and each arc traces the path of a conic.The Lambert problem is a twopoint BVP, and while it is limited in the fidelity of the gravitational model, it still enjoys broad use across the entire discipline.For example, preliminary orbit determination (OD) and uncorrelated track (UCT) association support critical tasks tied to the maintenance of the space object catalog and the space situational awareness (SSA) mission [3,4].Preliminary OD also supports satellite operators wishing to initiate a track of their satellite.Other applications include guidance, targeting, rendezvous, and maneuver planning efforts that are essential to the design of Earth and interplanetary missions.Practicing astrodynamicists benefit from publicly available source code that solves the Lambert problem, including Gooding's benchmark [5] Fortran implementation, as well as recent ones like that proposed by Izzo [6] or Russell [7], wherein seminal and competitive Lambert solution methods are also reviewed.
For decades, practitioners have been able to get by with the level of fidelity included in a traditional Lambert solver, but these solvers are nonetheless rendered less effective for some important applications.For example, consider the demands of active space debris removal (ASDR), which was studied by Cerf [8].This type of mission would typically be solved as a long-duration, multi-leg mission, often requiring a spacecraft to make many revolutions (revs) around the Earth.In ASDR, a spacecraft may need to rendezvous with many pieces of debris over the course of 100s of revolutions, where the unknown sequence of encounters is optimized.Under the hood, each rendezvous in turn requires the design of spacecraft maneuvers in the form of a guess, usually furnished by a Lambert solver.While the Lambert solver offers a good approximation for short orbit transfers and certain Earth-centered problems, over many revolutions, the effects of dynamical forces that are typically neglected accumulate, giving rise to large errors in a Lambert solution that subsequently stress or even break the optimizer.A rendezvous sequence generated from a standard Lambert algorithm can completely fall apart under the real dynamics, raising questions on how these issues can be mitigated.
Different strategies exist for addressing the described drawbacks that can be encountered with traditional Lambert solvers.Some may seek alternative optimizers that are less sensitive to the quality of a guess, but other strategies aim to lift the burden off the optimizer and place it squarely on the Lambert solver itself, i.e. improving the guess passed on to the optimizer without changing the optimization algorithm.The present work investigates the latter option, though it is generally worth considering both approaches since these architectures can fall in completely different regions on the speed versus robustness trade space.Computational speed and robustness are important because a Lambert solver may be called a billion times for one study.With typical Lambert runtimes around a microsecond, a billion calls could take about 20 minutes, and a 0.1% fail rate would lead to a million failures [7].Beginning with the early work of Andrus [9] and Engels and Junkins [10], the literature shows considerable interest in recent years in improving Lambert guesses [11][12][13][14][15][16], pursuing solutions to perturbed Lambert problems that include other phenomena in the BVP dynamics.The unperturbed Lambert problem is considered to be quite challenging, and the level of difficulty in solving the problem only increases with the addition of perturbations.For this reason, it is common to limit the fidelity to include only the most dominant perturbation [11,10], which for the Earth and many bodies is due to oblateness or the J 2 spherical harmonic coefficient.Figure 1 gives a notional illustration of the ballistic arc connecting two terminal positions over many revolutions when the object is influenced by a large equatorial bulge.Instead of tracing out an ellipse in a fixed orbital plane over multiple revolutions, the object's path is characterized by a distinct drift of the orbital plane.A few studies have looked at including additional perturbations [12][13][14][15][16], but they all use a classical, unperturbed Lambert solver to warm start their perturbed Lambert algorithms.
The perturbed Lambert problem can be solved in a number of ways.A straightforward approach may use a shooting method (SM) to differentially correct the unperturbed Lambert solution to target the desired final position.Shooting methods commonly use numerical integration in the trajectory propagation step, though advanced implementations leverage the speed of analytical propagators [12,17] and, with additional effort, an associated analytical state transition matrix (STM) [17], where the propagator and STM employed in the shooting method for these two studies are based on Vinti's asymmetric gravitational potential [18].Whether these shooting method components are computed numerically or analytically, the approach is still limited by the number of spacecraft revolutions undertaken in transit and can fail to converge when the number of revolutions is too high, numerical integration posing the additional disadvantage that the runtime increases with the number of revs.A multiple shooting method [17], which breaks long arcs into shorter ones, can remove this limitation, but tradeoffs can be expected between computational speed and robustness, and also the difficulty of implementation.A single shooting method can be combined with a continuation method that nudges J 2 from zero (unperturbed) up to 1.08 × 10 −3 Fig. 1 An object's ballistic arc in red (solid line) connecting two terminal positions in blue and green (dots) when influenced by Earth's oblateness.Oblateness is exaggerated for clarity (Earth-perturbed) to help circumvent the limitations on the number of spacecraft revolutions [11], but again, similar trade-offs persist.With these ideas in mind, the development of techniques to solve perturbed Lambert problems continues to be a major challenge and an active area of research.
The present study takes a markedly different approach to solving the J 2 -per- turbed Lambert problem.Instead of leveraging the two-body Lambert solution and adding perturbations on top of it, the symmetric Vinti potential [18][19][20] is used to define a novel unperturbed BVP that already includes the effects due to J 2 and a partial J 4 .This approach is inspired by the historically significant utility of exactly solvable problems, which typically betray a special structure or symmetry that can be exploited for gains in computational efficiency [21].Since the present problem is unperturbed, a solution may be developed without recourse to a shooting or other perturbation method, therefore holding promise for significant speedups in runtime, possibly with a more favorable trade in robustness than seen for other methods.A problem formulation and solution is sought and developed that is broadly applicable, not limited by the orbital regime or number of spacecraft revolutions.To achieve this goal, the study is broken down into two parts: the first part defines and sets up the BVP to be solved; the second part presents a specific method of solving the BVP and includes descriptions of the solution space and various performance results and examples.These ideas borrow heavily from universal Vinti orbit propagation methods concurrently developed by the author [22], solving the initial value problem (IVP) governed by Vinti's symmetric potential.This approach establishes techniques required to avoid computational difficulties commonly encountered in orbital mechanics and conveniently unifies solutions for bounded (elliptical) and unbounded (hyperbolic) Vinti trajectories.The analytical propagator developed by Biria [22] is notably validated against numerically integrated solutions to lend confidence to the correctness of the dynamics modeled in the BVP.Note that various solution methods packaged with publicly available source code exist for the Vinti IVP [18,[23][24][25][26].
A Lambert solution method presented in Bate et al. [1], notable for its robustness and applicability across all orbital regimes, uses universal variables (UVs) similar to those used by the author to solve Vinti's IVP [22].Therefore, much of the implementation of the author's UV Vinti solution may be directly reused in a congruent BVP solver if an approach is chosen that is similar to the one in Bate et al. [1].This choice allows the new BVP solver to inherit all of the benefits of the author's Vinti IVP formulation.Battin [27] notes that this particular unperturbed Lambert solution method was proposed by John Deyst of the MIT Instrumentation Laboratory, presumably in the early 1960s.The solution method of the current investigation essentially augments Deyst's algorithm to accommodate J 2 effects.Deyst's approach requires a root-solve of a certain system of equations.A main result of this paper is the conversion of the proposed BVP into a specific system of equations that aligns with those that Deyst identified for the two-body BVP.Then, a method is offered for iteratively solving this system of equations.Comparisons are made to the twobody BVP to illustrate the different ways that oblateness changes the solution space.Practitioners can benefit from this additional insight.Performance is also assessed through a number of examples, and strengths and pitfalls are identified.
For improved readability, an outline of the paper structure is offered before moving on.First, in Sect.2, the kinematic equations developed in prior work are reviewed and all quantities are defined.The next three sections respectively address the formulation, exploration, and solution of the BVP.Specifically, a novel formulation of a Lambert-type BVP under Vinti's potential is derived in Sect.3, the BVP solution space is explored in Sect.4, and a solution method for this BVP is proposed and evaluated in Sect. 5. Concluding remarks are offered in Sect.6.

Review of Kinematic Equations
Recall Vinti's symmetric gravitational potential [19,20] expressed in terms of oblate spheroidal (OS) coordinates , , as where is the gravitational parameter of the primary body, is the semiminor axis of the instantaneous oblate spheroid that intersects the spacecraft, is the right ascension, and is tied to latitude, approximately the sine of the declination in an Earth application.A free parameter c is fit to the oblateness coefficient, J 2 , of the spherical harmonic expansion as where R e is the equatorial radius.Refer to Vinti [19,20] for formal definitions of these quantities.The analytical solution to Vinti's integrable IVP can be stated classically as x = ̃ (t, t 0 , x 0 ) , where f is a vector-valued, nonlinear function representing a solution to the kinematic equations, such as offered by Vinti [20] or Getchell [26], t is the time, x = [r ⊤ v ⊤ ] ⊤ = [x y z ̇x ̇y ̇z] ⊤ is the Earth-centered inertial (ECI) state, r and v are the position and velocity vectors, and the subscript " 0 " denotes initial conditions (ICs).
The kinematic equations contained in f explicitly depend on anomalistic angles, but it is possible to develop alternative formulations with inherent advantages.For example, Biria [22] developed an alternative form of the kinematic equations that depends on the difference in anomalistic and other angles, and rigorous testing found that algorithms based on the new equations perform with heightened robustness in popular orbital regimes, among other benefits.The new form is referred to in this work as the "differential form" because the underlying equations depend on the difference in anomalistic and other angles instead of depending on the angles themselves.Letting Δt = t − t 0 be the time of flight, the differential form of Vinti's IVP can be stated as x = f(Δt, x 0 ) Biria [22].This ostensibly small adjustment to the form of the kinematic equations has large consequences.The refinement enables a novel formulation of the BVP governed by Vinti's potential, and it is considered foundational to the present work.
The differential form of the third-order kinematic equations [22] expressed in spheroidal UVs [28] is repeated here for convenience as (1) where subscripts indicating ICs have been changed from " 0 " to " 1 " to be consistent with typical BVP notation.The IC subscript convention applies only to 1 and 1 , where is tied to the dot product of the OS position and relative OS velocity [22,28].Equations ( 2), (3), and (4) govern anomalistic motion in the OS orbital plane, OS apsidal drift, and OS nodal drift, respectively.Next, over the remainder of the section, the large number of quantities introduced in Eqs.(2)(3)(4)) are all defined.
The solution of the IVP yields six constants of integration [18]: 1 or is the total energy, 2 is closely related to the total angular momentum, 3 is the polar component of angular momentum, = − 1 is the time of OS periapsis passage, = 2 is the argu- ment of OS periapsis, and Ω = 3 is the right ascension of the OS ascending node (spheroidal RAAN).The usual formulas [18] relate the j constants to a set of orbital elements termed the prime constants [18]: p 0 is the prime OS semi-latus rectum, 0 is the inverse of the prime OS semimajor axis a 0 , where 0 = −1∕a 0 for elliptical orbits, and S 0 = sin 2 i 0 , where i 0 is the prime OS inclination.When these elements appear without the "0" subscript, they denote a different set termed the mutual constants [18] that are related to the two quartics, F( ) and G( ) , noting that Q = sin I and a capital "I" denotes the mutual OS inclination.Specifically, p, , A 1 , and B 1 result from factoring the F( ) quartic, and Q and Q 1 result from factoring the G( ) quartic.The mutual OS eccentricity, e, is derived from p and .The factoring of the quartics also yields 1 and S 1 , two additional quantities defined by the relations 1 = 0 and S 1 S = S 1 Q 2 = S 0 , respectively [26].By convention, and because the mutual constants are used to set up the Vinti theory, references to OS elements in this work denote the mutual constants, not the prime constants, and elements should generally be interpreted as spheroidal or OS unless noted otherwise.From Getchell's formulation [26], the A n constants arise from the expansion of a quadratic in , and the C n , C 1n , and C 2n constants arise from (2) , different expansions leading to expressions in terms of Q 1 , where odd C n are zero, , C 1n and C 2n are given by d m are zero for odd m , , and d 6 = 5Q 3 1 ∕16.The time of flight and a number of important angular quantities are introduced in Eqs.(2-4) as well.First, note that the time of flight, Δt ≡ t 2 − t 1 , appears explic- itly on the left-hand side (LHS) in Eq. ( 2), where the Δ symbol denotes the differ- ence of two quantities or expressions, in this context marking how much a quantity changes during the time Δt .For example, the LHS quantities in Eqs. ( 3) and ( 4) denote two angular differences that accrue over Δt .The OS apsidal drift [22], Δ � , tracks the drift of the OS perifocal frame [28] in the OS orbital plane, and the OS nodal drift [22], ΔΩ � , tracks the drift of the OS equinoctial frame [23].Note that ′ is a different argument of OS periapsis [25] ( ′ ≠ ) and Ω � is a different spheroi- dal RAAN [29] ( Ω � ≠ Ω ).The right-hand side (RHS) in Eq. ( 2) contains additional angular variables.The W n and T n symbols are recursive functions defined by Getch- ell [26] and the computation of their differential forms is discussed extensively by Biria [22].Notably, ΔW 0 = Δf is the change in OS true anomaly and ΔT 0 = Δ that of the true argument of OS latitude.Related to the OS anomalistic angles are the UVs [28] in Eq. ( 2): x , where x = √ aΔE for elliptical orbits and E is the eccentric anomaly, ẑ = − x2 , and the Stumpff functions [1,28], C(ẑ) and S(ẑ).

Defining a BVP Under Vinti's Potential as a Universal System of Equations
A congruent IVP and BVP are governed by the same dynamical system, which in this work is generated by Vinti's symmetric, integrable potential.Consequently, when posing a boundary value problem instead of an initial value problem under the Vinti potential, the BVP may be formulated in part by rearranging many of the equations developed for the analytical universal Vinti propagator [22] into a new system of equations.To derive this system of equations, the general problem is first defined.From there, a good starting point is to define the orbital plane, and then to examine the OS Lagrange coefficients developed by Biria [28].The section concludes with procedures for enabling multiple-revolution (multi-rev) solutions and obtaining terminal velocity vectors after the root-solve.The root-solve itself is addressed in a later section.

Problem Definition
In general, the BVP is stated as: Given terminal position vectors r 1 and r 2 , the desired time of flight Δt * , the direction of motion (DOM), and the number of com- plete revolutions N, find the terminal velocity vectors v 1 and v 2 .Under Keplerian dynamics, the DOM is traditionally translated to a binary variable dk ∈ {−1, 1} that divides the solution space into two branches according to the range of the transfer angle, Δf k .The binary variable indicates short-way (SW) solutions ( 0 ≤ Δf k ≤ ) or long-way (LW) solutions (  < Δf k < 2 ) with dk = 1 or dk = −1 , respectively Rus- sell [7], where the subscript k denotes Keplerian motion.Alternatively, the retrograde factor, K ∈ {−1, 1} , another binary variable, can define the DOM instead, where K = 1 for direct orbits and K = −1 for retrograde orbits, noting from Russell [7] that if a Lambert algorithm accepts K as an input, then dk must be determined in the first step.While the DOM is still a required input under Vinti's potential, it is translated in a different way that is explained in a later section.To identify multi-rev solutions, N is used as a parameter, where the signed integer Ñ ∈ ℤ described by Russell [7] is adopted in this work to distinguish by its sign between short-period (SP) and long-period (LP) solutions.Using the same convention, Ñ < 0 corresponds to an SP solution, Ñ > 0 to an LP solution, and N = | Ñ| .Note that Δt , determined from the OS time of flight equation, is distinct from the input Δt * .The difference between the computed and desired times of flight, Δt − Δt * , in the end is root-solved as part of a larger root-solve procedure.

Mapping Terminal Positions to Oblate Spheroidal Geometry
Since Vinti's potential is defined in OS coordinates, both position vector boundary conditions (BCs) must be mapped to Vinti's OS geometry.This map is simple, where Vinti's equations [18] apply for i = 1, 2 as with r i = ‖r i ‖ , and the equation for the OS position unit vectors from Biria [28] applies as The particular definition of the OS position vectors leading to the above unit vectors is directly connected to the OS equinoctial basis vectors [23].

OS Orbital Plane Orientation
In contrast to the two-body Lambert problem, the given position vectors do not define the orbital plane.In fact, because the plane rotates, the OS position vectors still do not define the OS orbital plane.In other words, ̂ 1 is in the orbital plane at t 1 (6) . and ̂ 2 is in the orbital plane at t 2 .Choosing ̂ 1 as common to both the ECI and rotat- ing frames, rotate ̂ 2 about the polar axis by the OS RAAN drift, ΔΩ � , according to so that both vectors lie in the orbital plane, where the secondary index "1" denotes t 1 and "2" denotes t 2 .Now that ̂ 2 is in the same rotating reference frame as ̂ 1 , and both vectors lie in the orbital plane, the transfer angle subtended by ̂ 1 and ̂ 21 is actually the change in true argument of OS latitude, Δ , which can be computed from the dot product and cross product to retain accuracy.Specifically, cos Δ and sin Δ are obtained as respectively, and, using Eq. ( 9), basic trigonometry determines the angle Δ as where d is the SW/LW indicator for Δ , distinct from df , the SW/LW indicator for Δf with the same direction of motion.The indicators are +1 for SW ( 0 ≤ Δ(⋅) ≤ ) and −1 for LW (  < Δ(⋅) < 2 ).It makes sense to define the direction of motion by specifying the orbit as direct or retrograde with K first, since that determines the direction of the node, and then derive d from K and the cross product of the OS terminal positions, which is not a unit vector, as The traditional transfer angle, the change in OS true anomaly ( Δf = ΔW 0 ), and df can be determined as The mutual OS inclination can be obtained from the OS angular momentum direction, ŵ1 , which is a unit vector, as 38 Page 10 of 41 Note that I is not explicitly needed, only sin I and cos I , so that computing the arc- tangent may be avoided:

OS Lagrange Coefficients
Biria [28] defined and derived OS Lagrange coefficients, F , G , ̇F , and ̇G , under Vinti's potential in terms of various differential anomalistic angles.Equating the expressions for these coefficients in terms of Δf and x gives where h i =  2 i ̇fi for i = 1, 2 and ̇fi , the time derivative of the OS true anomaly, is defined by Biria and Russell [23] as which is also repeated in Eq. (42) for convenience.Adapting Deyst's approach [1], several useful relationships can now be obtained.Solving Eq. ( 16) for x yields Careful manipulation of Eq. ( 18) leads to a more useful alternative form of the equation.To derive this alternative equation, recall from Biria [28] that the OS conic equation can be written in terms of Δf as ( Then, rearrange Eq. ( 22) as divide Eq. ( 23) by p 1 to obtain and subtract 1∕ 1 from both sides of Eq. ( 24) to yield Next, manipulate the LHS of Eq. ( 18), multiplying the entire LHS by but canceling that operation out by multiplying both terms inside the square brackets by the reciprocal of Eq. ( 26), yielding Bringing the factor 1∕ √ p 1 inside the square brackets, canceling terms, and apply- ing some trigonometric identities in Eq. (27) gives Finally, recognizing that the RHS of Eq. ( 25) appears in Eq. ( 28), substitute Eq. ( 25) into Eq.( 28) and substitute the result into the LHS of Eq. ( 18) to obtain the desired alternative form of Eq. ( 18) as After substituting Eq. ( 21) into Eq.( 29), canceling h 2 ∕p on both sides, multiplying both sides by 1 2 , and defining the auxiliary variables Ã and ỹ , adapted from Bate et al. [1], as and ( 23) the resulting equation can be simplified to a more compact form as Note that Ã should be computed differently to avoid the indeterminate form for small values of Δf as When 1 + cos Δf ≤ 10 −2 , Ã should be computed from Eq. ( 30), where the threshold of 10 −2 is recommended by Russell [7] to mitigate the loss of precision that arises when cos Δf ≈ −1 .Substituting Eq. ( 31) into Eq. ( 21 yields a compact form for x in terms of ỹ and ẑ as Finally, Eq. ( 17) can be related to the time of flight by observing and substituting Eq. (35) into Eq.( 2) to obtain Note that Eq. (36) requires Δ and Δf explicitly in the secular terms, computed from inverse trigonometric functions, unlike the two-body Lambert problem that requires only cos Δf , not the angle itself.
Looking back over the preceding analysis, Eqs. ( 36), ( 3), ( 4), ( 14), ( 10), ( 13), (32), and (34) can be viewed as a system of eight equations in the eight unknowns ẑ , Δ � , ΔΩ � , I, Δ , Δf , ỹ , and x , respectively.Due to the nature of the specific functional rela- tionships between the variables, it is not required to iterate on all eight variables.To see why, consider first the following simpler system of three equations in three unknowns governing the two-body Lambert problem: which Bate et al. [1] state as Eqs.5.3-1, 5.3-2, and 5.3-3 in their book.The three unknowns are, using classical orbital elements, p k , a k , and ΔE k or, in universal form, ỹk , xk , and ẑk , recalling that ẑk = ΔE 2 k for elliptical orbits.Bate et al. [1] point out that iterating on a k does not lead to a straight forward root-solve algorithm because guessing a k does not uniquely determine p k or ΔE k .But iterating on either p k or ẑk does lead to a straight forward root-solve algorithm because guessing p k or ẑk uniquely determines the remaining two variables.Iterating on p k leads to the so- called p-iteration method [1]: guess the first variable, p k ; the second variable, a k , is determined from p k and the problem geometry, Δf k ; the third variable, ΔE k in the elliptical case, is determined from p k and a k , and the time of flight equation is used to check the trial value of p k .Alternatively, Bate et al. [1] describe the following universal approach credited to John Deyst: guess the first variable, ẑk ; the second variable, ỹk , is determined from ẑk and the problem geometry, Δf k ; the third variable, xk , is determined from ẑk and ỹk , and the time of flight equation is used to check the trial value of ẑk .Interested readers can refer to Bate et al. [1] for the details, but the key takeaway is that, due to the nature of the functional relationships, it is not required to iterate on all three variables in the two-body Lambert problem.Rather, if the iteration variables are chosen carefully, then a solution method can be devised that iterates on a single variable to satisfy the time of flight equation.
A very similar argument is applied to the system of eight equations in eight unknowns governing the oblate Lambert problem.If the iteration variables are chosen carefully, then a solution method can be devised that iterates on only three variables to satisfy the OS time of flight, OS apsidal drift, and OS nodal drift equations, which refer to Eqs. (36), (3), and (4), respectively.First, guess the three iteration variables ẑ , Δ � , and ΔΩ � The fourth variable, I, is determined from Eq. ( 14), depending on ΔΩ � and the problem geometry.The fifth variable, Δ , is determined from Eq. ( 10), also depending on ΔΩ � and the problem geometry.The sixth variable, Δf , is determined from Eq. ( 13), depending on Δ � and Δ .The seventh variable, ỹ , is determined from Eq. ( 32), depending on ẑ and Δf .The eighth variable, x , is determined from Eq. (34), depending on ẑ and ỹ .It follows that guessing ẑ , Δ � , and ΔΩ � uniquely determines the remaining five variables.Finally, the OS time of flight equation is used to check the trial value of ẑ , the OS apsidal drift equation is used to check the trial value of Δ � , and the OS nodal drift equation is used to check the trial value of ΔΩ � .Additional details are given in Sect. 5.
To summarize, as the two-body Lambert problem can be reduced from three equations in three unknowns to one equation iterating on ẑk [1], the present sys- tem can be reduced to three equations, Eqs.(36), (3), and (4), iterating on the three 38 Page 14 of 41 unknowns ẑ , Δ � , and ΔΩ � .It is not known if the system of equations can be further reduced, although three equations appears to be the minimum.These three equations can be viewed as a generalization of Lambert's equation to include J 2 , in one sense by generalizing the time of flight equation, which is the tradi- tional Lambert equation, and in another sense by adding two equations to the dynamical system to account for the drift of the node and apse line.These additional two equations may be considered present but degenerate in the classical Lambert problem, where Δ � = ΔΩ � = 0 , and nondegenerate in the BVP defined by the Vinti potential.In fact, as J 2 → 0 , (Δ � , ΔΩ � ) → 0 in Eqs.(3)(4), respectively, and Eq.(36) reduces to the traditional Lambert time-of-flight equation, √ Δt = x3 S + Ã√ ỹ , which is plainly embedded as the first two terms in Eq. (36).A consideration of the fundamental frequencies of the respective dynamical models hints at the observed properties described above.The two-body problem is fully degenerate [30], possessing only the one frequency, and this property manifests as a single Lambert root-solve function in the traditional Lambert problem.Vinti's potential, generating what Wiesel [30] describes as a full spectrum of frequencies, defines a nondegenerate dynamical system, and this property manifests as three OS Lambert root-solve functions in the congruent BVP.Remarkably, all of these results are obtained without perturbation methods.

Adapting the Equations to Enable Multi-Rev Solutions
To enable N-rev solutions, add 2N to Δf , modifying Eq. ( 13) to yield If Δf tot is below or above the valid range for N, then add or subtract 2 , respectively, for both Δf tot and Δ tot .With this approach, the search for ẑ is bounded between 4N 2 2 at the left end (LE) and (2N + 2) 2 2 at the right end (RE) for a given N, and the kinematic equations must use the total angles in the secular terms, meaning 2N < Δf tot ≤ (2N + 2) .Put another way, the range requirement of Δf is absolute, while that of Δ is relative to Δf .

Obtaining Terminal Velocity Vectors After the Root-Solve
The OS Lagrange coefficients can be rewritten in terms of Ã and ỹ , noting ̇F is not required, as The coefficients are defined in a rotating frame, which means the frame's rotation rates are required [23] for i = 1, 2 as (40) where Eq. ( 43) was derived in Biria [22] and are given in many references [23,26].Alternative, singular forms of Eq. ( 43) are not recommended because they contain 1 − 2 i in the denominator, which goes to zero when the spacecraft is on a pole ( i = 1 ).The resulting indeterminate forms and catastrophic cancellation must be resolved with special mitigation techniques.For example, a conservative heuristic defines 1 − 2 i ≤ 10 −2 as the condition where i is too close to unity and the subtraction operation loses too much precision.Other conditions may be more appropriate depending on the application, Biria and Russell [23] recommending the alternative condition 1 − | i | < 10 −7 with a different goal in mind.In any case, when the user-defined condition is satisfied, 1 − 2 i may be obtained by computing 2  2 − 2 3 first as Then, noting that Eq. ( 45) can be substituted into Eq.( 46) to obtain where . Now, the simplest terminal velocity component anticipated is ̇z1 =  1 η1 + ρ1  1 , which requires η1 , obtained from Alternatively, Q cos 1 can be determined by computing the sign of Q cos 1 sepa- rately from its magnitude as 38 Page 16 of 41 which can be multiplied by Similar to the two-body Lambert problem, numerical difficulties appear for N transfers in Δ because the position vectors do not geometrically define the inclina- tion.An indeterminate form arises in Eq. (48), while in Eq. ( 50), the signs of the signum arguments could be inaccurate, and Q, obtained from Eq. ( 15), also loses accuracy from the cross product.
From the definitions in Biria [28], the usual formulas can actually still be used to obtain the terminal velocities, but they must be generalized for rotating frames like Eq. 29 in Biria [22].Following that logic, the OS relative velocities in the perifocal frame can be computed as where the superscript P denotes the OS perifocal frame, all other vectors are in ECI coordinates, Rodrigues' formula is used to obtain the rotation matrix R(Δ � ) from an axis-angle representation as the position vector 20 denotes 2 as viewed in the perifocal frame, computed as the notation I 3×3 denotes the identity matrix in Eq. ( 52), and the subscript × the skew symmetric matrix equivalent for the cross product expressed as The axis-angle representation is convenient here because the angular momentum vector is easily adopted as the rotation axis, which directly enables the use of Δ � as the rotation angle, notably defining the rotation in terms of physical quantities.Adding the centripetal term to Eq. (51) gives inertial OS velocities for i = 1, 2 as where the angular velocity vectors are given by (51 OS RAAN at the terminal points are determined as and differencing the rotation rates in Eq. ( 42) yields an expression for the instantaneous apsidal rotation rate, ω′ i , as Finally, by rearranging Eq. 27 in Biria [22], the terminal inertial velocities in ECI coordinates, v i , can be computed as Equations (51-58) draw useful connections to the two-body Lambert problem, but they also present two issues.The appearance of the singular element Ω � i in Eq. (56) implies numerical difficulties for nearly equatorial orbits and motivates an alternative approach based on equinoctial elements.More importantly, the division by G in Eq. ( 51) implies an inability to compute velocities when half-rev transfers in Δf are encountered, specifically an ambiguity in the velocity direction in the OS perifocal frame, analogous on some level to the traditional Lambert problem.But since the OS perifocal frame is rotating, intuition suggests that the singularity in Eq. ( 51) is non-physical and merely the result of the formulation or approach.To verify this assertion, begin from considerations of the OS orbital inclination, as could equivalently be done for the two-body case.The OS inclination is defined in Eq. ( 14), which traces back to the cross product ̂ 1 × ̂ 21 that appears in Eqs.(10) and (11).When ̂ 1 and ̂ 21 are parallel or anti-parallel, then the OS rotating orbital plane, and thus the OS inclination, is undefined.But under Vinti's potential, the angle subtended by these vectors is Δ , which suggests that Δ , not Δf , is the transfer angle associated with half-rev singularities.Under Keplerian dynamics, Δ � = 0 , and so the half-rev singularity becomes associated with Δf instead, viewed here through the lens of a degenerate dynamical system.When considering Vinti's potential, then, there appears to be no obvious physical singularity associated with half-rev transfers in Δf , thus verifying the earlier assertion.If this physical singularity is labeled a half-rev singularity, regardless of the dynamics, then the observation could be made that, in summary, the half-rev singularity in Δf seems to have disappeared under Vinti's potential, appearing in Δ instead.Note that the full-rev singularity exists (56 in both variables, in Δ because of the same cross product and in Δf because of the division by ỹ , which is evident later in Eqs. ( 62) and (63).A logical next step is to pursue a different formulation devoid of singularities when Δf = .Such an alterna- tive approach is proposed next to avoid both stated Ω � i and Δf -half-rev issues.To avoid computing Ω � i , the velocities can be obtained entirely in terms of OS equinoctial elements instead.From Biria and Russell [23], compute the following: f1 , ĝ1 , ̇f 1 , ̇ĝ 1 via Eqs.125-128, cos L 1 , sin L 1 from Eq. 90-91, and Li from Sect. 9. Compute ρi from Eq. 21 in Vinti [29] and ηi = ψi Q cos  i .To avoid inverse trigonometric functions, obtain cos L 2 , sin L 2 from standard angle sum identities after computing ΔL from Eq. (70) in Appendix A in [22].Then propa- gate [22] p 11 , p 21 via Eq. (71)in Appendix A to obtain p 12 , p 22 , which determine f2 , ĝ2 , ̇f 2 , ̇ĝ 2 as before.Finally, compute v i by applying to each terminal point Eqs.133 and 134 in Biria and Russell [23].See Appendix A for computational details.Notice that as J 2 → 0 , all equations reduce to the two-body equations.

Properties and Comparisons of the Lambert Solution Space with Oblateness Effects
The examination of the solution space, which should be independent of the rootsolve technique, is divided into two parts: zero-rev and multi-rev.The simpler zero-rev case is explored first.Once basic properties are established, attention is drawn to the multi-rev case, which adds a significant layer of complexity.

Visualization Strategy and Example Setup
Under the Vinti potential, many of the properties observed under two-body dynamics continue to hold.Properties are interpreted from the Vinti and Keplerian Δt-ẑ curves [1], such as in the top panels in Fig. 2,where each point on the curve corresponds to a feasible Lambert-type transfer orbit between a fixed r 1 and r 2 in the com- puted transfer time Δt .Short-way transfers are depicted on the left and long-way transfers on the right.These transfers are automatically feasible for the Keplerian dynamics because the root-solve is one-dimensional.However, it is important to emphasize their feasibility for Vinti dynamics.To obtain this curve for Vinti dynamics, the three-dimensional root-solve of Eqs. ( 36) and (3-4) is collapsed onto one dimension by ensuring that the root-solve of Eqs.(3-4) for Δ � and ΔΩ � is converged for each ẑ on the curve (see the next section for root-solve details).This extra step enables the direct comparisons between Vinti and Keplerian dynamics in Fig. 2 and helps visualize the solution space.The terminal position vectors for Fig.

Zero-Rev Solution Space
Now that the analysis framework and example is set up, the top panels in Fig. 2 can be discussed in more depth.In this example, SW transfers are direct while LW ones are retrograde.The SW transfers, being direct orbits, have an inclination close to the value in Table 1.Ignoring the multi-rev cases for now, note that the SW and LW solutions are distinct, possessing the same general properties as observed for two-body dynamics.For either case, and noting the logarithmic scale on the vertical axis, Δt → ∞ as ẑ → (2) 2 , a limiting parabola with finite p (the orbit approaches a different limiting parabola when the asymptote is approached from the other direction).As ẑ is decreased from this value, it obtains the value for the well-known minimum-energy ellipse, which can be obtained analytically for Keplerian dynamics [2].The minimum-energy spheroidal ellipse is marked with a blue square and the Keplerian ellipse with a red "x".The two energy values are distinct, occurring for different TOF and ẑ values, and only appearing to be identical at this scale.While the minimum energy, min , retains the same value for all N and directions of motion under Keplerian dynamics [7] ( k,min = constant), this property does not hold under Vinti dynamics because Δf , and thus the chord length in the perifocal frame, changes with ẑ , N, and K.As ẑ is decreased further, the curves intersect the vertical axis when ẑ = 0 , corresponding to the parabolic transfer, and ẑ then obtains negative values for a range of hyperbolic transfers.The curves for SW transfers steeply intersect the horizontal axis at some negative ẑ as Δt → 0 , corresponding to a straight line ( p = ∞ ) relative to the OS orbital plane, while Δt → 0 asymp- totically for LW hyperbolic transfers, theoretically equaling zero for the limiting rectilinear hyperbola ( p = 0) [1], although noting that the formulation completely breaks down well before that limit [26] with the forbidden zone [31] boundary sitting at 210 ⪅ ⪅ 297 km for Earth.At this scale, there is almost no discernible dif- ference between the dynamical models, except for the LW hyperbolic transfers in the top panel in Fig. 2b, the curves visibly approaching the horizontal asymptote at different rates.This discrepancy is attributed to periapsis approaching Earth's center (  p < 235.2 km and r p k < 249.7 km when ẑ = −52.008rad 2 ), so that the effect of J 2 increases as ẑ decreases.Zooming in would reveal more discrepancies, and, as shown in a later section, these differences only become more exaggerated as N increases because the implied longer portion of the transfer time spent near periapsis allows errors to accumulate.Alternatively, the difference between the two curves, denoted as Δt err , can be plotted as in the second row of panels in Fig. 2, which clearly quantifies the discrepancies.The TOF error appears to blow up near the ẑ-full-rev asymptotes.As ẑ is decreased, Δt err for the long-way solution goes negative as the curve approaches the asymptote, the traditional Lambert solution overestimating the required transfer TOF by ≈ 17.0 sec at the LE.For the short-way solution, the TOF error curve exhibits a hook-like characteristic toward the left end, going slightly negative in a small region between ≈ [−4.54, −0.49] rad 2 , and otherwise remaining positive.The traditional Lambert solution underestimates the required SW transfer TOF by ≈ 5.6 sec at the LE ( ẑ = −13.0656rad 2 ) near where the Δt-ẑ curve intersects the horizon- tal axis.Note that the LE has the largest p for which the Vinti-Lambert algorithm converges, so p is not infinite ( p ≈ 266, 469 km, Δt ≈ 92.2 sec).
The next four rows of panels in Fig. 2 illustrate some of the characteristics of the other six main root-solved quantities: two transfer angles, two orbital drift angles, the orbital inclination, and the periapsis radius.Because the root-solve of Eqs.(3-4) is converged, the values of Δt , Δf , Δ , ΔΩ � , Δ � , I, and p can be directly read off the plots for a given ẑ .Results in Fig. 2 confirm basic expectations, that the first five root-solved quantities vary under Vinti dynamics while remaining constant under Keplerian dynamics, for which Δ k = Δf k is constant because ΔΩ k = Δ k = 0 and i k is constant because the orbital plane is inertially invariant. 1Periapsis, denoted p (spheroidal) or r p k (spherical), behaves similarly under either dynamical model.

Page 22 of 41
What immediately stands out is that the long-way Vinti transfers exhibit strong deviations away from the Keplerian transfers as the orbit becomes more hyperbolic.Evidently, because these transfers pass through periapsis and the periapsis radius is small, approaching 235 km at the left endpoint, J 2 considerably warps the inclina- tion, which shoots up 15 • to become more equatorial, and significantly contributes to the turning angle.Since the Vinti-Lambert solver was not intended to work for such a low periapsis radius, a few of the very hyperbolic transfers were more closely examined to assess whether the solutions were even remotely accurate or representative of the intended dynamics.Comparisons were made between the associated universal Vinti propagator and numerically integrated Vinti trajectories, and results indicated that the substantial increase in orbital drift and inclination still faithfully capture the Vinti dynamics to sub-kilometer accuracy.With a spacecraft on such trajectories traveling through a minimum radial distance close to 235 km, these scenarios are not practical, but the sub-kilometer level of accuracy is still considered sufficient to at least qualitatively characterize some of the corners of the solution space.This exercise is viewed more as an important test of the limitations of the formulation and the algorithm, and also a testament to the formulation's integrity.

Multi-Rev Solution Space
Figure 2 also reveals some preliminary, novel characteristics of multi-rev scenarios that appear even for very low N, limited to 1 or 2 revs in this data set.Qualitatively, curves for the Vinti transfer angles, drift angles, and inclination are observed to have either an arctangent ( arctan ) or arccotangent ( arccot ) shape for each N, appearing to approach distinct horizontal asymptotes as ẑ approaches either vertical asymp- tote.While the Lambert solutions approach the same limiting parabolas at the ẑ asymptotes for each N under Keplerian dynamics, J 2 causes the limiting parabolas to change slightly for each N (different p, I).When switching from the SW to LW solution space, the I curve, having the observed property I > i k for all examples considered, changes with increasing N from increasingly inclined arctan shapes to progressively more equatorial arccot shapes.Similarly, the maximum Δ and Δ � increases with N, switching from arctan shapes (SW case) to arccot shapes (LW case).The minimum Δf decreases with N, switching from arccot shapes (SW case) to arctan shapes (LW case), while the maximum |ΔΩ � | increases with N, maintaining an arccot shape regardless of the solution family but changing from negative to posi- tive when switching from direct to retrograde solutions.
Having covered some basics of the low-rev cases, the discussion is now focused on some general properties of multi-rev scenarios, using the low-rev results as illustrative examples.The multi-rev case adds a significant layer of complexity and poses interesting consequences to transfer design, particularly with regard to the existence of solutions and repercussions of incorrectly determining existence.An important property of the traditional Lambert solver is its ability to determine the existence of multi-rev transfers by identifying the minimum TOF, Δt min .As illustrated in the zero-rev case, the required transfer TOF does not agree in general between Vinti and Keplerian dynamics, and this discrepancy extends to Δt min for the multi-rev case ( ΔTOF ≡ Δt min − Δt k,min ≠ 0 ), where the traditional Lambert solver generally over- estimates or underestimates Δt min for an N-rev transfer, corresponding to a potential false negative or false positive, respectively, in terms of the existence of multi-rev solutions given Δt * .Existence can be determined from the Δt-ẑ curves or, alterna- tively, a rough idea may be obtained for low revs by referencing the second row of panels in Fig. 2. While ostensibly there is no difference in the TOF error trends in Fig. 2, a positive or negative TOF error at ẑΔt min , indicated by the red + symbols, can be interpreted as a false positive or false negative, respectively.For the short-way transfers, Δt err (ẑ Δt min ) < 0 for N = 1, 2 , the error being larger or more negative for N = 2 at ≈ −17.84 sec versus ≈ −4.58 sec for N = 1 .The long-way transfers exhibit the opposite trend: Δt err (ẑ Δt min ) > 0 for N = 1, 2 , the error being larger or more positive for N = 2 at ≈ 33.55 sec versus ≈ 14.03 sec for N = 1 .If obtained directly from the Δt-ẑ curves, the minimum TOF for short-way transfers is overestimated by 17.42 sec for N = 2 and 4.58 sec for N = 1 (false negative), while for long-way transfers it is underestimated by 33.93 sec for N = 2 and 14.03 sec for N = 1 (false positive).In this example, J 2 is thus observed to decrease Δt min for direct transfers and increase it for retrograde transfers, corresponding to a respective minimum TOF advance or delay, apparently depending on the orbital regime.It seems reasonable to suspect that the trend observed is directly caused by the orbital inclination, acting through the well-known apsidal and nodal drift phenomena and the secular effect on the mean anomaly, in addition to some interplay with periapsis distance and the length of time spent closer to periapsis.That is not to say that direct transfers always experience a minimum TOF advance and retrograde ones experience a delay; rather, the actual trends are thought to be much more nuanced.For instance, in this example, since the inclination is roughly 30 deg away from being equatorial in either case, the apsidal drift should impact the transfer time to a similar degree and in the same "direction" ( Δ � > 0 ) for either inclination, and therefore cannot be the dominant cause of the observed trend.By process of elimination, some combination of nodal drift and periapsis considerations are thought to be causing the trend, but it is not clear how.These trends hint at underlying physical mechanisms affecting the minimum TOF and are explored further in the next section.

Minimum TOF Advance or Delay with Practical Considerations
Preliminary notions of the physical mechanisms controlling trends in the existence of multi-rev solutions were probed in the previous section.It is shown in the following that the impact of accumulated errors from neglecting oblateness can be considerably more profound, even at Earth, and the impact on transfer existence exemplifies in different ways how neglecting J 2 can lead to errone- ous conclusions.Multiple examples are examined to investigate more broadly, though not exhaustively, how inclination can affect the solution space.

Very-High-Rev Example
An illustration of this problem for a 100-rev transfer is offered in Fig. 3 using a log scale as in Fig. 2. BCs are generated from Tables 1-2, leading to approximate terminal positions r 1 = [5558.706,3761.031,1990.162]⊤ km at one end and r 2 = [−3077.908,−514.361,−6285.426]⊤ km at the other.With r 1 and r 2 fixed, the Δt-ẑ slice of the solution space in Fig. 3a appears to retain the same basic shape with J 2 effects, but now with a more pronounced discrepancy.It is evi- dent from a glance that the two Lambert solvers predict a significantly different minimum TOF: the Keplerian Lambert solver indicates a minimum of 6.32 days while the more accurate Vinti one indicates a 6.66-day minimum, which is a difference of ≈ 8.1 hours.The minimum-time 100-rev short-way Vinti transfer has the following approximate orbital elements: p = 6173 km, a = 6916 km, e = 0.11 , I = 80.75 • , ΔE = 166.0• , Δf = 153.8• , Δ = 131.6 • , ΔΩ � = −8.2• , and Δ � = −22.2• .Since the transfer is only ≈ 9.25 deg away from being polar, the nodal drift should not significantly affect the transfer time, even for 100 revs.In fact, the plot clearly shows a time delay, implying that the apsidal drift ( Δ � < 0 ) is not only causing the delay, but is also great enough to overcome other effects.Ignoring other effects, the apsidal drift induces a time delay because, relative to the Keplerian solution, the transfer angle grows when Δ � < 0 , or roughly when 63.4 • < I < 116.6 • , and shrinks when Δ � > 0 , or roughly when I < 63.4 • or I > 116.6 • .It is straightforward to connect apsidal drift to a change in transfer angle because the drift is in the OS orbital plane, leading to the simple angle summation in Eq. ( 13).In contrast, ΔΩ � is not easily connected to a change in transfer angle, operating through Eqs.(8)(9)(10)(11)(12).In the present case, Δf k = 130.0 • and ΔE k = 154.1 • , the latter being cubed in Eq. ( 36) and larger with J 2 effects, leading to the longer TOF and 8.1-hour discrepancy.
Figure 3b highlights another interesting feature of Vinti dynamics that does not exist in the two-body case.As stated earlier, all transfers for the same BCs share the same minimum energy under Keplerian dynamics, but this property vanishes with J 2 effects.In this example, a 100-rev transfer not only requires a lot more time, but also much more energy ( Δ min ≈ 1.02 km 2 /sec 2 ) to execute compared to the two-body prediction, the Vinti energy curve appearing to reside inside a Keplerian energy envelope.The inverse was found to occur for BCs leading to a Δt min advance rather than a delay, where the Keplerian energy curve resides within a Vinti energy envelope and the Vinti transfer requires less energy ( Δ min < 0 ).Note that the envelope findings apply to the left end of the curves, not necessarily generalizing to the entire energy curve that includes regions where Δt → ∞ .In these regions, the inner curve was observed to pierce the "envelope", and it is not clear without further analysis whether this is an inaccuracy resulting from the third-order Vinti approximation, numerical inaccuracy due to large TOF, a different source of inaccuracy, or a property of the exact Vinti dynamics.These observations are based on a few numerical examples and warrant further study.

Multi-Rev Inclination Sweep
While the preceding results are illuminating, they still fall short of a broader mapping and characterization of the Lambert solution space with J 2 , which is better achieved with the inclination sweep presented in Fig. 4.These results emerge from a sweep across the domain of direct, short-way transfers with N = 20 revs, covering eight separate BVPs or BC pairs generated from the values in Tables 1 and 2. Various methods would be considered appropriate ways to choose these BVPs, and the specific thought process employed is discussed here for interested readers.Choosing 20 complete revs is somewhat arbitrary, but with a main goal of ensuring that inclination-induced J 2 effects would be visible in the results.Upon finding that such J 2 effects are easily observed with 20 revs, it follows that, for this type of investigation, a sweep at a moderate N like 20 revs is preferred to the very high N of the 100-rev example shown in Fig. 3.
With N selected, the next decision concerns how to constrain the geometry.Since the Keplerian solution space is being compared to the Vinti solution space, it would help, if possible, if the Keplerian solution space were constant while varying the inclination.For example, it seems that an appropriate approach could keep the in-plane geometry constant, under Keplerian dynamics, for the entire parameter sweep, to isolate the effects of varying I.However, there is a secondary goal, which is to demonstrate that, while the transfer angle is a simple function of the fraction of the orbit traversed when under Keplerian dynamics (and independent of I), this property does not hold under Vinti dynamics.Specifically, the inclination affects how long it takes to traverse through a certain orbit fraction, in addition to how long it takes to complete a certain number of revs (i.e., the orbital period varies).Now, a good way to accomplish both goals stated above is to use the same procedure to generate the BCs for Figs. 2 and 3, and then hold r 1 and the number of revs fixed for all cases while varying I.With some trial and error, choosing 20.4 revs led to favorable results that exposed interesting features of the solution space.The choice kept Δf k constant enough for the purposes of the first stated goal, spanning only ≈ 9 • , which is considered an insignificant change in the in-plane orbital geometry.Secondly, the choice demonstrates that holding the ICs and number of revs fixed at 20.4 revs over varying inclinations requires varying times of flight, and also results in varying transfer angles.In summary, an attempt is made to keep Δf k nearly the same for each case in Fig. 4, though this consistency in the in-plane geometry is balanced against a desire for Δt 0 to propagate each IC consistently through 20.4 revs under Vinti dynamics.For brevity, the sweep does not extend into the retrograde domain, which would include an additional eight BVPs, nor does it include the other half of the solution space for the selected BVPs, the long-way counterparts excluded from Fig. 4 that would also happen to be retrograde.Now that the BVPs are selected, the results in Fig. 4 are explored and discussed.Each column or subfigure of Fig. 4 is organized like Fig. 2 but without the periapsis panes, and Figure 4 as a whole can be analyzed by stepping from the left to right subfigures, i k increasing from ≈ 12 • to exactly 90 • .While some behaviors are main- tained from the low-rev examples in Fig. 2, others vary considerably.As i k is walked up to the polar case, the arctan behavior in I is retained, but notice the variations in I decreasing from a spread of ≈ 22 • in Fig. 4a (LE: p ≈ 4984 km; RE: p ≈ 1762 km) down to 0 • in Fig. 4h (LE: p ≈ 4524 km; RE: p ≈ 2691 km).The spread varia- tion is not only caused by the greater strength of J 2 close to the equator, but also its accumulated effect over many revs.The concentration of mass closer to the equator and away from the poles can be inferred from Fig. 4 in other ways, where the effect on Δf over 20 revs is markedly unintuitive.Specifically, while Δf is observed in Fig. 2 to vary monotonically with ẑ for low revs, Δf does not necessarily vary monotonically with ẑ for high revs.In Fig. 4a, b, for example, an arccot shape in the trend for Δf that may have been anticipated from Fig. 2 instead appears as if slightly corrupted by J 2 .In contrast, the trend in Δ retains the arctan shape seen in Fig. 2, the spread in Δ decreasing, as observed for I, from ≈ 24 • in Fig. 4a down to 0 • in Fig. 4h.The constant Δ observed for the polar case in Fig. 4h is not an intuitive result, and it is viewed as an instructive special case warranting further explanation.First, consider that the OS inclination is expected to be constant in the polar case because a polar transfer orbit is only possible with a 90 • inclination.To understand why Δ is constant in the polar transfer solution space, consider that ΔΩ � = 0 for all ẑ , leading to the condition ̂ 21 = ̂ 22 .It follows that Δ , the angle subtended by ̂ 1 and ̂ 21 , is equal to the angle subtended by ̂ 1 and ̂ 2 , analogous to how Δf k is the 38 Page 30 of 41 angle subtended by r 1 and r 2 under Keplerian dynamics.For polar Vinti transfers, then, the Δ transfer angle is solely dependent on the BCs or supplied geometry in the same way that Δf k is in general under Keplerian dynamics.A remarkable condi- tion follows: Δ ≈ Δf k for all ẑ .Of course, Δf still varies with ẑ at the polar inclina- tion because Δ � ≠ 0 .At i k ≈ 38 • (Fig. 4c), Δf is seen to vary the least, by ≈ 3 • .
For orbit drift angles, ΔΩ � retains an arccot shape that decreases in magnitude and spread as i k increases.In contrast, Δ � is observed to have a sort of inflection point (IP) captured in Fig. 4d.It does not coincide with a transition from positive to negative drift, but does seem to coincide with when ΔTOF is almost minimized, suggesting that some balance of J 2 effects occurs in this regime that may be con- nected to Δf obtaining an arctan shape after the IP.The apsidal drift retains an arctan shape before the IP (with a noticeable maximum almost coinciding with the minimum Δf in Fig. 4a) and an arccot shape after.Within ≈ 5 • of the IP, it is not known if the symmetric binomial shape with positive drift seen here is common for Vinti transfers in this regime.A deeper understanding of the IP may be obtained with further study.Equatorial trends were omitted for brevity.
The hypothesis that the inclination may serve as a major physical lever contributing to a minimum TOF advance or delay is further supported in Fig. 4, where strong correlations are observed between I and ΔTOF , with values given in each subfigure caption.A Δt min advance is observed for low i k , reaching nearly 58 min in Fig. 4a, transitioning through a minimum Δt min discrepancy near the i k 0 = 55 • case with ΔTOF ≈ 2 min, and then growing to a Δt min delay of ≈ 18 min for the polar case.It is suspected that ΔTOF = 0 at some i k between 38.4 • and 48.6 • , probably near 46 • , though no attempt is made to find this exact inclination where the estimate of Δt min agrees between the two dynamical models.Evidently, the potential for false negatives and false positives is substantial and one situation does not appear to be more common than the other.

Practical Considerations
Turning to some practical considerations, suppose, now referring to Fig. 3, that a 6.5day transfer is desired to satisfy mission requirements, indicated by the green horizontal dashed line.Here, a traditional Lambert solver erroneously claims a 100-rev transfer exists, a false positive.Without the insight of the Vinti-Lambert solver, it would be tempting to use the 6.5-day solutions as guesses in a perturbed Lambert algorithm, but Fig. 3 clearly shows that would be pointless.In fact, for a typical approach, such as a shooting method, the algorithm would waste valuable time trying to find a solution that does not exist, and the algorithm would ultimately fail for an unknown reason.These unexpected, undiagnosable failures can occur whenever ΔTOF > 0 , observable in Figs. 3 and 4. Turning to Fig. 4a, suppose a 31-hour transfer is desired.Here, a standard Lambert solver erroneously claims a 20-rev transfer does not exist, a false negative, and does not return a 20-rev solution to the shooting method, noting that N k,max may be several revs lower and the 20-rev case may not even be attempted in an algorithm that first identifies N k,max .Practically speaking, this situation is not catastrophic, but, when performing broad searches, for example, it does imply a statistically significant probability of overlooking good many-rev solutions that may save a large amount of fuel or be otherwise favorable.
The false negative scenario underscores the importance of developing good alternatives to using a traditional Lambert solution as a guess, because the two-body Lambert guess can literally skip over useful parts of the solution space without giving any clue to the user that much better transfers may exist.That being said, the results of the present study do point to potentially useful mitigation techniques, such as heuristics based on the inclination, that do not even require the practitioner to implement a Vinti-Lambert solver.Heuristics for the whole solution space would require more study, but something could be said for direct short-way transfers supported by Fig. 4. For example, in multirev cases, if an SM reaches a maximum number of iterations and i k > 45 • , the user can be informed that a transfer likely does not exist.If i k < 45 • and the two-body Lambert guess returns N k,max , the user can be warned that an alternative guess is recommended, if not executed, to determine if solutions are being overlooked for N > N k,max .These techniques can be refined and improved with further investigation, noting the i k thresh- old, which could be more of a region than a hard boundary, may depend on {c, N, Δf }.
The analytical BVP solver presented in the current work, which operates under Vinti dynamics, is seen in Figs. 3 and 4 to offer new insights that directly inform mission design, a way for practitioners to determine the existence of more realistic transfers and weed out infeasible ones.In simple diagrams, it also attributes to a physical cause some of the difficulties encountered in applications like ASDR, and it motivates, more generally, some of the risks assumed when neglecting J 2 in preliminary analyses.

A Method for Solving the BVP Defined by the Vinti Potential
Having converted the BVP to a system of equations and explored major parts of the solution space, an elementary method is now proposed and employed to solve that system of equations.While the system of equations can be solved in multiple ways, only one method is explored in the current work.

Initial Guesses for the Unknowns
There are multiple ways to obtain an initial guess, but one of the advantages of adapting Deyst's algorithm to the Vinti BVP is that the Keplerian orbital elements (KOEs) from a robust two-body BVP solver can be used directly to generate a guess.This approach stands distinctly apart from any existing initial guess techniques in perturbed Lambert solvers, which typically ingest the two-body position and velocity instead of the orbital elements, working in the highly nonlinear ECI space instead of the more linear orbital element space.As such, KOEs offer a decent estimate of the secular drift rates in Ω � and ′ if e k < 1 as (60) 38 Page 32 of 41 with zero rates otherwise.Multiplying these rates by Δt yields decent initial guesses for ΔΩ � and Δ � as The maximum of either i determines I min as arcsin i,max , so i k may be replaced with I min if i k < I min .For multi-rev cases, further improve the guess by computing I from Eqs. ( 8) and ( 14) and iterating a few times on all four equations, replacing i k with I in Eq. ( 60) and computing Δ � 0 and its rate outside the loop.When guessing ẑ , simply set ẑ0 = ẑk , where {ẑ, ẑk } have the same range for N > 0 (not the same minimum TOF) and the same upper bound with otherwise similar trends for N = 0 .While notably overcoming one of the main drawbacks of the universal approach, which lacks a good guess for ẑk , this method warrants some words of caution.For the zero-rev case, this guess tends to avoid the issue of guessing a value of ẑ that places p inside the Vinti forbidden zone [31], which would require a step rejection and some way of handling it, a new issue that does not exist under Keplerian dynamics.On the other hand, the effect of J 2 is very small over zero revs, so it is expected that ẑ ≈ ẑk when converged, implying a large reduction in the number of iterations and a significant speedup.Checks on the minimum value ẑmin for the short-way solution must be retained [1], ensuring that ỹ ≥ 0 , but note that ẑmin will generally be different between the dynamical models because nearly rectilinear Vinti trajectories passing through the forbidden zone have an unknown analytical representation, leading to the stricter condition  p > c .In practice, robustness was observed to remain high when  p > 2, 500 km.Note, too, that the value of Ã is updated on each iteration, so ẑmin can fluctuate: the same ẑ value rejected on one iteration could be valid on the next iteration if Ã is sufficiently different.For multi-rev scenarios, it seems possible to guess ẑ in the wrong bin, meaning that the user may request an LP solution, but inadvertently seed the guess with an SP ẑ value that is outside the bounds on ẑ .In this exam- ple, the value of ẑk for the LP solution, which is to the left of ẑk (Δt k,min ) , lies to the right of ẑ(Δt min ) , because the ẑ values at the minimum TOF can differ greatly between the two dynamical models.One mitigation technique in this scenario would be to reset ẑ0 to its value at Δt min within some tolerance.
Nonexistence of two-body solutions does not preclude the existence of Vinti or perturbed counterparts, as was demonstrated in an earlier section.If the two-body Lambert solver fails to find a solution for a desired N, even if the time of flight is below the minimum, a robust Vinti-Lambert solver must still pursue a solution with a different initial guess.Armellin et al. [11] noted the benefit of avoiding a Keplerian initial guess in some circumstances in favor of an alternative technique, but while their approach may be useful in this situation, this alternative is not explored in the current study and is left to future work.

An Elementary Vinti-Lambert Algorithm
After normalizing by r 1 + r 2 to obtain canonical units [7] with = 1 , the algorithm begins with the evaluation of Eqs. ( 6) and (7) for  i ,  i , ̂ i , where repeated constants (61) are also stored.Next, from the previous section, compute initial guesses for ΔΩ � , Δ � from Eqs. (60-61) after a single Keplerian Lambert call (or an alternative if the solver fails), noting that all variables have an implicit subscript j to denote iterates that is generally dropped to reduce notational clutter.Then, compute Δ , Δf , cos I , Q = sin I from Eqs. (8-15) and Ã from Eq. ( 33) or (30).The previous steps are all done before guessing ẑ .If N = 0 , guess ẑ = ẑk ; otherwise, N > 0 and the bot- tom of the Δt-ẑ curve must be found first.Once found, guess ẑ = 0.9ẑ(Δt min ) for the LP solution or ẑ = 1.1ẑ(Δt min ) for the SP solution.With all initial guesses in place, the root-solve begins.For N = 0 , a 1D Newton-Raphson method iterating on ẑ is used to root-solve Eq. ( 36) with analytical partial derivatives of the Keplerian Lambert problem [1], considered a good enough approximation; deriving partial derivatives of the Vinti-Lambert problem is considered outside the scope of this work.
As implemented, a truncated version of Eq. ( 36), Δt = (x 3 S + Ã√ ỹ)∕ √  1 , is used until Δẑ j < 1 to help with convergence (derivatives including √ 1 as well), after which Eq. ( 36) is used in its entirety, the derivatives subsequently including one extra term, A 1 ∕ √  1 ⋅ dx∕dẑ , one of the major contributors to O(J 2 ) terms.Before Eq. ( 36) and its approximate derivatives can be evaluated, though, a root-solve of Eqs.(3)(4) is required, performed as a nested root-solve via the method of successive approximations (MSA) [32] while notingthat a 2D MSA root-solve at this step could be a viable alternative.First, compute C(ẑ), S(ẑ) , followed by Eq. ( 32) and Eq.(34) for ỹ, x , respectively.If Δẑ j > 1 , the truncated Eq. ( 36) can be evaluated at this point, but typically after one iteration, Δẑ j < 1 , which requires many additional computa- tions.In this case, two nested loops are used, one checking convergence of Δ � , the inner one checking ΔΩ � , implicitly requiring that Δf , and subsequently Ã , be updated on each iteration.Specifically, with this setup, ỹ, x must be updated at the beginning of the Δ � loop.Then, compute At this point, check the forbidden zone condition and reject the ẑ step if violated.If not, from Biria and Russell [23], obtain A 1 , B 1 , 1 , p 0 , S 1 , and Q 1 = C 2 from Eqs. 143-148 and i = { , √ p 0 , h z } from Step 4 (p. 287).Next, follow steps in Biria's UV Vinti propagator [22] and ultimately ΔR 1 j , ΔR 2 j , ΔR 3 j , ΔN 1 j , ΔN 4 j , which enables the determination of Δ � from Eq. ( 3) and ΔΩ � from Eq. ( 4).Relative to thepropagation steps, a key difference in the Vinti-Lambert solver is that the initial Q cos 1 must be obtained from Eq. ( 48) and 1 , e sin f 1 , and e cos f 1 from where h 1 =  2  1 ̇f1 and Eqs. ( 42) and (17) give ̇f1 and G , respectively.Once the updated value of ΔΩ � is calculated, Δ, d are updated from Eqs. (8-12).Careful atten- tion must be paid to quadrant ambiguities.If d = −1 , then Δ = 2 + Δ so that 0 ≤ Δ < 2 initially.Then, compute Δf , df from Eq. ( 13), but adjusted so that (62) 38 Page 34 of 41 if Δf < 0 , 2 is added to Δf , Δ , and if Δf > 2 , 2 is subtracted from Δf , Δ , maintaining the correct relative difference between the two angles while ensuring 0 ≤ Δf < 2 .For multi-rev cases, apply Eq. (40) to compute the total angular dis- placements Δf tot , Δ tot .With the updated Δ , compute Q, cos I from Eq. ( 15).If ΔΩ � has not converged to within some tolerance, then return to the computation of A 1 , B 1 and repeat.Otherwise, exit the inner loop, compute Ã , and check Δ � against its value on the previous iteration.If Δ � has not converged to within some tolerance,  return to the beginning of the loop, compute ỹ, x , Eq. (62), and continue through the MSA loops.Otherwise, exit the loop, finally evaluate Eq. (36), and update ẑ .This entire process is repeated until ẑ converges to within some tolerance.
For N > 0 , the minimization of the Δt-ẑ curve is solved with a golden section search on ẑ , followed by a root-solve of Eq. (36) via the bisection method.Gradient- free root-solve methods offer both simplicity and robustness in this part of the algorithm.Future implementations employing more efficient root-solve algorithms and better initial guesses are expected to enjoy significant performance boosts.Regardless of N, after the root-solve converges, the algorithm concludes with the evaluation of Eqs.(41-50) and the equinoctial element approach detailed in Appendix A to obtain the terminal velocities.

Performance Evaluation
Performance of the above algorithm is evaluated for robustness, accuracy, and speed in Fig. 5.All measures of performance consider all combinations of SW/LW zerorev transfers and direct/retrograde orbits for each eccentricity regime (ER), which is a specific range of e k .Since 500,000 BC-Δt * pairs are run for each of 4 direction/ inclination combinations, and there are 5 ERs, each performance metric is measured from 10 million scenarios, constructed as follows.The ER number on the horizontal axis refers to the uniformly sampled ERs given in Table 3, chosen to broadly test and stress the algorithm, in order as: broad elliptical, e k → 1 from below, e k near unity, broad hyperbolic, and limited hyperbolic, the last one depending on direction.For each ER, the ICs in the last line of Table 1 are transformed to ECI coordinates to obtain r 1 , where the square brackets indicate uniformly sampled quantities and f k 0 is computed differently according to e k .If e k < 1 , then All performance metrics are then measured from this data set.Robustness is interpreted in Fig. 5a as the percentage of Lambert calls that converge when expected, so that the identification of problematic scenarios (full-rev, low r p , etc.) are treated as successes and do not reduce the robustness measure.Robustness is observed to be high, ranging between 99.22 and 99.99% over all ERs, averaging 99.76% overall.In Fig. 5b, c, accuracy is also compared side-by-side with a Keplerian UV Lambert solver as a benchmark, where average miss distance is used as a proxy for the accuracy of v 1 , and v 2 is assumed to have similar accuracy.Miss distance is defined here as the difference between the BC r 2 , taken as truth, and the predicted value of r 2 resulting from the analytical UV Vinti propagation of the initial state , where v 1 is the output of either the Vinti or Keplerian Lambert solver.Figure 5c considers all cases, including unreasonably long TOF, and is only presented for completeness; Fig. 5b is considered a more useful measure of accuracy that only includes cases with Δt * < 1 year, which, while still a large TOF limit, has the benefit of still stressing the algorithm.As such, Vinti-Lambert accuracy is found to be remarkably high, SW transfers boasting an average miss distance of ≈ 44 mm across all ERs.LW transfers exhibit lower accuracy, averaging ≈ 10 km miss distance across all ERs and 7.5 m ignoring ERs 2-3.This LW trend likely stems from the increased likelihood that transfers will pass through periapsis, noting that they must pass through periapsis for hyperbolic transfers, with nearly parabolic transfers appearing particularly stressing.Relative to a standard Keplerian Lambert solver, the Vinti-Lambert solver improves the accuracy by as much as 9 orders of magnitude on average in some eccentricity regimes.In ERs 2-3, the discrepancy may be caused by qualitative disagreement between physical models, one predicting an escape trajectory where the other predicts a bounded one.Note that because the integrals are only evaluated to O(J 3  2 ) , the entire Vinti transfer orbit solution is of that accuracy, including the orbital elements.This property stands in contrast to analytical solutions to the IVP, for which half of the orbital elements, p, e, and I, can be obtained exactly (to double-precision accuracy) despite the lower accuracy of the integral approximations.
Computational speed is evaluated on an HP EliteBook 830 G5 laptop computer with an Intel Core i5-8350U CPU operating at a base speed of 1.9 GHz with 16 GB of RAM.The maximum speed at which a single one of these cores is capable of operating is 3.6 GHz, and the CPU operated at this maximum speed during tests, where the UV Vinti-Lambert algorithm is benchmarked in runtime against a UV Keplerian Lambert algorithm.Both algorithms are implemented in Fortran 90 and compiled with Intel Fortran Compiler Classic for Windows in 64-bit mode with O3 and Qip optimization settings enabled.For the described BC-Δt * pairs, the Vinti-Lambert algorithm implemented in this study only takes 18.0 s on average to compute v 1 , v 2 .The Vinti-Lambert solver is found to be only 15.4 times slower than the Keplerian Lambert solver, which takes an average of 1.2 s on the same machine.Considering a typical universal Lambert call runs at ≈ 1 s and the Vinti-Lambert solver is running about 15 times slower, a two-body initial guess is considered a wise use of computational resources as the performance gains are thought to outweigh the cost of assigning, on average, ≈ 6.5-7% of the runtime to calculating the initial guess.Multi-rev runtimes were not assessed as rigorously, but a Lambert call for the 100-rev scenario in Fig. 3a (Tables 1 and 2) was measured to take ≈ 4 ms (LP), with runtimes generally seen to be on the order of milliseconds for any N. Based on published performance data, the Vinti-Lambert zero-rev algorithm is estimated to run at least 1,000 times faster than comparable perturbed Lambert solvers [11,14], depending on the method, though the other algorithms were not independently implemented.Note that implementing more efficient root-solve algorithms than those presented in this work, which prioritized simplicity and robustness, is likely to yield significant additional Vinti-Lambert speedups.Such improvements to the oblate Lambert algorithm are expected to lead to multi-rev runtimes on the order of zero-rev runtimes, with the only performance losses resulting from the need to obtain both LP and SP solutions, as well as the formal minimization step to find ẑ(Δt min ) .For additional perspective, note that the Vinti-Lambert solver imple- mented in this work for zero revs is only 4.3 times slower than the author's UV Vinti propagator [22], whose runtime is estimated to be near the best achievable Vinti-Lambert runtime.

Conclusions
The classical Lambert problem is generalized to include the physical effects of a celestial body's equatorial bulge, generalizing Lambert's equation to a system of three equations that together define a new boundary value problem (BVP).This novel BVP is governed by Vinti's symmetric gravitational potential, which causes the problem to remain unperturbed despite the inclusion of J 2 effects, enabling the reduction to a system of equations via arguments and the oblate spheroidal (OS) Lagrange coefficients.Geometric connections are drawn between the orbital reference frames of the two BVPs to show how generalizing the classical Lambert inertial frame to a rotating frame enables, in an approximate sense, the definition of the Vinti-Lambert problem.Like the recent third-order universal solution to Vinti's initial value problem, the O(J 3  2 ) Vinti-Lambert formulation leverages OS equinoctial elements, OS universal variables (UVs), and the differential form of the equations of motion, collectively working to avoid indeterminate forms, angle ambiguities, and singular elements for all orbital regimes excluding a region around the forbidden zone.The Vinti-Lambert solution space is visualized, discussed, and extensively mapped for the first time for both zero-and multi-rev cases, and then directly compared to the classical Lambert solution space, leading to multiple new insights that distinguish this BVP from the classical Lambert problem.In particular, the advance or delay of the multi-rev minimum time of flight is shown to be a key contributor to shooting method pitfalls, while at the same time pointing to how the pitfalls can be mitigated.Some preliminary understanding of the effects of J 2 on minimum energy are observed as well, the investigation generally leaving ample room for future solution space exploration, including the study of corner cases and bifurcations.
As part of this study, an elementary, iterative method is presented for obtaining a solution to the zero-rev and multi-rev oblate BVP.Performance is rigorously evaluated at Earth for the zero-rev case over a broad range of orbital regimes, demonstrating, on average, robustness of 99.8%, accuracy boosts between three and nine orders of magnitude relative to Keplerian dynamics for times of flight less than 1 year, and computational efficiency with runtimes around 18 s, only 15.4 times slower than a 38 Page 38 of 41 classical UV Lambert solver.These high performance levels stand to benefit numerous applications.More extensive testing of the multi-rev case is left to future work, though multi-rev runtimes were consistently observed to peak on the order of milliseconds for cases with any number of revs examined in this study.While some areas of potential improvement have been identified, this study presents a novel algorithm for the solution of the oblate Lambert problem demonstrating both efficiency and robustness in comparison to existing methods in the literature.

Appendix A: Equinoctial Approach to Compute Terminal Velocity Vectors
To avoid computing Ω � i , the velocities can be obtained entirely in terms of OS equinoctial elements instead.First, knowing ŵ1 from Eq. ( 14), compute p 11 , p 21 from Biria [22]  as Next, from Biria and Russell [23], compute the following: f1 , ĝ1 , ̇f 1 , ̇ĝ 1 via Eqs.24, 127, and 128, cos L 1 , sin L 1 from Eqs. 90-91, and Li from Sect. 9, where the equation numbers refer to that reference.For convenience, and with i = 1, 2 , these equations are repeated here in order as and Then, compute ηi from Eq. (48) or (50) and i and ρi from where the expression for i is new and that for ρi is derived from Biria [28].Note that, alternatively, η2 does not need to be computed from Eq. ( 48 does not need to be computed from Eq. (69).By using angle sum identities to determine Q cos 2 and e sin f 2 in Eq. (48) and Eq. ( 69), respectively, η2 can be determined directly from η2 = ψ2 Q cos  2 and the computation of 2 can be bypassed.Next, to avoid inverse trigonometric functions, obtain cos L 2 , sin L 2 from standard angle sum identities after computing ΔL from Biria [22] as Then, propagate p 11 , p 21 via Eq. 63in Biria in [22], repeated here for convenience as to obtain p 12 , p 22 , which determine f2 , ĝ2 , ̇f 2 , ̇ĝ 2 as before in Eqs.(65-66).Finally, compute v i by applying to each terminal point Eqs.133 and 134 in Biria and Russell [23], repeated here as ) ̇zi = ρi  i +  i ηi .

Fig. 2
Fig.2Comparison of Lambert solution spaces under Vinti and Keplerian dynamics from zero-rev to 2-rev transfers with fixed r 1 , r 2 .Minimum-energy ellipses are marked in blue squares(Vinti)  and red x's (Kepler)

Fig. 3 2 38
Fig. 3 Direct comparison of Lambert solution spaces underVinti and Keplerian dynamics for 100-rev short-way transfers at Earth with fixed r 1 , r 2 .The minimum predicted TOF for existence of solutions differs by ≈ 8 hours; a 6.5-day transfer exists under Keplerian dynamics, but not under the more realistic Vinti dynamics.The minimum predicted energy differs by ≈ 1 km 2 /sec 2

Fig. 4 3 38
Fig.4 Inclination sweep over Vinti and Keplerian Lambert solution spaces for direct 20-rev short-way transfers at Earth with r 1 , r 2 fixed for each scenario

Table 1
Initial osculating Keplerian orbital elements used to generate r 1 , r 2 for various examples 2 are generated by propagating the set of osculating Keplerian element ICs in Table 1 through the time of flight given in Table 2, or roughly 14,156 sec, under Vinti dynamics, resulting in the approximate terminal positions r 1 = [5372.789,4437.582,668.070] ⊤ km and r 2 = [−6903.967,−892.533,1480.227]⊤ km.All examples in the present work use Earth parameters = 3.986004415 × 10 5 km 3 /sec 2 , R e = 6378.137km, and J 2 = 1.082636022984 × 10 −3

Table 2
Initial TOF with the associated number of revs used to generate r 1 , r 2 for the examples in Table1