Revisiting Universal Variables for Robust, Analytical Orbit Propagation Under the Vinti Potential

To meet the growing complexity and demands of modern spacecraft missions, analytical solutions to initial value problems see continued use, typically supporting global searches of large trajectory design spaces. These efforts often employ universal two-body orbit propagators for their recognized speed and robustness, but many applications, like active space debris removal, would benefit from a comparable propagator with greater accuracy. Vinti propagators, which consider planetary oblateness, may serve this purpose, but existing Vinti solutions possess computational difficulties in certain orbital regimes. To mitigate these deficiencies, the present study develops and validates an analytical, third-order universal Vinti propagator free of computational difficulties by leveraging standard, oblate spheroidal (OS) universal variables and OS equinoctial orbital elements. Accuracy of the third-order approximation is assessed for multiple examples across an array of orbital regimes. Computational runtime is also evaluated, and performance is directly compared to the benchmark Vinti6 algorithm. On average, the Vinti propagator implemented in this work is only slower than a typical universal Kepler propagator by a factor of 4.0 and slower than Vinti6 by a factor of 1.8, but with greater robustness than the benchmark. The new form of the equations of motion also has favorable implications for the associated boundary value problem.


Introduction
The Vinti gravitational potential [1,2], like that of the two-body problem [3,4], can describe or capture the predominant natural motion of a small object traveling near a large gravitating body.The latter exactly describes satellite motion around a spherically symmetric body or point mass and the former exactly describes motion around an oblate spheroid, which better approximates the Earth's basic shape.The utility of these gravitational models derives largely from the fact that they admit analytical solutions [4][5][6] to their respective initial value problems (IVPs), furnishing a past or future state given specified initial conditions.The Kepler or two-body prediction problem has been studied for several centuries and is one of the most fundamental problems in astrodynamics [4].As a natural analogue of the Kepler problem for oblate spheroids and one of the few integrable problems in celestial mechanics, the Vinti problem warrants further study [7].These gravitational models have merit because, for many of the bodies in the solar system, such as Earth and other planets, these computationally efficient approximations of the actual physics are good enough to support numerous applications, either directly or as building blocks of a larger architecture [8][9][10][11][12].But the rising complexity and innovation in modern spacecraft missions, like active space debris removal (ASDR), often call for more accuracy than two-body gravity in preliminary trajectory designs.Replacing this model with Vinti's potential, which reduces the modeling error by three orders of magnitude at Earth [1], is an attractive way to better satisfy these accuracy requirements and provides the incentive to study Vinti's potential in the current work.
The importance of more accurate models like the Vinti potential motivates the continued refinement of analytical Vinti IVP solutions to improve performance metrics like computational speed and robustness.These metrics are important for complex preliminary mission design applications that require global searches of very large trade spaces, such as those encountered in ASDR, a time-dependent traveling salesman problem with moving nodes [13].While numerical integration schemes may be simpler than analytical solutions, and more accurate, the integrator may be called on the order of a billion times in an ASDR application, and the implied long runtimes can make a numerical approach unappealing or less viable.In contrast, a hallmark of implemented analytical solutions to both Vinti and Kepler IVPs, accordingly called Vinti and Kepler propagators, is the innately high speed at which they predict or compute a spacecraft's ballistic trajectory.State-of-the-art Kepler propagators have a runtime of roughly a microsecond and existing Vinti propagators are only five times slower on average [2], noting that those tested in that study also included J 3 through Vinti's asymmet- ric potential [14,15].A prior study by Biria and Russell [12] found for a lowaltitude ASDR application that a 400 times speed-up was attainable with a Vinti propagator, using the "Vinti6" algorithm [2], as compared to an eighth-order Runge-Kutta numerical integrator taking 50 steps per spacecraft revolution over a 10-day transit.Similar speed-ups are expected from the solution proposed in this paper, although algorithm robustness is prioritized over computational speed.Specifically, the novel universal Vinti solution proposed in the current work aims to retain the described high computational speeds while leveraging astrodynamics concepts known to promote robustness.
Having motivated the utility of Vinti propagators, it is worth noting that, like Kepler propagators, a variety of source codes are publicly available.As fundamental tools, many analytical propagator codes are found in textbooks.For Kepler propagators, undergraduate textbooks like Vallado [8] and Curtis [16] offer straightforward implementations.For Vinti propagators, Vinti's graduate textbook edited by Der and Bonavito [2] contains six source codes in multiple programming languages.More recently, Biria and Russell [11,12,17] developed enhanced versions that remove computational difficulties and indeterminate forms, substantially improving accuracy for bounded orbits, especially near the equatorial regime.Their source code is also publicly available online. 1niversal solutions to these IVPs, such as those archived by Der and Bonavito [2], are attractive because they should in theory possess no computational difficulties or limitations in validity, being valid for all orbital regimes.As such, while Wiesel [18] and Wright [10] contributed notable advancements to Vinti theory using actionangle variables and numerical methods, the current effort is instead focused on universal techniques.Universal Kepler propagators are ubiquitous, computer programs commonly adopting the implementation of Bate et al. [4] in standard universal variables (UVs), which are notable for their robustness and connection to the Lagrange coefficients.Prior work on universal Vinti propagators [2,19] use what Herrick [20] called "unified variables", which, while retaining accuracy near the parabolic orbital regime, still utilize classical orbital elements and are not connected to oblate spheroidal (OS) Lagrange coefficients.Lacking a standard-UV solution to Vinti's IVP, existing Vinti propagators must rely more on carefully chosen numerical algorithms to create and increase robustness.In contrast, a Vinti propagator based on standard UVs arguably inherits more robustness naturally.While such a formulation may still be paired with the same robust numerical methods, the incorporation of additional robustness by design implies an extra layer of robustness that enhances the entire propagator.
The main contribution of this paper is the development of an analytical thirdorder solution to the Vinti IVP using an analog of standard UVs adapted to OS geometry.Getchell [19] developed a universal third-order solution with OS "unified variables" (source code is available from Der and Bonavito [2]), but his solution does not remove all of the computational difficulties promised by the introduction of UVs.Note that UVs are not expected to mitigate computational issues for nearly rectilinear Vinti orbits, corresponding to the so-called forbidden zone [21] inside which an analytical solution is not known.Apart from nearly rectilinear orbits, the solution proposed in the current study is valid for all orbital regimes and is free of all computational difficulties, accurately handling circular, equatorial, polar, nearly parabolic, and hyperbolic orbits.To accomplish this result, the approach taken identifies desirable features unique to different analytical solution methods and blends those components together, unifying all of the individual benefits into a single algorithm.Drawbacks of the individual algorithms are circumvented and only the benefits are retained.Specifically, the proposed solution uses the definition of the node from Vinti's 1969 solution [22], techniques to evaluate the integrals and avoid indeterminate forms at zero energy from Getchell's 1970 solution [19], and finally, to avoid the computation of singular orbital elements, the solution adopts the standard definition of OS UVs [23] and leverages the OS equinoctial orbital elements (EOEs) developed recently by Biria and Russell [17].Results are validated against numerically integrated Vinti trajectories to assess the accuracy of the third-order solution, where the third-order approximation refers only to the evaluation of various integrals, not the factoring of the quartics, which is carried out to double-precision accuracy.The presented algorithm's runtime is also evaluated and compared to earlier findings, and its overall performance is directly compared to the benchmark Vinti6 propagator.
Another important aspect of the current study is the implication that it has for the closely related, congruent boundary value problem (BVP) governed by the Vinti potential.The reformulation of the equations of motion in terms of OS classical orbital element differences and coupling to OS equinoctial elements not only supports the solution of the IVP, but also enables the definition of the BVP, which is explored in related work by the author [24,25].In the BVP application, Biria essentially uses the results of the present work to generalize Lambert's equation to a system of equations that can be solved iteratively, avoiding targeting algorithms and the need for shooting methods that are commonly employed in perturbed Lambert solvers.

Kinematic Equations
Vinti's symmetric potential [1,6] is written in terms of the oblate spheroidal (OS) coordinates , , as where is the gravitational parameter of the central body, is the OS coordinate equal to the semiminor axis of the instantaneous oblate spheroid, is the right ascension, and is the OS coordinate tied to latitude, approximately the sine of declination in an Earth application.The parameter c is fit to the oblateness term, J 2 , of the spherical harmonic expansion as c 2 = R 2 e J 2 , where R e is the equatorial radius.Formal definitions of these quantities are available in many references [1,6].The analytical solution to Vinti's IVP can be stated as x = f(t, t 0 , x 0 ) , where 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, t is the time, and the subscript " 0 " denotes initial conditions (ICs) or the initial value of a quantity. (1) Using Getchell's definitions [19] of the R j and N j integrals for j = 1, 2, 3 , Vinti's classical element solution [6] to the IVP in terms of right ascension is given by where R j are defined as N j are defined as and j and j are the constants of integration with = − 1 , = 2 , and Ω = 3 .These quantities are all defined in other references, but for convenience, note that 1 is the energy integral, 3 is the polar component of the angular momentum, 2 is similar to the total angular momentum and obtained from separating the Hamilton-Jacobi equation, is the time of OS periapsis passage, is the argument of OS periapsis, and Ω is the right ascension of the OS ascending node (spheroidal RAAN) [12].Further note that F( ) and G( ) are the quartics that must be factored to obtain the classical OS orbital elements a, e, and I = arcsin Q , which are the OS semimajor axis, eccentricity, and inclination, respectively.For universal variables, a is not directly pursued; instead, when factoring F( ) , two other quantities are obtained: p, the OS semilatus rectum, and , the inverse of a, where = −1∕a for elliptical orbits so that  < 0 for elliptical orbits ( 0 ≤ e < 1 ) and  > 0 for hyperbolic orbits ( e > 1 ).Orbital elements with a "0" subscript denote the prime constants [2]  obtained directly from the j Jacobi constants.For the lower limits of integration, p is the OS periapsis radius and Getchell's lower limit on the N j integrals is zero when J 3 = 0 in the potential.In the following, elements should generally be interpreted as spheroidal or OS unless noted otherwise.
(2) Going forward, it is helpful to continue to combine the contributions of Vinti [22] and Getchell [19] because their solutions possess important but different features that enable the current work.To exploit them, the solutions must be merged in some way.Specifically, many of Getchell's results are desirable because he removes indeterminate forms in the R j integrals.It is convenient to use his expressions for N j as well.Vinti's solution in terms of a slowly-varying OS RAAN [22], Ω � , is desirable because it ena- bles the coupling of OS equinoctial elements to universal variables, which seems necessary for a solution devoid of angle ambiguities and singular orbital elements.Getchell defined a quantity N 4 without explanation, but it is actually the slow, nonsensitive part of N 3 , which means it can be used to write the kinematic equations in terms of Ω � as where, if N 4 g is Getchell's definition of N 4 , then N 4 ≡ 3 N 4 g in Eq. ( 9).This redefini- tion is made so that N 3 and N 4 have consistent definitions.Using Getchell's notation, it is straightforward to show that where 1 and 2 are quantities defined by Getchell [19] and not related to , the true argument of OS latitude.Equation (10) relates to a different OS RAAN ( Ω � ≠ Ω ) and has not previously appeared in the literature.
To express the solution in terms of universal variables [23], a different form of the kinematic equations is required that depends on the difference in anomalistic and other angles.To obtain such a form, it is assumed that Eq. ( 9) is computed at two different times, t and t 0 , and the corresponding equations are then differenced.Since j are con- stants of the motion, the result is where Δ is shorthand for taking the difference of scalars ( Δt = t − t 0 ) or long expres- sions like R j and N j .The R j integrals are ultimately a function of anomalistic angles like true and eccentric anomaly, f and E, respectively, and the N j integrals a function of = f + � , where ′ is a different OS argument of periapsis ( ′ ≠ ).When differenced, ΔR 1 ≡ ΔR 1 (x, Δf ) , ΔR j ≡ ΔR j (Δf ) for j = 2, 3 , and ΔN j ≡ ΔN j (Δ ) for j = 1, 2, 3, 4 , where x is the OS universal variable identified by Biria [23].Notice that the time of flight (TOF), Δt , appears explicitly in Eq. (11).In this new form, Vinti's IVP is now stated as x = f(Δt, x 0 ).

Universal Spheroidal Time of Flight Equation
The first kinematic equation in Eq. ( 11) is the OS time of flight equation.After a considerable amount of algebra, and adopting Izsak's and Vinti's conic parameterization [5,6], this equation can be written to O(J 3  2 ) as ( 9) with Getchell's notation, where Biria [23] defines 0 and the UVs x and ẑ for OS geometry (under Vinti's potential), and also explains how the Stumpff functions, C(ẑ) and S(ẑ) , may be adapted to OS UVs.Equation (12), which directly relates TOF to the standard OS UVs for the first time, represents a significant development that is central to all of the present work, for both the IVP and BVP.While Eq. 36 in Getchell [19] requires the computation of the singular orbital element 1 before the root-solve of various orbital elements, Eq. ( 12) of this paper circumvents the computation of 1 and instead prescribes the root-solve of the orbital element differences discussed in the context of Eq. (11).Note in Eq. ( 12) that x = √ aΔE for elliptical orbits and ẑ = − x2 .The symbol W k represents a special recursive function defined by Getchell [19], seeded with , and C 6 = 3Q 2 1 ∕8 .Then, only even T k matter, and the even Q k T k can be expressed as by manipulating Getchell's equations.Substituting for Q sin = , Eq. ( 13) can be simplified as and Eq. ( 14) can be differenced at two times to obtain Δ(Q k T k ) as Computational details for ΔW k and ΔT k are discussed in Appendix A.

Universal Spheroidal Argument of Periapsis or Apsidal Drift Equation
The second kinematic equation in Eq. ( 11) is the OS argument of periapsis equation, which can be stated as or, using Getchell's notation [19], to O(J 3  2 ) as ( 12) but the drift in the OS argument of periapsis, Δ � , is not apparent or accessible in this form.Note in Eq. ( 17) that a number of simplifications have been made assuming J 3 = 0 .Now, to extract Δ � from Eq. ( 17), there are multiple ways to proceed.The simplest option is to substitute Δ = Δf + Δ � into the isolated Δ term and rearrange Eq. ( 17) so that Δ � is on the left-hand side as or, defining the right-hand side (RHS) in Eq. (18) as g2 , identify a recursive relation- ship as For convenience, define Δ Ñ2 as the modified ΔN 2 integral in Eq. ( 18) with Δf in the secular term.When the constants and Δf are known, Eq. ( 18) can be used to iteratively solve for Δ � .If the derivative of g2 Δ � with respect to Δ � has an absolute value less than unity, then the method of successive approximations will converge [26].To verify this sufficient condition, the derivative is evaluated as Each derivative in Eq. ( 20) contains terms either equal to factors of Q 2 or 2 , or the product of multiple trigonometric functions, which all have a maximum absolute value of unity (see Appendix C).Each of these terms has a coefficient of . Therefore, if evaluated from left to right, the terms have maximum values on the order of 10 −3 , 10 −6 , and 10 −9 for the Earth, respectively.Since |g � 2 Δ � | < 1 , the method of successive approximations will converge, and because the derivative is small, with |g � 2 Δ � | < O J 2 , it will converge rapidly.The maximum error of the nth approximation can be stated as , which means that for an Earth application with a worst-case initial guess error of 1 radian, the fourth iteration will have an error of roughly 10 −12 in the approximation of the root Δ � .In practice, the initial error |Δ � 0 − Δ � | is orders of magnitude smaller than unity, even for a many-revolution (many-rev) orbit, e.g. for a representative low-Earth orbit (LEO) 20-rev scenario inclined 30 degrees, the observed initial error in apse line drift is only 2.6 × 10 −3 radians or about 0.15 degrees.( 17) Note that alternatives to Eq. ( 18) may be obtained.For example, define the secular coefficient for the R 2 integral as C R 2 and subtract C R 2 Δ � from both sides of Eq. ( 17) to obtain where ΔW k p represents the periodic part of each ΔW k term.Solving Eq. ( 21) for Δ � leads to For convenience, define Δ R2 as the modified ΔR 2 integral in Eq. (22) with Δ in the secular term.

Universal Spheroidal RAAN or Nodal Drift Equation
The third kinematic equation in Eq. ( 11) is the spheroidal RAAN equation, which in general enables the computation of the spheroidal RAAN drift as To O(J 3  2 ) , Eq. ( 23) can be expressed with Getchell's notation as where , and d 6 = 5Q 3 1 ∕16 , which means that C 1k and C 2k are nonzero for all k.Observe that using Vinti's potential with J 3 = 0 has resulted in a number of simplifications and will also cause the C 1k terms to cancel with the C 2k terms for odd k.

Universal Vinti Orbit Propagator
The new form of the kinematic equations established in the previous section enables the development of a novel, universal Vinti orbit propagator, but several additional components are required to furnish a complete algorithm.It is helpful to have in Fig. 1 an outline of the author's computational procedure for the propagator, showing how to initialize the algorithm and how to use UVs and OS equinoctial elements in concert.The computational procedure is discussed in detail in the following sections.

On the Spheroidal Equinoctial Elements
It seems a nonsingular solution based on the standard UVs, as presented here, would not be possible without the definition of OS equinoctial orbital elements [17], because a nonsingular coordinate transformation between the ECI frame and an orbital frame is required before and after the root-solve.A new, simple method based on vectors is proposed for the computation of the OS equinoctial elements p 1 and p 2 , which physically represent the components of the OS ascending node vector [17].The OS position vector in ECI coordinates is given by Fig. 1 Proposed and validated computational procedure for a universal Vinti propagator and taking the time derivative gives the OS inertial velocity as where ρ and ̇ρ denote the OS position unit vector and its inertial velocity, respectively, the subscripts x, y denote x and y components, and the superscript N denotes the Newtonian or ECI frame.Letting the superscript R denote the rotating frame attached to the OS orbital plane with angular velocity = [0 0 Ω� ] ⊤ , the OS velocity relative to the R frame attached to the OS orbital plane is given by where R ρ is constrained to lie in the orbital plane by definition.The OS specific angular momentum vector relative to the OS orbital plane can then be computed as and the OS equinoctial unit vector ŵ parallel to R h as ŵ = R h∕‖ R h‖ in ECI coordi- nates.Note that while ‖ R h‖ can be computed directly from Eq. ( 30), consideration of its physical meaning as the total angular momentum normal to the OS orbital plane points to a simple expression for the magnitude of R h as which can be verified from Eq. (30).Finally, p 1 and p 2 can be computed as using the same formulas as Danielson et al. [27] with K as the retrograde factor.The guidelines for computing Ω� offered by Biria and Russell [17] may be used in Eq. ( 29), although a better alternative is proposed here that is free of singularities at the poles.
One of the two exact expressions for Ω� put forward by Biria and Russell [17] has the following form: where Q 1 is O(J 2 ) and determined from factoring the G( ) quartic.The 3 ∕(1 − 2 ) expression is indeterminate near the poles, but the term it multiplies on the right in Eq. ( 33) can be manipulated.First, rewrite the right term inside the brackets as Then, rationalize the numerator and simplify to obtain so that the problematic 1 − 2 term can be canceled out of the denominator.Substi- tuting the result in Eq. ( 35) back into Eq.( 33) yields a final expression for Ω� , analo- gous to those for ̇f or ψ , that is exact, always nonsingular, and independent of any method for evaluating the integrals, given by When obtaining Ω� from OS orbital elements, an alternative expression for 3 may be used to give where c 2 = R 2 e J 2 and Q 1 = −c 2 0 ∕(p 0 S 1 ) have been used to write Ω� in a more familiar form.Equation (37) indicates that Ω� is O(J 2 ) and proportional to cos I as expected.All quantities on the right-hand side are constant except for and , which cause the rate to vary over time.Interestingly, since 0 = 0 for parabolic orbits, the term on the right inside the brackets does not contribute to Ω� in that orbital regime: The various connections to standard concepts for two-body dynamics evident in Eqs.(30) and (32) highlight the role and utility of OS vectors in a way not previously shown.While the approach by Biria and Russell [17] establishes useful relationships for computing p 1 and p 2 in certain contexts, the physical insight and sim- plicity offered by the above algorithm is considered an improvement over previous methods.Furthermore, the introduction of an exact, nonsingular expression for Ω� in Eq. ( 36) finally removes all singularities from the OS equinoctial coordinate transformation, greatly simplifying the transformation and making it easier to use.

Well-Behaved Exact Forms for the Mean Frequencies
Let C R j be Getchell's three secular coefficients for the respective R j integrals and C N j be his three secular coefficients for the N j integrals (see Appendix B).From here, defining a mean-motion-like quantity as for convenience, the anomalistic and draconitic mean frequencies can be written respectively as and From Eqs. ( 40) and (41), the secular rate for OS argument of periapsis can be obtained as Similarly, assuming 3 is factored out of C R 3 and C N 3 , the secular rate for OS RAAN can be determined as where 3 = √ p 0 1 − Q 1 S 0 cos I has been used to make the dependence on cos I explicit.

Initial Guesses for a One-Dimensional, Two-Tier, Nested Root-Solve
Taking an approach analogous to Bate et al. [4], good initial guesses can be derived depending on the orbital regime.By choosing the definition of Vinti UVs proposed by Biria [23], the arguments presented by Bate et al. to derive various initial guesses for x also follow for Vinti dynamics.New arguments are presented to derive initial guesses for Δ � , which is zero under Keplerian dynamics.
For elliptical orbits [23], x = √ aΔE , which means x = 2 √ a after one orbital period in the sense of one revolution of an anomalistic angle.Following Bate et al. [4], note that where t p is the anomalistic orbital period.Let t p = 2 ∕2 1 = 1∕ 1 and solve for x to obtain Substituting 2 1 from Eq. (40) into Eq. ( 45 and noting that a good initial guess for x in the elliptical regime can finally be expressed as Using mean frequencies or Eq.(42), a good guess for Δ � in the elliptical regime (  < 0 ) is given by For hyperbolic orbits with a large change in hyperbolic eccentric anomaly, the arguments by Bate et al. [4] can be applied to obtain the approximations: where "e" in this section refers to the exponential function.Recall the universal TOF equation in Eq. ( 12) and neglect not only the  0 x term, as Bate et al. did, but also O(J 2 ) and smaller terms, to give Substituting Eq. (49) into Eq.(50) gives and solving for √ −ẑ leads to after taking the logarithm.Noting from Bate et al. [4] that x is positive when Δt is positive, the initial guess for x when dealing with hyperbolic orbits can finally be stated as Over the span of a typical hyperbolic orbit, ′ will not change much, so the drift Δ � can be guessed as For nearly parabolic orbits, the same assumptions are made for Vinti dynamics as for Keplerian dynamics.Accordingly, an appropriate guess for x in this orbital regime is xguess = 0 , which may be applied for nearly parabolic orbits in the bounded or unbounded case, perhaps with a tolerance of || < 10 −5 , depending on how many terms are retained for the approximations of the C and S functions in this regime [4].When nearly parabolic orbits are escape trajectories ( ≥ 0 ), then Δ � guess = 0 is a good guess, as for hyperbolic orbits.
Striking similarities are observed for the initial guesses of x between Keplerian and Vinti dynamics.While the guess for x is different between elliptical, nearly par- abolic, and hyperbolic regimes, the guess for Δ � is binary: zero for escape trajec- tories ( ≥ 0 ), but with a dependence on mean frequencies for bounded trajectories (  < 0 ).When employed in a root-solve procedure, these initial guesses for x and Δ � will speed up convergence and increase robustness of the root-solve in a univer- sal Vinti orbit propagator. (49 29 Page 16 of 38

Initialization and Root-Solve Procedure
An outline of the initialization and root-solve procedure employed in this study and depicted in Fig. 1 is offered here.For the majority of the algorithm's initialization phase, the procedure developed in Biria and Russell [17] is used to convert an initial ECI state to initial OS equinoctial elements and compute all the intermediate constants.
In other words, the first eight steps in the flow chart follow the procedure in that reference, except that Ω� is computed from Eq. ( 36) and p 1 0 , p 2 0 from Eqs. (26-32), the improved algorithm developed in an earlier section.Next, since the current propagator uses universal variables, Getchell's method [19] is used to compute the coefficients of the six integrals and certain initial values, with subtle differences to avoid indeterminate forms.For the ΔR j integrals, the coefficients A k are computed from Getchell's equa- tions [19] and the initial values e k V k 0 are computed from Eqs. (72-74) in Appendix A. Additional constant coefficients generated by the recursive functions W k in Eq. ( 69) and e k V k are computed and stored for later use at this stage, noting that e j is factored out of each coefficient such that the exponent matches the number of trigonometric functions multiplied together in a given term.The e j factors are discarded from the coefficients and accounted for by grouping them with the appropriate terms in V k .This bookkeep- ing process avoids divisions by e that could cause issues for nearly circular orbits when e ≈ 0 .For the ΔN j integrals, the coefficients C k , C 1k , C 2k and initial values Q k T k 0 are required.Expressions for coefficients are given explicitly just before Eq. ( 13) to compute C k for ΔN 1 , explicitly in Eq. ( 17) for ΔN 2 , and in Eq. ( 25) to compute C 1k , C 2k for ΔN 4 .The Q k T k 0 terms can be computed from Eqs. ( 14) and (68).After computing all of the necessary constants, the next step is to initialize the rootsolve.Of the three kinematic equations, the root-solve only involves Eq. ( 12) for x and Eq.(18) for Δ � .Equation (24) for ΔΩ � is decoupled and is not used until after the root- solve.The root-solve is initialized by guessing x and Δ � using the methods developed in the previous section.Guess x from Eq. (47) in the elliptical case or Eq.(53) in the hyperbolic case, or set xguess = 0 in the parabolic case.Guess Δ � from Eq. (48) in the elliptical case or Eq.(54) in the hyperbolic or parabolic case.
A desired root-solve procedure follows the calculation of the two initial guesses described above, where the one-dimensional root-solve depicted in Fig. 1 searches for the value of the universal variable x that solves Eq. ( 12).It is one-dimensional in the sense of Vinti's solution [6] or that of Getchell [19], where the equation governing the OS argument of periapsis, ′ , is used in an intermediate step.In those references, with J 3 = 0 , the equation for ′ is used to sequentially obtain the value of .In the current approach, where the value of Δ � is unknown but a good initial guess is available, a straightforward method of successive approximations is applied to Eq. (18) to solve for Δ � at each iteration of the root-solve on x .This approach is taken while acknowledg- ing that better alternatives may exist.For the described outer loop, a Newton-Raphson algorithm combined with a bisection method [28] is used, a technique that Press et al. [28] describe as a fail-safe alternative to the less robust Newton-Raphson method.While robustness may improve with a variable order Laguerre's method [9,29], its implementation in this framework is left to future work.For iteration j = 0 , set xj = xguess and Δ � j = Δ � guess , and compute ẑj , C j , and S j from the definitions described in Biria [23], Δf j from Eq. (84), and Δ j from As a safeguard, a quadrant check is performed to ensure sgn Δf j = sgn(Δt) , which could be violated near the full-rev boundary if the derivatives are inaccurate; Δf j is appropriately shifted to the range 0, sgn(Δt)2 as necessary.The quantities ΔR 1 j and ΔN 1 j are required to calculate the right-hand side of Eq. (12).Compute ΔW k j from Eq. ( 70) and ΔR 1 j from Also compute ΔR 2 j in the outer loop as To obtain ΔN 1 j , a converged value of Δ � j is required, which is furnished by the method of successive approximations.First, for iteration i = 0 , set Δ � ji = Δ � j and Δ ji = Δ j , and compute Δ(Q k T k ) ji from Eqs. ( 15) and (65-68).Then, compute Δ � j(i+1) from Eq. ( 18).If |Δ � j(i+1) − Δ � ji | > 10 −12 or a desired tolerance, apply Eq. (55) to update Δ ji to Δ j(i+1) , set i = i + 1 , and return to Eq. ( 15).Otherwise, Δ � j(i+1) has converged to within the desired tolerance, which was found to take only a few iterations in practice since the initial guess is a good estimate.Use Eqs. ( 55) and (15) to update Δ ji and Δ(Q k T k ) ji one more time, respectively, and then compute ΔN 1 j as and the right-hand side of Eq. ( 12) as Δt j = ΔR 1 j + c 2 ΔN 1 j .The Newton update Δx j can be obtained from where and its derivative (55) Δ j = Δf j + Δ � j .
(56)  86), (92), and (93), with computational details on the first derivatives given in Appendix C, to be calculated in the outer loop.Using a simple stopping criterion like the one in Curtis [16], if the update |Δx j | > 10 −12 or a desired tolerance, set j = j + 1 , update xj , ẑj , C j , S j , Δf j , and return to Eq. (55) to update Δ j and repeat the process, including the evaluation of Eq. ( 18) as part of the method of successive approximations.Otherwise, xj has converged to within the desired toler- ance and the root-solve is terminated.Then, compute ΔΩ � from Eqs. (24)(25).Note that any desired one-dimensional root-solve on x could be substituted into the above process, but it must be ensured that any time x is updated, such as right before ter- minating the root-solve, all derived quantities that flow from x need to be updated as well to be consistent, including Δf , Δ , ΔW k , Δ(Q k T k ) , and the evaluation of Eq. ( 18) as part of the method of successive approximations.
At this stage in the propagator, differential OS orbital elements have been obtained for the classical OS elements that vary with time: Δf and UVs for the anomalistic motion, and Δ � and ΔΩ � to describe the total drift of the trajectory within the rotat- ing OS orbital plane and the total drift of the OS orbital plane itself, respectively.The intermediate determination of the secular drift and evolution of the system [6] has been circumvented, and singular elements have been bypassed.If indeterminate calculations are to be avoided, a nonsingular orbital reference frame is required to map the state back to the ECI frame, and the OS equinoctial reference frame can be employed again for this purpose.

Propagating OS Equinoctial Elements
At the conclusion of the root-solve and after solving for ΔΩ � , three of the relevant equi- noctial elements can be propagated, which amounts to adding ΔL to the initial true lon- gitude L 0 as and appropriately rotating the OS ascending node vector by ΔΩ � as Then, with the final Ω� obtained from Eq. ( 36) and ρ and η from Eqs. 21 and 17 in Vinti [22], respectively, as where e sin f and Q sin are computed from angle sum identities and ̇f and ψ from Biria and Russell [17], the transformation to ECI coordinates is almost complete.Finally, the steps outlined in Biria and Russell [17] can be followed to compute the (61) remaining final scalar time derivative, L , final OS equinoctial unit vectors, f and ĝ , their time derivatives, ̇f and ̇ĝ , and ultimately the ECI position and velocity.

Mitigating Degradation of Accuracy in Elliptical Many-Rev Scenarios
In elliptical many-rev scenarios, Δt , x , ẑ , and root-solved angular quantities can grow very large if a naive approach is adopted that does not leverage knowledge of the orbital period, incurring computational errors that grow over time and become large for ultra-high-rev scenarios.To avoid this degradation of accuracy, a robust version of the universal Vinti orbit propagator should reduce the TOF and angular quantities to a range within a single orbital period.Such a reduction is complicated to achieve under the Vinti potential, but it is possible and has been implemented in this work.
Under Keplerian dynamics, it is trivial to apply the modulo operation to the orbital period to achieve the desired reduction [4].When this approach is adapted to the Vinti propagator developed in this study, however, the orbital "period" is not constant; it is slightly different for each rev, where the (anomalistic) period is defined here as the time it takes for |Δf | to go from 0 to 2 for a particular rev.The desired quantity then becomes the time Δt N required to complete N revs, which can be found iteratively.The iteration required is nowhere near as complex as the root-solve depicted in Fig. 1, mainly because Δ � N is the only unknown, with Δf N = sgn(Δt)2 N , xN = sgn(Δt) √ a2N , ẑN = (2N) 2 , C(ẑ) N = 0 , and S(ẑ) N = 1∕ẑ N .The propagator uses the method of successive approximations to iteratively solve Eq. (18) for Δ � N with the initial guess Δ � N = sgn(Δt) ω� s t p N guess , and then, once converged, it is trivial to compute Δt N from Eq. ( 12).Computational details, including a method for guessing N, are provided in Appendix D. The implied small increase in compute time is considered a small sacrifice in exchange for accuracy preservation.

Examples: Universal Propagator Accuracy, Speed, and Robustness
Having introduced, in full, a novel Vinti orbit propagator, the performance of the presented algorithm must be assessed in terms of accuracy, runtime, and robustness.First, accuracy is evaluated for a small set of carefully chosen orbital regimes of interest.Then, computational speed is measured from a large set of random initial conditions, from which robustness statistics are also extracted.Finally, additional performance comparisons are made to the benchmark Vinti6 algorithm [2].

Accuracy
The accuracy of Vinti dynamics has been evaluated against spherical harmonics techniques in multiple prior studies [11,30,31] and is not examined further in this 29 Page 20 of 38 work.The modeling of non-conservative forces under Vinti's potential has been explored as well, Vinti [32] and Tong and Wu [33] taking a generalized approach while Sherrill [34] and Watson et al. [35] focused on atmospheric drag.While perturbations due to non-conservative forces are important for certain applications, they are outside the scope of this work.Instead, the present work focuses on evaluating the accuracy of the third-order approximation of the R j and N j integrals in differ- ent orbital regimes and the computational runtime of the proposed universal Vinti propagator.
The propagator's accuracy in four different orbital regimes is examined in Fig. 2 to validate the approach and implementation, with ICs for each scenario given in Table 1.Other parameters are = 3.986004415 × 10 5 km 3 /sec 2 , R e = 6378.137km, and J 2 = 1.082636022984 × 10 −3 .Depending on the orbital regime, the trajectories are initialized with Keplerian mean anomaly, M k , or true anomaly, f k .The examples include a geostationary orbit (GEO) (Fig. 2a) and a Molniya orbit (Fig. 2b) taken directly from Appendix C in Der and Bonavito [2], in addition to two examples of inclined escape trajectories: an exactly parabolic case (Fig. 2c) and a retrograde hyperbolic case (Fig. 2d).Figure 2 shows the 3D trajectory and log-scale position error for each example.The position error metric is the magnitude of the vector error between the third-order analytical solution and the numerically integrated equations of motion in ECI coordinates [36].Numerical integration is performed in MATLAB using ode45, a variable-step fourth-order Runge-Kutta method, with a 2.224 × 10 −14 relative accuracy tolerance and a 1.0 × 10 −20 absolute tolerance.In the Vinti propagator, the tolerance on Δ � is set to 10 −15 rad in both the primary and full-rev root-solve algorithms, to extract as much accuracy as possible, and x is iterated upon until the solution cannot be further improved in double precision.
The accuracy and reliability of the propagator is found to be high across all orbital regimes, with errors in Fig. 2 seen to be within about 10 −8 km over a range of flight times for escape trajectories and within roughly 10 −6 km for long many-rev propaga- tions up to 20 days.OS Lagrange coefficients are used to help assess and monitor convergence [23].The secular error growth observed may be attributed to two factors: 1) the computation of approximate orbital element coordinates (the momenta are computed to double-precision accuracy), meaning the representation of the Vinti trajectory in orbital element space does not exactly match its representation in initial ECI coordinates; 2) the secular terms in the propagation step are only accurate to the third order in J 2 .Additional testing showed that the present algorithm generally matches the accuracy of the Vinti6 algorithm [2] when the many-rev accuracy retention measure is applied.These details are presented in Section 4.3 along with some runtime comparisons.
In accordance with the goals of this work, the long propagation times described above are chosen to assess the algorithm's performance and identify any limitations, not to demonstrate its performance in practice.For example, while a Vinti GEO propagation of a few days may be useful in preliminary mission design, time spans in a catalog building application may be limited to a day or so and is subject to nonphysical constraints like data processing schedules.That being said, the practicality of long Vinti propagations in GEO depends not only on the application, but also on how the algorithm is used.As noted by Getchell [19], for comparisons to perturbation theories or numerical integration, fitting the Vinti solution to an ephemeris would extract the greatest performance from the algorithm, where principal differences would be due to neglecting equatorial obliquity.
The results above also shed light on some important qualitative differences that can occur between various dynamical models.In particular, Fig. 2c is a notable example of a set of ICs that leads to a qualitative discrepancy between the Kepler and Vinti propagators.Two-body dynamics predict a bounded orbit with e k = 0.99998 (see Table 1 for the exact value), while Vinti dynamics pre- dict an escape trajectory with e = 1.0, which is exactly parabolic.Other tests identified cases where the inverse occurs, where the Keplerian dynamics predict an escape trajectory and the Vinti dynamics predict a bounded one.For objects orbiting oblate bodies, the Vinti theory offers a quick way to evaluate boundedness more accurately, of practical importance for preliminary, interplanetary mission design that typically starts with two-body dynamics.When considering orbit insertion at Saturn, for example, such oblateness effects would be even more pronounced, and employing a universal Vinti propagator would imply a more efficient orbit design process.

Speed and Robustness
Computational speed and robustness are 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 universal Vinti propagator is benchmarked in runtime against a universal Kepler propagator.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.ICs and times of flight are given in Table 2: 1,536,000 ICs spanning a range of orbital regimes are propagated for Δt values ranging from 1 sec to 1 day, specifically for [1; 10; 100; 1,000; 10,000; 86,400] sec, totaling 9,216,000 function calls.For each of r p k , e k , and i k , 40 sam- ples are randomly chosen from a uniform distribution, except that the exact values r p k = 2,500 km, e k = (0, 1), and i k = (0, 90, 180) deg are explicitly sampled.The value of Ω k is explicitly set to 60 deg for all ICs because RAAN intuitively should not affect convergence properties or runtime; the values of k and f k are similarly not random and are explicitly given in Table 2, overall ensuring the inclusion of potential stress cases associated with the following orbit types: nearly circular and/ or exactly equatorial (direct and retrograde), exactly polar, initialized exactly on a pole, and nearly parabolic (both elliptical and hyperbolic).
For the described times of flight and ICs, the Vinti algorithm described above only takes ≈ 4.15 μs on average to compute a set of final ECI position and velocity vectors.The tolerance on Δ � for this data set is reduced to 10 −13 rad in the main root-solve because the additional accuracy is not considered worth the trade in runtime.The Vinti propagator produces zero errors and is found to be only ≈ 3.98 times slower than a comparable universal Kepler propagator based on the Bate et al. [4] formulation, which takes an average of ≈ 1.04 μs on the same machine and also pro- duces zero errors.As before, both algorithms are set up to iterate until the solution cannot be further improved in double precision.Runtimes observed here are consistent with Der and Bonavito's findings [2], where the Vinti propagators they examined, which add J 3 to the baseline Vinti dynamics via Vinti's asymmetric potential [14,15], run five times slower than a Kepler propagator on average.It should be noted that other known universal Vinti propagators [2] require a function call to a universal Kepler propagator to obtain an initial guess for the independent UV in the Vinti algorithm.In contrast, the initial guesses implemented in the current algorithm are simple closed-form expressions that can be computed quickly.Assuming an average 1.04-μs Kepler call, the current Vinti propagator is saving about (1 − 4.15 μs∕5.19 μs) × 100 = 20% in runtime relative to an identical imple- mentation using a Kepler initial guess, which highlights the value of the proposed initial guess technique.Note that the baseline Kepler propagators are not identical between the current study and those performed by Der and Bonavito [2].As such, conclusions on the relative computational speeds of the current algorithm and propagators like Vinti6 cannot be drawn with the given information.The direct runtime comparisons required to assess the efficiency of the Vinti6 algorithm are explored in the following section.

Comparisons to the Vinti6 Algorithm
For completeness, the presented Vinti propagator is directly compared to the Fortran version of Vinti6 [2] in terms of accuracy, robustness, and computational efficiency.Accuracy and robustness are compared in Fig. 3 using the same four examples defined in Table 1.While Vinti6, which is based on Vinti's asymmetric potential [14], is capable of modeling both J 2 and J 3 , the J 3 effects are omit- ted (set to zero) to increase the fairness of the accuracy and robustness comparisons, the number of mathematical computations in the code remaining unchanged.Note that, before performing the comparisons, the author corrected an error in the original Vinti6 Fortran code located inside the root-solve, where the code originally exits the loop upon convergence without updating and subsequently the T k functions.The algorithm error is eliminated by removing line 614 and inserting it immediately after line 628 in the code. 2 Evidently, this correction is only needed in the Fortran version of Vinti6 and not in the C version.Now, to obtain accuracy assessments for Vinti6, the error metric discussed in the previous section is simply applied to the Vinti6 outputs and that error is plotted side-by-side with the error in Fig. 2. Der and Bonavito [2] note that Kepler1, the default Kepler propagator in Vinti6, may be replaced with any UV Kepler propagator, and so it is replaced with the Kepler propagator described in the previous section to make the comparisons as fair as possible.This choice to replace the Kepler algorithm promotes a fair comparison from the perspective of both robustness, since its robustness is demonstrated in the previous section, and runtime, since the measured runtime is recorded in the previous section.The panels in Fig. 3 are arranged in the same order as the panels in Fig. 2 to make comparisons easier.Focusing first on Fig. 3b, d, which are the Molniya and hyperbolic cases, respectively, Vinti6 is observed to agree almost exactly with the presented algorithm.This level of agreement is expected because the approximations in both algorithms are correct to O(J 3  2 ) and both leverage some form of universal variables.However, in Fig. 3a, c, which are the GEO and parabolic cases, respectively, there are some clear discrepancies.Over the chosen time spans, while the presented Vinti propagator maintains millimeter-level or better accuracy, the Vinti6 errors are consistently on the order of a few kilometers for these two examples, representing two orbital regimes that should be easily handled by universal techniques, but which are not handled well by Vinti6.Further investigation suggested that the large Vinti6 errors result from the algorithm converging to poor solutions for the GEO and parabolic examples.Based on the GEO results and a few additional quick tests, the accuracy degradation is found to exist for nearly circular and/or equatorial orbits and is worse for the exactly equatorial scenario.Further study is required to understand the cause of the parabolic discrepancy.The findings not only indicate that the presented Vinti propagator is a significant improvement over a well-known benchmark, but they also suggest that the coupling of UVs to equinoctial orbital elements is a necessary theoretical development for enhancing the robustness of analytical Vinti orbit propagators.Given the unequal levels of robustness observed between the presented Vinti propagator and Vinti6, it may be considered unfair to perform a runtime comparison of the two algorithms.The lack of robustness in one algorithm raises a number of challenges for runtime comparisons and raises questions on the usefulness of making the assessment.Generally, such a measurement of the computational speed of Vinti6 could be misleading.One area of concern is that when Vinti6 converges to poor solutions, the algorithm may converge more quickly or slowly than average, depending on the situation, which could unfairly benefit or hurt its performance.Further study would be required to determine if the convergence to poor solutions is predictable.Nevertheless, it is considered worthwhile to still evaluate the average Vinti6 runtime for archival purposes.For the times of flight and ICs specified in Table 2, Vinti6 takes an average of ≈ 2.30 μs to compute final ECI state vectors on the same machine described in Section 4.2.While Vinti6 can potentially compute a solution in 55% of the time as the presented Vinti propagator on average, the presented orbit propagator is much more robust.The practitioner must decide whether trading robustness for speed is acceptable for their application, keeping in mind that the speedup offered by Vinti6 may be less in practice because the runtime measured in this study includes Vinti6 solutions that failed to converge correctly.The speed and robustness of Vinti6 also has a strong dependence on the speed and robustness of the underlying UV Kepler propagator, which must be taken into consideration when deciding which algorithm to use.

Conclusions
A universal analytical solution is developed for Vinti's symmetric, unperturbed dynamical problem, which uses oblate spheroidal (OS) geometry to add planetary oblateness to the dynamics.Standard OS universal variables are meshed with OS equinoctial orbital elements to eliminate all singularities and computational difficulties from the solution, rendering it universally valid for bounded and unbounded trajectories at any inclination, including polar and direct or retrograde equatorial orbits.The forbidden zone associated with nearly rectilinear orbits, a rarely encountered region in practice, remains the only regime for which an analytical solution is not known.Central to this approach, a generalized universal Kepler equation is developed, and a two-tier, nested root-solve is proposed to solve it, where steps are taken to guarantee that the simpler inner root-solve converges in a few iterations.Interestingly, the entire root-solve is viewed as operating in a rotating OS orbital reference frame, and the equinoctial elements enable a nonsingular mapping between this frame and the inertial frame.The algorithm design permits forward or backward propagation.
The presented Vinti solution is obtained with notable simplifications and improvements in accuracy relative to earlier solutions.In particular, the coordinate transformations between inertial position and velocity and OS equinoctial elements are now exact, where they previously relied on approximations to evaluate the transformations near the poles.Exactness is enabled by the derivation of an exact, nonsingular representation of the rate at which the OS ascending node drifts over time.The computation of OS ascending node vector components from inertial coordinates is also simplified using a new approach based on vectors that streamlines the algorithm.Otherwise, the accuracy of the solution is as high as previous analytical solutions, where the quartics are factored to double-precision accuracy and the integrals are evaluated to an accuracy of O(J 3  2 ) .A technique to preserve solu- tion accuracy in many-rev scenarios is also offered.Ultimately, with these mechanisms in place, the new Vinti propagator is demonstrated to be more robust than the benchmark Vinti6 with J 3 = 0 , especially in the nearly circular, equatorial, and parabolic regimes.When applied at Earth, the implemented Vinti algorithm, which includes novel, simple initial guess techniques that do not require a Kepler solution, is only slower than a universal Kepler propagator by a factor of 4.0 on average.The new Vinti algorithm is slower than Vinti6 by an average factor of 1.8, but the increased runtime comes with substantial gains in robustness for popular orbital regimes.
The formulation of a solution to Vinti's dynamical problem that combines OS universal variables and equinoctial elements implies a level of robustness and efficiency, both expected and observed, that can be passed on to related applications.Depending on the primary body's shape, preliminary mission design that typically relies on two-body approximations can benefit from the increased accuracy, enjoying a minor sacrifice in compute time relative to two-body propagation and a significant speed boost relative to numerical integration that includes J 2 effects.The solution also prescribes singularity-free partial derivatives and an analytical, universal state transition matrix that is nonsingular for bounded and unbounded orbits in any nondegenerate orbital regime.Finally, the new form of the kinematic equations lends itself to the definition of a boundary value problem as a system of equations, which will be explored in a follow-on study.

Appendix A: Computational Details for Kinematic Equations
Recall Eq. ( 15) as and observe that For the ΔW k terms, note first from Getchell [19] that if W 0 = f and W 1 = (f + eV 1 )∕p , then Equation (69) can be differenced at two times to obtain ΔW k as where powers of eccentricity are grouped with the V k terms in the same way that powers of Q are grouped with the T k terms to avoid angle ambiguities.Specifically, recalling that Getchell [19] writes the recursive function V k , seeded with V 0 = f and V 1 = sin f , as the e k V k terms for k > 0 can be expressed as where To better see how computations of ΔW k may be carried out, it is convenient to split ΔW k into secular and periodic parts.Due to the recursive generation of terms defined in Eqs. ( 70) and (72), the secular terms contain Δf with various constant coefficients, while the periodic terms have the form with different constant coefficients.
(68)   (80) where the required derivatives of ΔR 2 and ΔN 2 , respectively, can be obtained explic- itly as and The simplest way to evaluate the RHS derivatives in Eq. ( 89) is to note that taking the derivative of ΔW k with respect to Δf recovers the original integrand that gener- ates the recursive function.The derivatives can be stated in terms of the integrand, found in Eq. 24 of Getchell [19], as Equation (91) is the simplest expression for computing the derivative, and it can be further verified if desired by evaluating Eqs.(104) and (105).The RHS partial derivatives in Eq. (90) can be evaluated with Eqs.(96-97).Finally, multiplying Eq. (88) by Eq. (86) yields the desired result as With Eq. (86) and Eq.(92) established, the first derivative of the OS time of flight equation, Eq. ( 12), can be pursued with respect to the independent variable, x .The first derivative can be stated as or it can be simplified to  (93) where has been substituted into the first term in Eq. ( 12) using the relation from Biria [23], Eq. (91) has been substituted, and other derivatives can be computed as and Noting dΔT 0 ∕dΔ = 1 , Eqs. ( 96) and (97) can be written explicitly for even k up to k = 6 as and As an alternative to verifying Eq. ( 91), take the derivative of Eq. ( 70) with respect to Δf to give (94) Taking the first derivative of each secular term will leave only the various constant secular coefficients found from Eqs. ( 70) and (72).For the periodic terms, taking the derivative of Eq. ( 75) with respect to Δf yields which after multiplying by the associated constant periodic coefficients from Eqs. ( 70) and (72) results in the desired derivatives of the periodic terms.Sum the secular and periodic terms to obtain the total derivative.While the above derivatives are technically correct, as the exact derivatives of the approximate equations being root-solved, these derivatives can be approximated very well at Earth by the exact derivatives of the actual R j and N j integrals to yield extensive simplifications.For the N j integrals, the simplifications result from the direct application of elliptic integrals.The N 1 and N 2 integrals can be expressed in terms of the incomplete elliptic integrals of the first and second kind as [2] where are the incomplete elliptic integrals of the first and second kind, respectively, and the constant modulus q in Vinti's notation is related to the parameter k 1 in Getchell's notation [19] as Although the derivative dΔN j ∕dΔ is desired for j = 1, 2 , it is equivalent to dN j ∕d because d ∕dΔ = 1 and the integrals evaluated at t 0 are not a function of Δ .Substituting Eq. (107) into Eq. ( 106 and taking the derivative of the resulting equations for N j with respect to yields simple nonsingu- lar expressions as (106) F( , q) = ∫ 0 1 − q 2 sin 2 − 1 2 d ; E( , q) = ∫ 0 1 − q 2 sin 2 after applying Leibniz rule and some algebraic manipulation to eliminate indeterminate forms in dN 1 ∕d for nearly parabolic orbits.The R j integrals are treated simi- larly, expressed in different forms as where the integrals on the first line in Eqs. ( 109) and ( 110) are taken from Getchell [19] and those on the second line are obtained by substituting Getchell's differential relationships for , into the original integrals as needed, noting for Getchell's X variable that d X = dx .
The R 1 integral is written with respect to x because dΔt∕dx is ultimately desired, while that of R 2 is written with respect to f because dΔR 2 ∕dΔf is required in Eq. (92).Note that dΔR 1 ∕dx is equivalent to dR 1 ∕dx and dΔR 2 ∕dΔf to dR 2 ∕df for reasons analogous to those cited for the N j integrals and because df ∕dΔf = 1 .Tak- ing the derivatives of Eqs. ( 109) and (110) with respect to x and f , respectively, yields simple nonsingular expressions as and Finally, the derivative of Eq. ( 12) can be computed as by substituting the appropriate expressions from Eqs. (108), ( 112), ( 113), (86), and (92), the latter simplifying to (109) (110)

29
Page 18 of 38 is obtained from Eqs. (

Fig. 2
Fig. 2 Propagated trajectories at Earth (left column) and log-scale position errors (right column) for four orbital regimes.The new analytical Vinti propagator agrees well with numerically integrated Vinti trajectories

Fig. 3
Fig.3 Comparisons between the presented Vinti propagator and Vinti6 using log-scale position errors for four orbital regimes.The new analytical Vinti propagator consistently performs well, while Vinti6 fails for the GEO and parabolic cases
f ; dΔ e k V k p dΔf = e k cos k f − (k − 1)e k sin 2 f cos k−2 f , k ≥ 2,

Table 1
Initial osculating Keplerian orbital elements for Fig. 2 examples a The GEO example is taken from Appendix C, Scenario IV in Der and Bonavito [2]; full numerical values used for quantities truncated above are a k = 42,164.169613508km b The Molniya example is taken from Appendix C, Scenario III in Der and Bonavito [2]; full numerical values used for quantities truncated above are a k = 26,628.1361947432km, e k = 0.741696641081651, M k = 144.0088647361997deg c The Parabolic and Hyperbolic examples are initialized with the osculating Keplerian periapsis radius, r p k , instead of the semimajor axis.For the parabola, the full numerical value of the Keplerian eccentricity is e k = 0.99997885396048585

Table 2
Initial osculatingKeplerian orbital elements for the assessment of computational runtime, using each of the following six Δt values:

a
The "min" function and f k∞ = arccos(−1∕e k ) are only invoked fore k ≥ 1 .If f k = f k∞ , then f k is reset to 5 − f k∞ in degrees