Paolo Izzo
email: izzogianpaolo@gmail.com

Lamberto Dell'elce

Pini Gurfil

Aaron J Rosengren

Daniel J Scheeres

Nonsingular vectorial reformulation of the short-period corrections in Kozai's oblateness solution

Keywords: Kepler problem, Oblateness, Nonsingular elements, Averaging, Artificial satellite theory

We derive a new analytical solution for the first-order, short-periodic perturbations due to planetary oblateness and systematically compare our results to the classical Brouwer-Lyddane transformation. Our approach is based on the Milankovitch vectorial elements and is free of all the mathematical singularities. Being a non-canonical set, our derivation follows the scheme used by Kozai in his oblateness solution. We adopt the mean longitude as the fast variable and present a compact power-series solution in eccentricity for its short-periodic perturbations that relies on Hansen's coefficients. We also use a numerical averaging algorithm based on the fast-Fourier transform to further validate our new mean-to-osculating and inverse transformations. This technique constitutes a new approach for deriving short-periodic corrections and exhibits performance that are comparable to other existing and well-established theories, with the advantage that it can be potentially extended to modeling non-conservative orbit perturbations.

Introduction

In the brief span of time after the launch of Sputnik, a whole succession of analyses was devoted to the problem poised by the drag-free motion of an artificial satellite about an oblate planet, employing almost every known perturbation method. Although, in a sense, the problem is a classic one that also occurred among the natural satellites, in the applications of artificial satellite motion it was necessary to obtain a more general, detailed, and accurate solution. The most intricate and notable investigations were presented by [START_REF] Brouwer | Solution of the problem of artificial satellite theory without drag[END_REF], [START_REF] Garfinkel | The orbit of a satellite of an oblate planet[END_REF], and [START_REF] Kozai | The motion of a close Earth satellite[END_REF] in the celebrated 1959 issue of the Astronomical Journal. These authors treat the first-and second-order secular perturbations, as well as the firstorder short-periodic (related to the satellite's mean motion) and long-periodic (related to the evolution of the argument of the perigee) perturbations of the orbital elements, where order here refers to the oblateness parameter.

Astronomical experience amply bears out the notion of separation of perturbing effects into periodic and secular variations and the distinction between fast and slow time variables concerning the motion of satellites [START_REF] Lara | A new radial, natural, higher-order intermediary of the main problem four decades after the elimination of the parallax[END_REF]. The device of canonical transformations employed by [START_REF] Brouwer | Solution of the problem of artificial satellite theory without drag[END_REF] and [START_REF] Garfinkel | The orbit of a satellite of an oblate planet[END_REF] permits a deeper understanding of the difference between periodic and secular perturbations and provides a systematic procedure for the inclusion of higher-order terms (Kozai 1962b). On the other hand, the technique used by [START_REF] Kozai | The motion of a close Earth satellite[END_REF], known formally as the method of analytic continuation [START_REF] Scheeres | Orbital Motion in Strongly Perturbed Environments: Applications to Asteroid[END_REF], while applicable to any set of orbital elements, canonical or otherwise, can be cumbersome beyond the first order. The Poincaré-von Zeipel method of canonical transformations, though still in use [START_REF] Nie | Semi-analytical model for third-body perturbations including the inclination and eccentricity of the perturbing body[END_REF], was significantly generalized by [START_REF] Hori | Theory of general perturbations with unspecified canonical variables[END_REF] and [START_REF] Deprit | Canonical transformations depending on a small parameter[END_REF] based on Lie series; the latter becoming the sine qua non of modern perturbation theory in celestial mechanics [START_REF] Deprit | The main problem of artificial satellite theory for small and moderate eccentricities[END_REF][START_REF] Kaufman | First order semianalytic satellite theory with recovery of the short period terms due to third body and zonal perturbations[END_REF]).

Kozai's gravity solution to the artificial satellite problem became the basis for the simplified general perturbations theory, SGP, that would be supplanted by SGP4, an analytical solution for satellite short-term prediction using the two-line element sets that instead has its roots in Brouwer's gravitational theory [START_REF] Hoots | History of analytical orbit modeling in the U. S. Space Surveillance System[END_REF]). [START_REF] Brouwer | Solution of the problem of artificial satellite theory without drag[END_REF] developed his original solution in Delaunay variables, the canonically-conjugate counterpart of the Keplerian orbital elements used by [START_REF] Kozai | The motion of a close Earth satellite[END_REF], which, like the classical ones themselves, are singular for circular and equatorial orbits. The mathematical singularities associated with zero eccentricity and vanishing line of nodes that plague these solutions can be removed by reformulating them in nonsingular variables, such as those of Poincaré [START_REF] Lyddane | Small eccentricities or inclinations in the Brouwer theory of the artificial satellite[END_REF][START_REF] Breakwell | On error bounds and initialization in satellite orbit theories[END_REF], Hill [START_REF] Izsak | A note on perturbation theory[END_REF][START_REF] Aksnes | On the use of the Hill variables in artificial satellite theory: Brouwer's theory[END_REF], or the equinoctial set [START_REF] Gim | State transition matrix of relative motion for the perturbed noncircular reference orbit[END_REF][START_REF] Alfriend | Spacecraft Formation Flying: Dynamics. Control and Navigation[END_REF][START_REF] Le Fèvre | Compliance of disposal orbits with the French Space Operations Act: the good practices and the STELA tool[END_REF]),1 or Euler-parameter-based elements [START_REF] Alfriend | Spacecraft Formation Flying: Dynamics. Control and Navigation[END_REF]. A simplified form of the Lyddane-modified Brouwer theory, optimized in SGP4 for the rapid propagation of satellite ephemerides of the space object catalog (SOC) [START_REF] Hoots | Reformulation of the Brouwer geopotential theory for improved computational efficiency[END_REF], is the basis for many tracking and prediction operations. Moreover, SGP4, or its deep space equivalent, SDP4, must be used to convert the mean, orbit-averaged, TLEs of the SOC into an osculating set of elements for use in special perturbation theory to obtain more accurate predictions [START_REF] Levit | Improved orbit prediction using two-line elements[END_REF].

The principal observable features due to Earth's oblateness are a secular precession of the orbit plane about the polar axis and a steady motion of the major axis in the moving orbit plane. Well-known applications of this particular problem are the two inclined and elliptical orbit systems (12-hr Molniya and 24-hr Tundra), placed at the reputed critical inclination of ∼63.4 • so that the apsidal precession freezes on average. The critical inclination is an intrinsic singularity in the artificial satellite theories of our forebears [START_REF] Coffey | The critical inclination in artificial satellite theory[END_REF], and hitherto remains a fundamental issue in all modern reformulations of the main problem [START_REF] Breakwell | On error bounds and initialization in satellite orbit theories[END_REF][START_REF] Aksnes | On the use of the Hill variables in artificial satellite theory: Brouwer's theory[END_REF]Lara 2015a, b). The orbit-averaging technique fundamentally involves removing all terms that depend on the fast-varying mean anomaly, thus retaining the (secular and long-period) mean-element motion. As noted by [START_REF] Kozai | The motion of a close Earth satellite[END_REF], the long-periodic perturbations of the first order come from the terms of the second order, and the singularity associated with the critical inclination only arises when such long-period effects are retained. Nevertheless, the distinction between secular and mean elements has often been muddled in the literature, and tabulated mean-to-osculating transformations [START_REF] Gim | State transition matrix of relative motion for the perturbed noncircular reference orbit[END_REF][START_REF] Schaub | Analytical Mechanics of Space Systems[END_REF]) apply Brouwer's full periodic corrections. This approach will inevitably lead to significant errors near the critical inclination, but can be properly amended by neglecting the long-periodic terms [START_REF] Breakwell | On error bounds and initialization in satellite orbit theories[END_REF].

An important application of the mean-to-osculating transformation is in the computation of the nominal osculating orbits from the frozen-orbit conditions determined in mean-elements space [START_REF] Gurfil | Motion near frozen orbits as a means for mitigating satellite relative drift[END_REF]. Frozen orbits correspond to equilibria for the averaged equations of motion, and, in the oblateness model, occur when secular effects due to even zonal harmonics are canceled by the long-period perturbations of the odd harmonics. While modern formulations compute frozen orbits directly from the non-averaged equations, based on the underlying quasi-periodic structure of librations around these mean equilibria, explicit analytical solutions form the starting point for the numerical optimization process. Recovering the short-periodic effects needed for this initialization can be troublesome at the critical inclination, whose precise frozen-orbit location is very sensitive to model truncation [START_REF] Lara | Exploring sensitivity of orbital dynamics with respect to model truncation: The frozen orbits approach[END_REF].

Here, we present a new formulation of the mean-to-osculating and inverse conversions for first-order oblateness perturbations based on the Milankovitch elements (Rosengren andScheeres 2013, 2014). We use the direct method of [START_REF] Kozai | The motion of a close Earth satellite[END_REF], as further elucidated by [START_REF] Scheeres | Orbital Motion in Strongly Perturbed Environments: Applications to Asteroid[END_REF], and present an explicit analytical short-period correction in vector form that is valid for all elliptical orbits. We adopt the mean longitude as the fast variable and present a compact power-series solution in eccentricity for its short-periodic perturbations that can be truncated to achieve the necessary accuracy. We establish a numerical averaging approach based on the fast-Fourier transform [START_REF] Uphoff | Numerical averaging in orbit prediction[END_REF][START_REF] Ely | Transforming mean and osculating elements using numerical methods[END_REF] and make detailed comparisons between our vectorial solution and the classical Brouwer-Lyddane (BL) theory. For the latter, we adapt the more streamlined formulas presented in [START_REF] Schaub | Analytical Mechanics of Space Systems[END_REF] and [START_REF] Gim | State transition matrix of relative motion for the perturbed noncircular reference orbit[END_REF].

Problem formulation

Analytical averaging

The basic idea in orbit-averaging methods is to obtain approximate equations for the system evolution that contain only slowly changing variables by exploiting the presence of a small dimensionless parameter that characterizes the size of the perturbation. The tacit assumption is that the perturbing forces are sufficiently weak so that these approximate mean equations of motion can be used to describe the secular and long-period orbital evolution. The perturbation equations in celestial mechanics, relating the time variation of the orbit parameters to the perturbing accelerations, in Gauss or Lagrange form, are nonlinear, nonautonomous, firstorder differential equations: [START_REF] Alfriend | Spacecraft Formation Flying: Dynamics. Control and Navigation[END_REF][START_REF] Scheeres | Orbital Motion in Strongly Perturbed Environments: Applications to Asteroid[END_REF]

ẋ = g(x, t), ( 1
)
123 in which g(x, t) is assumed to be T -periodic in time t. Equation ( 1) is trivially solved when = 0, yielding the integrals (Keplerian elements) in the unperturbed problem. The method of averaging consists of replacing Eq. ( 1) by the averaged autonomous system [START_REF] Sanders | Averaging Methods in Nonlinear Dynamical Systems[END_REF])

ẋ = g(x), g(x) = 1 T T 0 g(x, t) dt, (2)
where the average is performed over time, and it is understood that x in the integrand is to be regarded as a constant during the averaging process. The basis for this approximation is the averaging principle, which states that in the general, non-resonant case, the short-period terms removed by averaging cause only small oscillations that are superimposed on the long-term solution described by the averaged system.

Comparison between numerical integrations and the mean solution will in general show a divergence between the two as a result of an inconsistent choice of initial conditions, as depicted in Fig. 1. This offset can be understood by decomposing the osculating elements into mean and short-period components [START_REF] Kozai | The motion of a close Earth satellite[END_REF][START_REF] Scheeres | Orbital Motion in Strongly Perturbed Environments: Applications to Asteroid[END_REF]:

x(t) = x(t) + x sp (t).
(3) Differentiating, we can obtain an approximate equation governing the short-period dynamics:

ẋsp (t) = ẋ -ẋ ∼ = g(x, t) -g(x), (4)
in which the mean state x has replaced the corresponding osculating elements x in the dynamical equations and where the small parameter has been omitted without lack of generality. From Eqs. (3) and ( 4), the short-periodic perturbations can be obtained as [START_REF] Kozai | The motion of a close Earth satellite[END_REF])

x sp (t) = dx sp (t) -dx sp = [g(x, t) -g(x)] dt - 1 T T 0 [g(x, t) -g(x)] dt 2 .
(5)

Accordingly, it can be seen that x sp = 0, which is also the implicit assumption in the averaging process. The interpretation of this result is that given an initial condition for a state x 0 = x(t 0 ), the mean equations have to be initialized at a value provided by Eqs. ( 3) and ( 5) as

x 0 = x 0 -x sp (t 0 ), ( 6
)
to have the averaged dynamics track the true evolution more closely. The removal of time, or analogously of the mean anomaly, requires computing the quadrature of functions depending implicitly on this variable through the true anomaly. The time averaging is performed over a periodic motion having a period much shorter than the time that characterizes the evolution of the dynamical system; this periodicity necessarily implies that averaging is taken for elliptical orbits. Given a quantity g(x, M) representing the right-hand side of the equations of motion, defined as a function of the dimensionless time variable M (the mean anomaly) in addition to the other orbital elements given by x, the average, Eq. ( 2), can be redefined as where the orbital elements x are held constant in the integration. Although the average is defined with respect to mean anomaly, it is often more easily calculated by means of the true or eccentric anomaly, f and E, respectively, using the differential relationships

g(x) = 1 2π 2π 0 g(x, M) dM, (7)
dM = r a dE = r 2 ab d f , (8)
in which a and b = a √ 1e 2 are the semi-major and semi-minor axes, respectively, and e is the eccentricity, yielding the equivalent forms for averaging:

g(x) = 1 2π 2π 0 g(x, M) dM = 1 2πa 2π 0 g(x, E) r dE = 1 2πab 2π 0 g(x, f ) r 2 d f . (9)
Note that r can be expressed in terms of f and E as

r = ⎧ ⎨ ⎩ a(1 -e 2 ) 1 + e cos f = H 2 /μ 1 + e cos f , a(1 -e cos E). ( 10
)
Here, H is the specific angular momentum and μ is the gravitational parameter.

The Milankovitch elements consist of the two fundamental vectorial integrals of motion, namely, the eccentricity vector e and angular momentum vector H. These vectors can be parameterized in terms of the Keplerian elements relative to an inertial frame:

H = H ĥ = H (sin i sin Ω x -sin i cos Ω ŷ + cos i ẑ), (11)
e = eê = e (cos ω cos Ω -cos i sin ω sin Ω) x +(cos ω sin Ω + cos i sin ω cos Ω) ŷ + sin i sin ωẑ , ( 12
)
where i is the inclination, Ω is the right ascension of the ascending node, and ω is the argument of periapsis. Because of the orthogonality constraint (i.e., H • e = 0), we need a sixth scalar element to fully define an orbit [START_REF] Roy | Studies in the application of recurrence relations to special perturbation methods III. Non-singular differential equations for special perturbations[END_REF]. Adopting the mean longitude l = ω + Ω + M, the non-averaged equations of motion for an arbitrary disturbing 123 acceleration a d can be stated in 'dyadic' form2 as [START_REF] Battin | An Introduction to the Mathematics and Methods of Astrodynamics[END_REF][START_REF] Rosengren | On the Milankovitch orbital elements for perturbed Keplerian motion[END_REF]:

ė = 1 μ v • r -H • a d = g e (x, t), ( 13a
) Ḣ = r • a d = g H (x, t), ( 13b
) l = n + - e μ(1 + √ 1 -e 2 ) H (ê • r)r + (r + p)(ê • v) θ - 2 na 2 r + r • ẑ H (H + H • ẑ) H • a d , = n + g l (x, t), ( 13c
)
where n 2 = μ/a3 , p = H 2 /μ, and θ = ĥ • r. Note that Eq. (13c) consists of terms that collect the contributions due to the three components in the radial, transverse, and normal directions of the disturbing acceleration and is valid for elliptical orbits in which e < 1.

The position and velocity vectors, r and v, may be expressed as

r = r cos f ê + sin f ê⊥ , ( 14
) v = μ H -sin f ê + (e + cos f ) ê⊥ , ( 15
)
where ê = e/e, ê⊥ = ĥ • ê, and ĥ = H/H . The Gauss equations have time-varying terms multiplying the accelerations involving the true anomaly, and thus they must each be averaged separately. The mean evolution of these elements can be computed as

ė = 1 2π 2π 0 ė dM = g e (x), (16a)
Ḣ = 1 2π 2π 0 Ḣ dM = g H (x), (16b) l = 1 2π 2π 0 l dM = n + g l (x). ( 16c
)
The approximate short-period equations of motion for each element can then be formulated by subtracting Eq. ( 16) from Eq. ( 13), while holding all orbital elements but f constant. [START_REF] Kozai | The motion of a close Earth satellite[END_REF], in his solution employing the classical elements a, e, i, Ω, ω, and M, obtained the non-averaged, mean, and short-period equations of motion using the Lagrange planetary equations. We note that a more general procedure has been outlined herein.

Particularly, for e and H, following Eq. ( 5), we have

de sp (t) = g e (x, t) -g e (x) dt, ( 17a
) d H sp (t) = g H (x, t) -g H (x) dt, (17b)
so that e sp (t) = de sp (t)de sp , (18a)

H sp (t) = d H sp (t) -d H sp . ( 18b
)
Special care, however, must be taken in the case of the mean longitude due to the presence of the mean motion appearing without any factor in Eq. (13c) [START_REF] Kozai | The motion of a close Earth satellite[END_REF]. Expanding the osculating n into a first-order Taylor series about the mean elements, and following Eq. ( 4), we can write the approximate equation governing the short-period dynamics as

lsp (t) = l -l ∼ = g l (x, t) -g l (x) + ∇ e n(x) • e sp (t) + ∇ H n(x) • H sp (t), ( 19
)
where

∇ e n(x) = ∂n ∂ e x=x = - 3μ 2 H 3 1 -e 2 e, (20)
∇ H n(x) = ∂n ∂ H x=x = - 3μ 2 H 5 (1 -e 2 ) 3/2 H. ( 21
)
Accordingly,

dl sp (t) = g l (x, t) -g l (x) dt + ∇ e n(x) • e sp (t) dt + ∇ H n(x) • H sp (t) dt, (22)
and

dl sp = 1 T T 0 g l (x, f ) -g l (x) dt 2 + ∇ e n(x) • 1 T T 0 e sp (t) dt 2 + ∇ H n(x) • 1 T T 0 H sp (t) dt 2 . ( 23
)
As a result, the short-periodic perturbations of l are given by l sp (t) = dl sp (t)dl sp .

(24)

Numerical averaging

The conversion between mean and osculating elements can be obtained numerically, as previously shown by [START_REF] Walter | Conversion of osculating orbital elements into mean elements[END_REF], [START_REF] Uphoff | Numerical averaging in orbit prediction[END_REF][START_REF] Ely | Transforming mean and osculating elements using numerical methods[END_REF]. Here, we will exploit a numerical implementation of the near-identity transformation in Eq. (3) between averaged and osculating Milankovitch elements to validate our subsequent analytical developments. The mapping is obtained by numerically solving

∂ x sp ∂ t = [g (x, t) -g (x)] , T 0 x sp (x, t) dt = 0. (25)
The boundary condition in Eq. ( 25) guarantees that the oscillations of x sp (t) are unbiased with respect to x. The formal solution is then [START_REF] Sanders | Averaging Methods in Nonlinear Dynamical Systems[END_REF])

x sp (x, t) = - i n ∞ k=1 c k (x) k e iknt , ( 26
)
where c k (x) are the Fourier coefficients of g (x, t), and i is the imaginary unit.

123

In this study, we numerically approximate Eq. ( 26) by truncating the series and by using the fast-Fourier transform (FFT) algorithm as discussed in [START_REF] Ely | Transforming mean and osculating elements using numerical methods[END_REF], under the assumption that g is analytic on the continuation of nt in a non-vanishing complex strip. This assumption guarantees that the Fourier coefficients exhibit exponential decrease as a function of k. It results in a numerical averaging scheme that we call "FFT transformation" from now on. This scheme yields a first-order approximation of the motion, similarly to the Brouwer-Lyddane or Milankovitch ones. For the mean orbit longitude, we incorporate the additional terms in Eq. ( 19), resulting from the Taylor series of n, into Eq. ( 25), before performing the numerical quadrature.

Analytical short-period correction for oblateness perturbations

The quadrupolar (i.e., J 2 -truncated) disturbing function arising from an oblate planet can be stated in a general vector expression as [START_REF] Scheeres | Orbital Motion in Strongly Perturbed Environments: Applications to Asteroid[END_REF])

R = μJ 2 R 2 2r 3 1 -3(r • p) 2 , ( 27
)
where J 2 is the second zonal harmonic coefficient, R is the mean equatorial radius of the planet, and r = r/r from Eq. ( 14). Note that the Earth's spin axis p is assumed to be fixed in inertial space, and, as such, is aligned with ẑ.

The perturbing acceleration is then given by ∂R/∂ r as

a = - 3μJ 2 R 2 2r 5 1 -5(r • p) 2 r + 2(r • p) p . ( 28
)
Accordingly, following Eq. ( 13), the perturbation equations can be stated as

ė = 3J 2 R 2 2r 5 1 -5(r • p) 2 H • r -2(r • p)( v • r -H) • p , ( 29a
) Ḣ = - 3μJ 2 R 2 r 5 (r • p) r • p, (29b)
l = n + 3μJ 2 R 2 2r 5 1 μ(1 + √ 1 -e 2 ) H (e • r)(1 -3(r • p) 2 ) + 2(r + p)(e • v)(r • p)( θ • p) + 2r 2 na 2 (1 -3(r • p) 2 ) - 2(r • p) 2 (H • p) H (H + H • p) . ( 29c
)
The averaged equations of motion for the eccentricity and angular momentum vectors are given by [START_REF] Ward | On the secular variations of the elements of satellite orbits[END_REF] and [START_REF] Rosengren | Long-term dynamics of high area-to-mass ratio objects in high-Earth orbit[END_REF]. Averaging Eq. ( 29c) directly requires computing the quadrature of various dyadics and higher rank tensors of the dynamical variables, the details of which we omit. The mean equations can be stated as

ė = - 3n J 2 R 2 4 p 2 1 -5( ĥ • p) 2 ĥ + 2( ĥ • p) p • e, (30a)
Ḣ = 3Hn J 2 R 2 2 p 2 ( ĥ • p) ĥ • p, (30b)
l = n + 3n J 2 R 2 4 p 2 1 -e 2 3( ĥ • p) 2 -1 + 5( ĥ • p) 2 -2( ĥ • p) -1 , ( 30c
)
where the bar operator is omitted from the elements because there is no ambiguity in what follows, i.e., all variables are averaged variables. Note that Eq. (30c) can be seen as the sum of the classical secular precession rates arising from planetary oblateness, l = Ṁ + ω + Ω, where

Ṁ = n + 3n J 2 R 2 4 p 2 1 -e 2 3 cos 2 i -1 , ω = 3n J 2 R 2 4 p 2 5 cos 2 i -1 , Ω = - 3n J 2 R 2 2 p 2 cos i. ( 31
)
Following the procedure outlined in Sect. 2.1 and detailed in Appendix B, the shortperiodic perturbations in the eccentricity vector, angular momentum vector, and mean longitude can be stated as

e sp (t) = - 3J 2 R 2 2 p 2 I 2 + 2(I 1 + eI I 11 + I I I 111 -I I I 122 -5I V 122 )(ê • p)(ê ⊥ • p) -(2 I I I 112 + 5 I V 112 )(ê • p) 2 + (2 I 2 + 2e I I 12 + 2 I I I 112 -5 I V 222 )(ê ⊥ • p) 2 ê -I 1 -Me + 2( I 2 -e I I 12 -I I I 112 + I I I 222 -5 I V 112 )(ê • p)(ê ⊥ • p) + (3Me/2 + 2I 1 + 2I I I 122 -5I V 111 )(ê • p) 2 + (3Me/2 -2eI I 22 -2I I I 122 -5I V 122 )(ê ⊥ • p) 2 ê⊥ + Me(ê ⊥ • p)( ĥ • p) ĥ -2e I I 12 (ê • p) + I I 22 (ê ⊥ • p) p , ( 32a
)
H sp (t) = - 3H J 2 R 2 p 2 I I 12 (ê • p)( ĥ • p) + (I I 22 -M/2)(ê ⊥ • p)( ĥ • p) ê -(I I 11 -M/2)(ê • p)( ĥ • p) + I I 12 (ê ⊥ • p)( ĥ • p) ê⊥ + I I 11 -I I 22 + I I 12 (ê ⊥ • p) 2 -(ê • p) 2 ĥ , ( 32b
)
l sp (t) = 3J 2 R 2 2 p 2 e 1 + √ 1 -e 2 I 1 -3 I V 111 (ê • p) 2 + 2 I V 112 (ê • p)(ê ⊥ • p) + I V 122 (ê ⊥ • p) 2 + 2 I I I 222 -I I I 112 + I V 222 -I V 112 (ê • p)(ê ⊥ • p) + 2 (I I I 122 + I V 122 ) (ê • p) 2 -(ê ⊥ • p) 2 -2 3 1 -e 2 + ĥ • p 1 + ĥ • p I I 11 (ê • p) 2 + 2 I I 12 (ê • p)(ê ⊥ • p) + I I 22 (ê ⊥ • p) 2 + 2 1 -e 2 I 0 - 1 2 M 1 -e 2 3( ĥ • p) 2 -1 + 5( ĥ • p) 2 -2( ĥ • p) -1 + 3eH np 2 1 -e 2 2 IV + e II (ê • p)(ê ⊥ • p) + I 2 1 + 2(ê ⊥ • p) 2 -5 I V 112 (ê • p) 2 -5 I V 222 (ê ⊥ • p) 2 -2 I I I 112 + e I I 12 (ê • p) 2 -(ê ⊥ • p) 2 + 6H np 2 (1 -e 2 ) 3/2 II(ê • p)(ê ⊥ • p) -I I 12 (ê • p) 2 -(ê ⊥ • p) 2 , ( 32c
)
where Roman numerals I 1 , I 2 , I I 11 , . . . designate various functions of true anomaly, I 2 , I I 12 , . . . represent the difference of these trigonometric expressions from their averaged values, I 1 , . . . are indefinite integrals of the previous core functions, and I 2 , . . ., represent differences between the doubly-integrated expressions, their averaged values, and the previous core averages. While somewhat cumbersome, the adopted notation mirrors the derivation of Appendix B and is otherwise systematic and methodical. The needed results pertaining to the solution, Eq. ( 32), are given by

I 0 = f + e sin f , ( 33a
)
I 1 = 1 12 12e f + (12 + 9e 2 ) sin f + 6e sin 2 f + e 2 sin 3 f , ( 33b
)
123

I I 11 = 1 12 6 f + 9e sin f + 3 sin 2 f + e sin 3 f , ( 33c
)
I I 22 = 1 12 6 f + 3e sin f -3 sin 2 f -e sin 3 f , ( 33d
)
I I I 111 = 1 96 (36e f + 72 sin f + 24e sin 2 f + 8 sin 3 f + 3e sin 4 f ) , ( 33e
)
I I I 122 = 1 96 (12e f + 24 sin f -8 sin 3 f -3e sin 4 f ) , ( 33f
)
I V 111 = 1 240 180e f + 30(6 + 5e 2 ) sin f + 120e sin 2 f + 5(4 + 5e 2 ) sin 3 f +15e sin 4 f + 3e 2 sin 5 f , ( 33g
)
I V 122 = 1 240 60e f + 30(2 + e 2 ) sin f -5(4 + e 2 ) sin 3 f -15e sin 4 f -3e 2 sin 5 f , ( 33h
) I = - 1 12 (12 + 3e 2 )(cos f -X 0,1 0 ) + 6e(cos 2 f -X 0,2 0 ) + e 2 (cos 3 f -X 0,3 0 ) , ( 34a
)
I I = - 1 12 3e(cos f -X 0,1 0 ) + 3(cos 2 f -X 0,2 0 ) + e(cos 3 f -X 0,3 0 ) , ( 34b
)
I I I = - 1 96 24(cos f -X 0,1 0 ) + 12e(cos 2 f -X 0,2 0 ) + 8(cos 3 f -X 0,3 0 ) + 3e(cos 4 f -X 0,4 0 ) , ( 34c
)
I I I = - 1 96 72(cos f -X 0,1 0 ) + 12e(cos 2 f -X 0,2 0 ) -8(cos 3 f -X 0,3 0 ) -3e(cos 4 f -X 0,4 0 ) , ( 34d
)
I V = - 1 240 30(2 + e 2 )(cos f -X 0,1 0 ) + 60e(cos 2 f -X 0,2 0 ) + 5(4 + 3e 2 )(cos 3 f -X 0,3 0 ) + 15e(cos 4 f -X 0,4 0 ) + 3e 2 (cos 5 f -X 0,5 0 ) , ( 34e
)
I V = - 1 240 30(6 + e 2 )(cos f -X 0,1 0 ) + 60e(cos 2 f -X 0,2 0 ) -5(4 -e 2 )(cos 3 f -X 0,3 0 ) -15e(cos 4 f -X 0,4 0 ) -3e 2 (cos 5 f -X 0,5 0 ) , ( 34f
) I = - 1 12 ∞ k=1 1 k (12 + 3e 2 )C 0,1 k + 6eC 0,2 k + e 2 C 0,3 k sin k M, ( 35a
)
I I = - 1 12 ∞ k=1 1 k 3eC 0,1 k + 3C 0,2 k + eC 0,3 k sin k M, ( 35b
)
I I I = - 1 96 ∞ k=1 1 k 24C 0,1 k + 12eC 0,2 k + 8C 0,3 k + 3eC 0,4 k sin k M, ( 35c
)
I V = - 1 240 ∞ k=1 1 k 30(2 + e 2 )C 0,1 k + 60eC 0,2 k + 5(4 + 3e 2 )C 0,3 k + 15eC 0,4 k + 3e 2 C 0,5 k sin k M, ( 35d
)
I V = - 1 240 ∞ k=1 1 k 30(6 + e 2 )C 0,1 k + 60eC 0,2 k -5(4 -e 2 )C 0,3 k -15eC 0,4 k -3e 2 C 0,5 k sin k M, ( 35e
)
123

II = - 1 6 ∞ k=1 1 k 3eS 0,1 k + 3S 0,2 k + eS 0,3 k cos k M, (36a) IV = - 1 48 ∞ k=1 1 k 6(2 + e 2 )S 0,1 k + 36eS 0,2 k + (28 + 9e 2 )S 0,3 k + 18eS 0,4 k + 3e 2 S 0,5 k cos k M , ( 36b
)
where all intermediate terms are given in Appendix B. They are obtained by also using the auxiliary formulas given in Appendix A. Thus, following Eq. ( 6), given the initial osculating state (e 0 , H 0 , l 0 ), the mean equations of motion, Eqs. ( 30), have to be initialized as

⎧ ⎨ ⎩ e 0 = e 0 -e sp (t 0 ), H 0 = H 0 -H sp (t 0 ), l 0 = l 0 -l sp (t 0 ). ( 37
)

Validation of the Milankovitch scheme

In this section, we consider an extended grid of initial conditions and test the developed Milankovitch formulation against both BL and the fully numerical transformation of Sect. 2.2. For Brouwer-Lyddane, we have verified that the more streamlined formulas presented in [START_REF] Schaub | Analytical Mechanics of Space Systems[END_REF] have been correctly transcribed according to their original sources, excepting the missing sin(2ω) factor in the long-period terms of M, ω, and Ω, which has only been noted in the recent erratum of the latest edition of this widely used monograph.3 Furthermore, rather than using Lyddane's adhoc modification, [START_REF] Gim | State transition matrix of relative motion for the perturbed noncircular reference orbit[END_REF] developed a new theory based on Brouwer's generating function that uses equinoctial elements. While still invalid at the critical inclination, [START_REF] Gim | State transition matrix of relative motion for the perturbed noncircular reference orbit[END_REF] concluded from various numerical simulations that their method produces reasonable results within 0.25 • of this small divisor. Being inherently rooted in Brouwer's theory, it was expected at the outset that [START_REF] Lyddane | Small eccentricities or inclinations in the Brouwer theory of the artificial satellite[END_REF] and [START_REF] Gim | State transition matrix of relative motion for the perturbed noncircular reference orbit[END_REF] would yield the same overall degree of accuracy. Nevertheless, for the sake of completeness, and because, superficially, it is not apparent that the formulas of [START_REF] Gim | State transition matrix of relative motion for the perturbed noncircular reference orbit[END_REF] are mathematically equivalent to those systematized in Schaub and Junkins ( 2018), 4 we have also extended our numerical campaign to include a comparison between these as well (see Appendix C). Figure 2 shows a numerical confirmation of the validity of the Milankovitch formulation for satellites of the Sun-synchronous and Molniya type. The procedure was to first convert the initial osculating orbit into its corresponding mean elements using the developed formulas. These initial osculating and mean states were then propagated according to their dynamics described by Eqs. ( 29) and (30). In the former case, the results are equivalent to, though , 30 • ) (bottom) were converted into their corresponding mean elements using the developed Milankovitch formulation, and each set was propagated according to Eq. ( 29) (denoted "simulated osculating" dynamics in the legend) and Eq. (30) (denoted "simulated mean" dynamics), respectively. The osculating (cyan) and mean (yellow) trajectories were recovered (purple circles and red diamonds) at equal time steps using the aforementioned transformations, taking the respective simulated dynamics as input generally more accurate than, a simple Cowell integration in Cartesian space. The time histories of these evolutions at various subintervals were subsequently used as input to the respective osculating-to-mean and mean-to-osculating transformations in order to recover the aforementioned simulated trajectories.

Projecting both the simulated and recovered osculating evolutions into the radial, alongtrack, and cross-track frame, we use the norm root mean square (RMS) of the positional error as a means to assess the accuracy of the transformation. To keep the presentation brief, we do not consider the velocity errors or other metrics. Figure 3 shows the results of this process for Brouwer-Lyddane, Milankovitch, and the FFT transformations, respectively, for Sun-synchronous-like orbits and nearly critically inclined, semi-synchronous ones. The dependence of the resulting errors on the choice of orbit orientation angles is also highlighted by Fig. 3.

For completeness and further validation, Fig. 4 compares the aforementioned Brouwer-Lyddane transformation which includes the long-period terms, with a BL mean-to-osculating implementation that omits them. These test cases are simulated in the trusted semi-analytical were converted to their corresponding mean states using each respective transformation and subsequently propagated following Eq. 30. At every time step, the simulated mean trajectory was converted to osculating according to each transformation and compared against the simulated dynamics of Eq. 29. The norm RMS of the difference between the recovered and simulated positions (in km) over five orbital periods varied between 0.0556 and 20.9714 for Brouwer-Lyddane, 0.0666 and 0.3114 for Milankovitch, and 0.0533 and 20.9796 for the FFT transformation propagator, STELA, using a J 2 -only model. While slight discrepancies with STELA are to be expected due to different platforms, mean-element equations, and astronomical constants used in the respective simulations, the results are in good quantitative agreement. Importantly, the long-period terms in the full Brouwer theory do not cause appreciable changes in the recovered solutions throughout the whole of orbital phase space excepting a narrow band centered around the critical inclination.

Figure 5 shows error maps corresponding to two different initial eccentricities for 500 × 500 grids of initial inclinations and semi-major axes. Each grid point, together with (M, ω, Ω) = (0, 0, 0), was used to form the full osculating element state vector and propagated for five orbital periods. The same procedure outlined in the previous paragraph was used to characterize the accuracy of each transformation, where the colorbar in Fig. 5 , and where the initial mean anomaly, perigee, and node angles were all set to zero. The colorbar represents the norm RMS of the difference between the recovered and simulated positions over five orbital periods, according to each formulation. The colorbar limit of each map was set to the maximum error found in the Milankovitch scheme to provide a better contrast. Grid points leading to unphysical errors or to values that exceed this limit are represented in white. on the colorbar of each map in this numerical campaign so as to more clearly highlight the differences among the various transformations.

Figure 6 shows error maps in the semi-major axis-eccentricity plane for four different initial inclinations. For all three transformations, the error is highest at the boundary of allowable orbits, above which the perigee altitude would equal the Earth's radius.

The remaining slice of the action-like element space is given in Fig. 7, which shows how the errors manifest in the (i, e) plane for two values of initial semi-major axes (representing LEO satellites at 800 km altitude and GEO birds).

From this simulation campaign, we can conclude that the Milankovitch method consistently agrees with the classical Brouwer-Lyddane solution with each of the approaches 123 Fig. 7 Error maps in the inclination-eccentricity plane using the Brouwer-Lyddane (top), Milankovitch (middle), and numerical (bottom) transformations, respectively. Each panel samples an equidistant grid of 40 thousand initial osculating (i, a) values, for initial semi-major axes of R + 800 km (left) and a GEO (right), and where the initial mean anomaly, perigee, and node angles were all set to zero. The colorbar represents the norm RMS of the difference between the recovered and simulated positions over five orbital periods, according to each formulation. The colorbar limit of each map was set to the maximum error found in the Milankovitch scheme to provide a better contrast. Grid points leading to unphysical errors or to values that exceed this limit are represented in white. From left to right, the (maximum, mean) errors (in km) were (1.2992, 0.2601) and (94.1968, 1.9108) for Brouwer-Lyddane, (1.1378, 0.2973) and (3.0803, 0.0955) for Milankovitch, and (0.7508, 0.2230) and (45.3521, 1.7884) for the numerical scheme showing a slightly better accuracy in specific orbital regimes. Differences in positional error residuals remain in the order of a few kilometers or less.

Discussion

Our Milankovitch formulation is closed-form in both the eccentricity and angular momentum vectors. We have provided a general series solution for the mean longitude based on Hansen coefficients, which can be taken to any desired order of accuracy5 , and shown that our scheme performs in agreement with standard BL theories. We note that our formulation bypasses the critical inclination because we do not consider the long-periodic terms that arise from a second-order perturbation treatment. On this account, we also omitted these in Brouwer's theory in Fig. 4 to emphasize that the long-periodic terms are negligible away from the critical inclination in the time-scale of interest.

Being a non-canonical set of elements, our derivation followed the approach used by [START_REF] Kozai | The motion of a close Earth satellite[END_REF], as further elucidated by [START_REF] Scheeres | Orbital Motion in Strongly Perturbed Environments: Applications to Asteroid[END_REF]. We note, however, that, like [START_REF] Gim | State transition matrix of relative motion for the perturbed noncircular reference orbit[END_REF], we could have merely adopted Brouwer's generating function and computed the short-periodic corrections using Poisson-bracket operations. Nevertheless, we chose to present an independent derivation that, as a positive outcome, results to be applicable to other perturbations for which a convenient generating function is not always available, see [START_REF] Shen | Propagation and transformation of mean elements at geostationary orbits[END_REF]. Having provided a general Kozai-like scheme, not rooted in canonical perturbation theory, our approach can be potentially extended to non-conservative forces, such as solar radiation pressure and atmospheric drag, in addition to treating lunisolar third-body gravity and other predominant perturbations. Future work will be devoted to the inclusion of long-periodic terms in our transformation and to the application of this method to modeling non-conservative perturbations.

Conclusions

We have developed a new mean-to-osculating and inverse transformation based on the Milankovitch elements with the mean longitude as the fast variable, which is valid for all eccentricity values smaller than one. An extensive numerical campaign was used to validate the vectorial transformation over orbital phase-space grids tailored to the relative distribution of cataloged Earth satellites and debris.

A Auxiliary Formulas

We require the double quadrature of f , sin m f , cos m f , where m is an integer. The traditional approach has resorted to more or less elaborate expansions of various functions into sums of periodic terms, which depend chiefly on Taylor's series and Fourier's theorem, mainly because the integrals of these functions cannot be obtained conveniently in any other way [START_REF] Brumberg | Analytical Techniques of Celestial Mechanics[END_REF]. While detailed explicit tabulations for these quantities may be found in Cayley's and Newcomb's classical tables, these results may be deduced from the Hansen coefficients [START_REF] Hughes | The computation of tables of Hansen coefficients[END_REF]. Namely,

sin m f = ∞ k=1 S 0,m k (e) sin k M, ( 38a
) cos m f = C 0,m 0 (e) + ∞ k=1 C 0,m k (e) cos k M, ( 38b
)
where the coefficients S n,m k and C n,m k for integers k, n, and m are power series in the eccentricity, starting with degree |m -k|, and are related to Hansen coefficients X n,m k by

C n,m 0 = X n,m 0 , ( 39
) C n,m k = X n,m k + X n,m -k , ( 40
) S n,m k = X n,m k -X n,m -k . ( 41
)
123

For k = 0, exact analytical expressions exist for the zero-order Hansen coefficients X n,m 0 for all values of n and m; in particular, we recall a special result for the zero-order Hansen coefficients, derived by Kozai (1962a):

X 0,m 0 = 1 2π 2π 0 cos(m f ) dM = (-e) m (1 + m √ 1 -e 2 ) (1 + √ 1 -e 2 ) m , m = 1, 2, 3, . . . (42)
For k = 0, the analytical expressions for them do not terminate and, consequently, the series have to be truncated at some particular order in the eccentricity.

Hansen's coefficients can be expressed as series invoking Bessel and hypergeometric functions from which several formulae of recurrence can be derived that greatly facilitate their calculation [START_REF] Challe | Disturbing function and analytical solution of the problem of the motion of a satellite[END_REF][START_REF] Giacaglia | A note on Hansen's coefficients in satellite theory[END_REF]:

k X n,m k = ne 2η X n-1,m-1 k -X n-1,m+1 k + mηX n-2,m k , ( 43
)
(1 -e 2 ) 2 X n,m k = 1 - e 2 2 X n+2,m k + e(1 -e 2 ) X n+1,m+1 k + X n+1,m-1 k - e 2 4 X n+2,m+2 k + X n+2,m-2 k . ( 44
)
A useful series representation to seed the recursion, which dates back to Hansen, is given by [START_REF] Proulx | Series representations and rational approximations for Hansen coefficients[END_REF])

X n,m k = (1 + β 2 ) -(n+1) (-β) |k-m| ∞ =0 L (n-m--s-1) +s k 1 + β 2 L (a+m--t-1) +t k 1 + β 2 β 2 , ( 45
)
where the functions L ( j) i (x) are the generalized Laguerre polynomial defined as

L ( j) i (x) = ∞ =0 (-1) i + j i -l x ! , ( 46
)
and

s = [|k -m| + (k + m)] /2, ( 47
) t = [|k -m| -(k + m)] /2. ( 48
)
In the preceding expressions, η = √ 1e 2 and β = (1 -η)/e. Also needed in our derivation is the classical equation of the center, which is given by

f = M + 2 ∞ k=1 1 k J k (ke) + ∞ n=1 β n (J k-n (ke) + J k+n (ke)) sin k M, = M + ∞ k=1 φ k sin k M, ( 49
)
where J k (x) are Bessel functions of the first kind.

B Derivation of the short-periodic perturbations in the Milankovitch elements

The procedure involves first writing down an approximate differential equation that governs the short-period dynamics for each variable (e, H, l), according to Eq. ( 4). The first indefinite time integration of Eq. ( 5) requires the quadrature of dyadics and higher rank tensors of the dynamical variable r. In particular, we can expose all of the fast-variable terms of Eq. ( 13) and can complete the needed core integrals by changing the dependent variable from time to true anomaly or mean anomaly:

dt = r 2 H d f = 1 n dM. ( 50
)
Specifically, we have

1 r 3 dt = 1 H 1 r d f = μ H 3 (1 + e cos f ) d f = μ H 3 I 0 , ( 51
)
For the vector term, we have

I = r r 5 dt = 1 H r r 2 d f = μ 2 H 5 (1 + e cos f ) 2 (cos f ê + sin f ê⊥ ) d f = μ 2 H 5 (I 1 ê + I 2 ê⊥ ). ( 52
)
Similarly, for the dyadic terms, we find

I I = r r r 5 dt = 1 H r r r d f = μ H 3 (1 + e cos f ) cos 2 f êê + cos f sin f (êê ⊥ + ê⊥ ê) + sin 2 f ê⊥ ê⊥ d f = μ H 3 I I 11 êê + I I 12 (êê ⊥ + ê⊥ ê) + I I 22 ê⊥ ê⊥ , ( 53
)
where

I I 11 = (1 + e cos f ) cos 2 f d f = 1 12 6 f + 9e sin f + 3 sin 2 f + e sin 3 f , ( 54
)
I I 12 = (1 + e cos f ) cos f sin f d f = - 1 12 3e cos f + 3 cos 2 f + e cos 3 f , ( 55
)
I I 22 = (1 + e cos f ) sin 2 f d f 123 = 1 12 6 f + 3e sin f -3 sin 2 f -e sin 3 f . ( 56
)
In the same way, we have for the first triadic terms: Finally, for the last triadic term, we have

I I I = r r r r 6 dt = 1 H r r r r d f = μ H 3 (1 + e cos f ) cos 3 f êêê + cos 2 f sin f (êêê ⊥ + êê ⊥ ê + ê⊥ êê) + cos f sin 2 f (êê ⊥
I V = r r r r 7 dt = 1 H r r r r 2 d f = μ 2 H 5 (1 + e cos f ) 2 cos 3 f êêê + cos 2 f sin f (êêê ⊥ + êê ⊥ ê + ê⊥ êê) + cos f sin 2 f (êê ⊥ ê⊥ + ê⊥ êê ⊥ + ê⊥ ê⊥ ê) + sin 3 f ê⊥ ê⊥ ê⊥ d f = μ 2 H 5 I V 111 êêê + I V 112 (êêê ⊥ + êê ⊥ ê + ê⊥ êê) + I V 122 (êê ⊥ ê⊥ + ê⊥ êê ⊥ + ê⊥ ê⊥ ê) + I V 222 ê⊥ ê⊥ ê⊥ . ( 58
)
where

I 0 = (1 + e cos) d f = f + e sin f , (59a)
I 1 = (1 + e cos f ) 2 cos f d f = 1 12 12e f + (12 + 9e 2 ) sin f + 6e sin 2 f + e 2 sin 3 f , ( 59b
)
I 2 = (1 + e cos f ) 2 sin f d f = - 1 12 (12 + 3e 2 ) cos f + 6e cos 2 f + e 2 cos 3 f , ( 59c
)
I I I 111 = (1 + e cos f ) cos 3 f d f = 1 96 36e f + 72 sin f + 24e sin 2 f + 8 sin 3 f + 3e sin 4 f , ( 59d
)
I I I 112 = (1 + e cos f ) cos 2 f sin f d f = - 1 96 24 cos f + 12e cos 2 f + 8 cos 3 f + 3e cos 4 f , ( 59e
)
I I I 122 = (1 + e cos f ) cos f sin 2 f d f 123 = 1 96 12e f + 24 sin f -8 sin 3 f -3e sin 4 f , ( 59f
)
I I I 222 = (1 + e cos f ) sin 3 f d f = - 1 96 72 cos f + 12e cos 2 f -8 cos 3 f -3e cos 4 f , ( 59g
)
I V 111 = (1 + e cos f ) 2 cos 3 f d f = 1 240 180e f + 30(6 + 5e 2 ) sin f + 120e sin 2 f + 5(4 + 5e 2 ) sin 3 f + 15e sin 4 f + 3e 2 sin 5 f , ( 59h
)
I V 112 = (1 + e cos f ) 2 cos 2 f sin f d f = - 1 240 30(2 + e 2 ) cos f + 60e cos 2 f + 5(4 + 3e 2 ) cos 3 f + 15e cos 4 f + 3e 2 cos 5 f , ( 59i
)
I V 122 = (1 + e cos f ) 2 cos f sin 2 f d f = 1 240 60e f + 30(2 + e 2 ) sin f -5(4 + e 2 ) sin 3 f -15e sin 4 f -3e 2 sin 5 f , (59j)
I V 222 = (1 + e cos f ) 2 sin 3 f d f = - 1 240 30(6 + e 2 ) cos f + 60e cos 2 f -5(4 -e 2 ) cos 3 f -15e cos 4 f -3e 2 cos 5 f . ( 59k
)
We also require the averages of these core integrals, which can easily be deduced by exploiting the Hansen coefficients (Eq. ( 38)) and equation of the center (Eq. ( 49)) expansions:

I 0 = 1 2π 2π 0 ( f + e sin f ) dM = π, ( 60a
)
I 1 = 1 12 1 2π 2π 0 12e f + (12 + 9e 2 ) sin f + 6e sin 2 f + e 2 sin 3 f dM = πe, ( 60b
)
I 2 = - 1 12 1 2π 2π 0 (12 + 3e 2 ) cos f + 6e cos 2 f + e 2 cos 3 f dM = - 1 12 (12 + 3e 2 )X 0,1 0 + 6eX 0,2 0 + e 2 X 0,3 0 , ( 60c
)
I I 11 = 1 12 1 2π 2π 0 6 f + 9e sin f + 3 sin 2 f + e sin 3 f dM = 1 2 π, ( 60d
)
123

I I 12 = - 1 1 2π 2π 0 3e cos f + 3 cos 2 f + e cos 3 f dM = - 1 3eX 0,1 0 + 3X 0,2 0 + eX 0,3 0 , ( 60e
)
I I 22 = 12 1 2π 2π 0 6 f + 3e sin f -3 sin 2 f -e sin 3 f dM = 1 2 π, ( 60f
)
I I I 111 = 96 1 2π 2π 0 36e f + 72 sin f + 24e sin 2 f + 8 sin 3 f + 3e sin 4 f dM = 3 8 πe, ( 60g
)
I I I 112 = - 1 1 2π 2π 0 24 cos f + 12e cos 2 f + 8 cos 3 f + 3e cos 4 f dM = - 1 24X 0,1 0 + 12eX 0,2 0 + 8X 0,3 0 + 3eX 0,4 0 , ( 60h
)
I I I 122 = 96 1 2π 2π 0 12e f + 24 sin f -8 sin 3 f -3e sin 4 f dM = 1 8 πe, (60i)
I I I 222 = - 1 1 2π 2π 0 72 cos f + 12e cos 2 f -8 cos 3 f -3e cos 4 f dM = - 1 72X 0,1 0 + 12eX 0,2 0 -8X 0,3 0 -3eX 0,4 0 , ( 60j
)
I V 111 = 1 240 1 2π 2π 0 180e f + 30(6 + 5e 2 ) sin f + 120e sin 2 f + 5(4 + 5e 2 ) sin 3 f + 15e sin 4 f + 3e 2 sin 5 f dM = πe, (60k)
I V 112 = - 1 1 2π 2π 0 30(2 + e 2 ) cos f + 60e cos 2 f + 5(4 + 3e 2 ) cos 3 f + 15e cos 4 f + 3e 2 cos 5 f dM = - 1 30(2 + e 2 )X 0,1 0 + 60eX 0,2 0 + 5(4 + 3e 2 )X 0,3 0 + 15eX 0,4 0 + 3e 2 X 0,5 0 , ( 60l
)
I V 122 = 1 240 1 2π 2π 0 60e f + 30(2 + e 2 ) sin f -5(4 + e 2 ) sin 3 f -15e sin 4 f -3e 2 sin 5 f dM = 1 4 πe, ( 60m
)
I V 222 = - 1 1 2π 2π 0 30(6 + e 2 ) cos f + 60e cos 2 f 123 -5(4 -e 2 ) cos 3 f -15e cos 4 f -3e 2 cos 5 f dM = - 1 240 30(6 + e 2 )X 0,1 0 + 60eX 0,2 0 -5(4 -e 2 )X 0,3 0 -15eX 0,4 0 -3e 2 X 0,5 0 . ( 60n
)
The previous results are sufficient to derive the short-periodic perturbations in the eccentricity and angular momentum vectors, according to Eqs. ( 18a) and (18b), respectively. The mean longitude, however, requires the computation of additional integrals due to the mean motion Taylor-series terms, as detailed in Eqs. ( 22) and ( 23). These can be carried through to the core integrals, where again we can use Eqs. ( 38) and (49). In particular, using the mean anomaly as the dependent variable, we have

I 1 = 1 12 12e f + (12 + 9e 2 ) sin f + 6e sin 2 f + e 2 sin 3 f dM = 1 12 6eM 2 - ∞ k=1 1 k 12eφ k + (12 + 9e 2 )S 0,1 k + 6eS 0,2 k + e 2 S 0,3 k cos k M , ( 61a
)
I 2 = - 1 12 (12 + 3e 2 ) cos f + 6e cos 2 f + e 2 cos 3 f dM = - 1 12 (12 + 3e 2 )C 0,1 0 + 6eC 0,2 0 + e 2 C 0,3 0 M + ∞ k=1 1 k (12 + 3e 2 )C 0,1 k + 6eC 0,2 k + e 2 C 0,3 k sin k M , (61b)
I I 11 = 1 12 6 f + 9e sin f + 3 sin 2 f + e sin 3 f dM = 1 12 3M 2 - ∞ k=1 1 k 6φ k + 9eS 0,1 k + 3S 0,2 k + eS 0,3 k cos k M , ( 61c
)
I I 12 = - 1 12 3e cos f + 3 cos 2 f + e cos 3 f dM = - 1 12 3eC 0,1 0 + 3C 0,2 0 + eC 0,3 0 M + ∞ k=1 1 k 3eC 0,1 k + 3C 0,2 k + eC 0,3 k sin k M , ( 61d
)
I I 22 = 1 12 6 f + 3e sin f -3 sin 2 f -e sin 3 f dM = 1 12 3M 2 - ∞ k=1 1 k 6φ k + 3eS 0,1 k -3S 0,2 k -eS 0,3 k cos k M , ( 61e
)
I I I 111 = 1 96 36e f + 72 sin f + 24e sin 2 f + 8 sin 3 f + 3e sin 4 f dM = 1 96 18eM 2 - ∞ k=1 1 k 36eφ k + 72S 0,1 k + 24eS 0,2 k + 8S 0,3 k + 3eS 0,4 k cos k M , ( 61f
)
123

I I I = - 1 96 24 cos f + 12e cos 2 f + 8 cos 3 f + 3e cos 4 f dM = - 1 96 24C 0,1 0 + 12eC 0,2 0 + 8C 0,3 0 + 3eC 0,4 0 M + ∞ k=1 1 k 24C 0,1 k + 12eC 0,2 k + 8C 0,3 k + 3eC 0,4 k sin k M , ( 61g
)
I I I = 1 96 12e f + 24 sin f -8 sin 3 f -3e sin 4 f dM = 1 96 6eM 2 - ∞ k=1 1 k 12eφ k + 24S 0,1 k -8S 0,3 k -3eS 0,4 k cos k M , ( 61h
)
I I I = - 1 96 72 cos f + 12e cos 2 f -8 cos 3 f -3e cos 4 f dM = - 1 96 72C 0,1 0 + 12eC 0,2 0 -8C 0,3 0 -3eC 0,4 0 M + ∞ k=1 1 k 72C 0,1 k + 12eC 0,2 k -8C 0,3 k -3eC 0,4 k sin k M , ( 61i
)
I V = 1 240 1 2π 2π 0 180e f + 30(6 + 5e 2 ) sin f + 120e sin 2 f + 5(4 + 5e 2 ) sin 3 f + 15e sin 4 f + 3e 2 sin 5 f dM = 1 240 90eM 2 - ∞ k=1 1 k 180eφ k + 30(6 + 5e 2 )S 0,1 k + 120eS 0,2 k + 5(4 + 5e 2 )S 0,3 k + 15eS 0,4 k + 3e 2 S 0,5 k cos k M , ( 61j
)
I V = - 1 240 30(2 + e 2 ) cos f + 60e cos 2 f + 5(4 + 3e 2 ) cos 3 f + 15e cos 4 f + 3e 2 cos 5 f dM = - 1 240 30(2 + e 2 )C 0,1 0 + 60eC 0,2 0 + 5(4 + 3e 2 )C 0,3 0 + 15eC 0,4 0 + 3e 2 C 0,5 0 M + ∞ k=1 1 k 30(2 + e 2 )C 0,1 k + 60eC 0,2 k + 5(4 + 3e 2 )C 0,3 k + 15eC 0,4 k + 3e 2 C 0,5 k sin k M , ( 61k
) I V = 1 240 60e f + 30(2 + e 2 ) sin f -5(4 + e 2 ) sin 3 f -15e sin 4 f -3e 2 sin 5 f dM = 1 240 30eM 2 ∞ k=1 1 k 60eφ k + 30(2 + e 2 )S 0,1 k 123 -5(4 + e 2 )S 0,3 k -15eS 0,4 k -3e 2 S 0,5 k cos k M , ( 61l
)
I V 222 = - 1 240 30(6 + e 2 ) cos f + 60e cos 2 f -5(4 -e 2 ) cos 3 f -15e cos 4 f -3e 2 cos 5 f dM = - 1 240 30(6 + e 2 )C 0,1 0 + 60eC 0,2 0 -5(4 -e 2 )C 0,3 0 -15eC 0,4 0 -3e 2 C 0,5 0 M + ∞ k=1 1 k 30(6 + e 2 )C 0,1 k + 60eC 0,2 k -5(4 -e 2 )C 0,3 k -15eC 0,4 k -3e 2 C 0,5 k sin k M . ( 61m
)
Note that integrands are the same as those of the previous averages I 1 , I 2 , I I 11 , . . ., in Eq. ( 60), but they differ in that the former are definite while the later are indefinite. Finally, we require the averages of these doubly-integrated expressions:

I 1 = 2 3 π 2 e, ( 62a
)
I 2 = - 1 12 (12 + 3e 2 )C 0,1 0 + 6eC 0,2 0 + e 2 C 0,3 0 π , ( 62b
)
I I 11 = 1 3 π 2 , ( 62c
)
I I 12 = -1 12 3eC 0,1 0 + 3C 0,2 0 + eC 0,3 0 π , (62d)

I I 22 = 1 3 π 2 , ( 62e
)
I I I 111 = 1 4 π 2 e, ( 62f
)
I I I 112 = - 1 96 24C 0,1 0 + 12eC 0,2 0 + 8C 0,3 0 + 3eC 0,4 0 π , ( 62g
)
I I I 122 = 1 12 π 2 e, ( 62h
)
I I I 222 = - 1 96 72C 0,1 0 + 12eC 0,2 0 -8C 0,3 0 -3eC 0,4 0 π , ( 62i
)
I V 111 = 1 2 π 2 e, ( 62j
)
I V 112 = - 1 240 30(2 + e 2 )C 0,1 0 + 60eC 0,2 0 + 5(4 + 3e 2 )C 0,3 0 + 15eC 0,4 0 + 3e 2 C 0,5 0 π , ( 62k
)
I V 122 = 1 6 π 2 e, ( 62l
)
I V 222 = - 1 240
30(6 + e 2 )C 0,1 0 + 60eC 0,2 0 -5(4e 2 )C 0,3 0 -15eC 0,4 0 -3e C 0,5 0 π .

(62m)

With these ingredients at hand, we now step through the derivation of the short-period corrections. Eccentricity vector: Subtracting Eq. (30a) from Eq. (29a), holding the slowly varying orbital elements constant, and integrating, according to Eq. (17a), gives the short-periodic perturbations in e as

de sp (t) = 3J 2 R 2 2 p 2 I 1 ê⊥ + 2(ê • p) ĥ • p -I 2 ê -2(ê ⊥ • p) ĥ • p -5I V 111 (ê • p) 2 ê⊥ + 5I V 112 (ê • p) 2 ê -2(ê • p)(ê ⊥ • p)ê ⊥ + 5I V 122 2(ê • p)(ê ⊥ • p)ê -(ê ⊥ • p) 2 ê⊥ + 5I V 222 (ê ⊥ • p) 2 ê -2eI I 11 (ê • p)(ê ⊥ • p)ê + 2eI I 12 (ê • p) p -(ê • p)(ê ⊥ • p)ê ⊥ -(ê ⊥ • p) 2 ê + 2eI I 22 (ê ⊥ • p) p -(ê ⊥ • p) 2 ê⊥ + 2 (ê • p) 2 -(ê ⊥ •
+ 1 2 M 1 -5( ĥ • p) 2 ĥ + 2( ĥ • p) p • e , ( 63
)
The mean value of this deviation with respect to the mean anomaly is given by π 1 -5( ĥ • p) 2 ĥ + 2( ĥ • p) p • e , (

de sp = 3J 2 R 2 2 p 2 I 1 ê⊥ + 2(ê • p) ĥ • p -I 2 ê -2(ê ⊥ • p) ĥ • p -5I V 111 (ê • p) 2 ê⊥ + 5I V 112 (ê • p) 2 ê -2(ê • p)(ê ⊥ • p)ê ⊥ + 5I V 122 2(ê • p)(ê ⊥ • p)ê -(ê ⊥ • p) 2 ê⊥ + 5I V 222 (ê ⊥ • p) 2 ê -2eI I 11 (ê • p)(ê ⊥ • p)ê + 2eI I 12 (ê • p) p -(ê • p)(ê ⊥ • p)ê ⊥ -(ê ⊥ • p) 2 ê + 2eI I 22 (ê ⊥ • p) p -(ê ⊥ • p) 2 ê⊥ + 2 (ê • p) 2 -(ê ⊥ •
We then have from Eq. (18a), after simplifying, that

e sp (t) = 3J 2 R 2 2 p 2 I 1 ê⊥ + 2(ê • p) ĥ • p -I 2 ê -2(ê ⊥ • p) ĥ • p -5I V 111 (ê • p) 2 ê⊥ + 5 I V 112 (ê • p) 2 ê -2(ê • p)(ê ⊥ • p)ê ⊥ + 5I V 122 2(ê • p)(ê ⊥ • p)ê -(ê ⊥ • p) 2 ê⊥ + 5 I V 222 (ê ⊥ • p) 2 ê -2eI I 11 (ê • p)(ê ⊥ • p)ê + 2e I I 12 (ê • p) p -(ê • p)(ê ⊥ • p)ê ⊥ -(ê ⊥ • p) 2 ê + 2eI I 22 (ê ⊥ • p) p -(ê ⊥ • p) 2 ê⊥ + 2 (ê • p) 2 -(ê ⊥ •
I V 112 = I V 112 -I V 112 , ( 66e
)
I V 222 = I V 222 -I V 222 . (66f)
Equation ( 65) can be cast into a simpler form that exposes the perifocal frame components, as given in Eq. (32a).

Angular momentum vector: In the same vein, following Eq. ( 17b), upon integrating, we have The result is then

d H sp (t) = - 3H J 2 R
H sp (t) = - 3H J 2 R 2 p 2 p • I I 11 êê + I I 12 (êê ⊥ + ê⊥ ê) + I I 22 ê⊥ ê⊥ + 1 2 M ĥ ĥ • p. (69)
Equation ( 69) can also be written in a non-dyadic form as shown in Eq. (32b) Mean longitude: We first note that

g l (x, t) -g l (x) dt = 3J 2 R 2 2 p 2 e 1 + √ 1 -e 2 I 1 -3 I V 111 (ê • p) 2 + 2I V 112 (ê • p)(ê ⊥ • p) + I V 122 (ê ⊥ • p) 2 + 2 (I I I 222 -I I I 112 + I V 222 -I V 112 ) (ê • p)(ê ⊥ • p) + 2 (I I I 122 + I V 122 ) (ê • p) 2 -(ê ⊥ • p) 2 -2 3 1 -e 2 + ĥ • p 1 + ĥ • p I I 11 (ê • p) 2 + 2I I 12 (ê • p)(ê ⊥ • p) + I I 22 (ê ⊥ • p) 2 + 2 1 -e 2 I 0 - 1 2 M 1 -e 2 3( ĥ • p) 2 -1 + 5( ĥ • p) 2 -2( ĥ • p) -1 , ( 70
) and 1 T T 0 g l (x, t) -g l (x) dt 2 = 3J 2 R 2 2 p 2 e 1 + √ 1 -e 2 I 1 -3 I V 111 (ê • p) 2 + 2I V 112 (ê • p)(ê ⊥ • p) + I V 122 (ê ⊥ • p) 2 + 2 I I I 222 -I I I 112 + I V 222 -I V 112 (ê • p)(ê ⊥ • p) + 2 I I I 122 + I V 122 (ê • p) 2 -(ê ⊥ • p) 2 -2 3 1 -e 2 + ĥ • p 1 + ĥ • p I I 11 (ê • p) 2 + 2I I 12 (ê • p)(ê ⊥ • p) + I I 22 (ê ⊥ • p) 2 123 + 2 1 -e 2 I 0 - 1 2 π 1 -e 2 3( ĥ • p) 2 -1 + 5( ĥ • p) 2 -2( ĥ • p) -1 . (71)
Furthermore,

e sp (τ ) dτ = 3J 2 R 2 2np 2 I 1 ê⊥ + 2(ê • p) ĥ • p -( I 2 -I 2 M) ê -2(ê ⊥ • p) ĥ • p -5 I V 111 (ê • p) 2 ê⊥ + 5( I V 112 -I V 112 M) (ê • p) 2 ê -2(ê • p)(ê ⊥ • p)ê ⊥ + 5 I V 122 2(ê • p)(ê ⊥ • p)ê -(ê ⊥ • p) 2 ê⊥ + 5( I V 222 -I V 222 M)(ê ⊥ • p) 2 ê -2e I I 11 (ê • p)(ê ⊥ • p)ê + 2e( I I 12 -I I 12 M) (ê • p) p -(ê • p)(ê ⊥ • p)ê ⊥ -(ê ⊥ • p) 2 ê + 2e I I 22 (ê ⊥ • p) p -(ê ⊥ • p) 2 ê⊥ + 2 (ê • p) 2 -(ê ⊥ •
+ 1 4 M 2 1 -5( ĥ • p) 2 ĥ + 2( ĥ • p) p • e , ( 72
) H sp (t) dt = - 3H J 2 R 2 np 2 p • I I 11 êê + ( I I 12 -I I 12 M)(êê ⊥ + ê⊥ ê) + I I 22 ê⊥ ê⊥ + 1 4 M 2 ĥ ĥ • p, (73)
We also have that

1 T T 0 e sp (t) dt 2 = 3J 2 R 2 2np 2 I 1 ê⊥ + 2(ê • p) ĥ • p -( I 2 -I 2 π) ê -2(ê ⊥ • p) ĥ • p -5 I V 111 (ê • p) 2 ê⊥ + 5( I V 112 -I V 112 π) (ê • p) 2 ê -2(ê • p)(ê ⊥ • p)ê ⊥ + 5 I V 122 2(ê • p)(ê ⊥ • p)ê -(ê ⊥ • p) 2 ê⊥ + 5( I V 222 -I V 222 π)(ê ⊥ • p) 2 ê -2e I I 11 (ê • p)(ê ⊥ • p)ê + 2e( I I 12 -I I 12 π) (ê • p) p -(ê • p)(ê ⊥ • p)ê ⊥ -(ê ⊥ • p) 2 ê + 2e I I 22 (ê ⊥ • p) p -(ê ⊥ • p) 2 ê⊥ + 2 (ê • p) 2 -(ê ⊥ •
+ 1 3 π 2 1 -5( ĥ • p) 2 ĥ + 2( ĥ • p) p • e , ( 74
) 1 T T 0 H sp (t) dt 2 = - 3H J 2 R 2 np 2 p • I I 11 êê + ( I I 12 -I I 12 π)(êê ⊥ + ê⊥ ê) + I I 22 ê⊥ ê⊥ + 1 3 π 2 ĥ ĥ • p, (75)
123

The short-periodic perturbations for mean longitude can thus be stated as 4π 2 ) 1 -5( ĥ • p) 2 ĥ + 2( ĥ

l sp (t) = dl sp (t) -dl sp (t) = 3J 2 R 2 2 p 2 e 1 + √ 1 -e 2 (I 1 -I 1 ) -3 (I V 111 -I V 111 )(ê • p) 2 + 2(I V 112 -I V 112 )(ê • p)(ê ⊥ • p) + (I V 122 -I V 122 )(ê ⊥ • p) 2 + 2 I I I 222 -I I I 222 -I I I 112 + I I I 112 + I V 222 -I V 222 -I V 112 + I V 112 (ê • p)(ê ⊥ • p) + 2 I I I 122 -I I I 122 + I V 122 -I V 122 (ê • p) 2 -(ê ⊥ • p) 2 -2 3 1 -e 2 + ĥ • p 1 + ĥ • p (I I 11 -I I 11 )(ê • p) 2 + 2(I I 12 -I I 12 )(ê • p)(ê ⊥ • p) + (I I 22 -I I 22 )(ê ⊥ • p) 2 + 2 1 -e 2 (I 0 -I 0 ) - 1 2 (M -π) 1 -e 2 3( ĥ • p) 2 -1 + 5( ĥ • p) 2 -2( ĥ • p) -1 - 3J 2 R 2 2np 2 3μ 2 H 3 1 -e 2 ( I 1 -I 1 ) ê⊥ + 2(ê • p) ĥ • p -( I 2 -I 2 -I 2 (M -π)) ê -2(ê ⊥ • p) ĥ • p -5( I V 111 -I V 111 )(ê • p) 2 ê⊥ + 5( I V 112 -I V 112 -I V 112 (M -π)) (ê • p) 2 ê -2(ê • p)(ê ⊥ • p)ê ⊥ + 5( I V 122 -I V 122 ) 2(ê • p)(ê ⊥ • p)ê -(ê ⊥ • p) 2 ê⊥ + 5( I V 222 -I V 222 -I V 222 (M -π ))(ê ⊥ • p) 2 ê -2e( I I 11 -I I 11 )(ê • p)(ê ⊥ • p)ê + 2e( I I 12 -I I 12 -I I 12 (M -π)) (ê • p) p -(ê • p)(ê ⊥ • p)ê ⊥ -(ê ⊥ • p) 2 ê + 2e( I I 22 -I I 22 ) (ê ⊥ • p) p -(ê ⊥ • p) 2 ê⊥ + 2 (ê • p) 2 -(ê ⊥ •
• p) p • e • e + 3H J 2 R 2 np 2 3μ 2 H 5 (1 -e 2 ) 3/2 p • ( I I 11 -I I 11 )êê + ( I I 12 -I I 12 -I I 12 (M -π ))(êê ⊥ + ê⊥ ê) + ( I I 22 -I I 22 )ê ⊥ ê⊥ + 1 12 (3M 2 -4π 2 ) ĥ ĥ • p • H. ( 76
)
Simplifying gives

l sp (t) = 3J 2 R 2 2 p 2 e 1 + √ 1 -e 2 I 1 -3 I V 111 (ê • p) 2 + 2 I V 112 (ê • p)(ê ⊥ • p) + I V 122 (ê ⊥ • p) 2 + 2 I I I 222 -I I I 112 + I V 222 -I V 112 (ê • p)(ê ⊥ • p) + 2 (I I I 122 + I V 122 ) (ê • p) 2 -(ê ⊥ • p) 2 -2 3 1 -e 2 + ĥ • p 1 + ĥ • p I I 11 (ê • p) 2 + 2 I I 12 (ê • p)(ê ⊥ • p) + I I 22 (ê ⊥ • p) 2 + 2 1 -e 2 I 0 - 1 2 M 1 -e 2 3( ĥ • p) 2 -1 + 5( ĥ • p) 2 -2( ĥ • p) -1 + 3eH np 2 1 -e 2 2 IV + e II (ê • p)(ê ⊥ • p) + I 2 1 + 2(ê ⊥ • p) 2 -5 I V 112 (ê • p) 2 -5 I V 222 (ê ⊥ • p) 2 -2 I I I 112 + e I I 12 (ê • p) 2 -(ê ⊥ • p) 2 + 6H np 2 (1 -e 2 ) 3/2 II(ê • p)(ê ⊥ • p) -I I 12 (ê • p) 2 -(ê ⊥ • p) 2 , ( 77
)
123

where (78g)

I 2 = I 2 -I 2 -I 2 (M -π) = - 1 12 ∞ k=1 1 k (12 + 3e 2 )C 0,1 k + 6eC 0,2 k + e 2 C 0,3 k sin k M, ( 78a

C Numerical comparisons between Brouwer-Lyddane and Gim-Alfriend

Figures 8,9,10,11 show comparisons between the Brouwer-Lyddane and Gim-Alfriend formulations, as concerns the error maps of Sect. 4.

Fig. 8

Error maps in the inclination-semi-major axis plane using the Brouwer-Lyddane (top panels) and Gim-Alfriend (bottom panels) transformations, respectively. Each panel samples an equidistant grid of initial osculating (i, a) values, for initial eccentricities of 0.01 (top-left), 0.2 (top-right), 0.69 (bottom-left), and 0.74 (bottom-right), and where the initial mean anomaly, perigee, and node angles were all set to zero. The colorbar represents the norm RMS of the difference between the recovered and simulated positions over five orbital periods, according to each formulation. The colorbar limit was set to the maximum error found between the Milankovitch and FFT schemes (cf. Figs. 5), as the Brouwer-Lyddane formulation becomes singular near the critical inclination. Grid points leading to unphysical errors or to values that exceed this limit are represented in white. Clockwise starting from the top left: Brouwer-Lyddane recorded (maximum, mean) errors (in km) of (1.0626, 0.1967), (68.9371, 0.1804), (6.7332, 4.3022), and (1689.0290, 19.5805); and Gim-Alfriend recorded errors of (1.3476, 0.1967), (65.1025, 0.0728), (3.4344, 0.0914), and (1470.5632, 18.7458) 123 Fig. 9 Error maps in the semi-major axis-eccentricity plane using the Brouwer-Lyddane (top panels) and Gim-Alfriend (bottom panels) transformations, respectively. Each panel samples an equidistant grid of initial osculating (a, e) values, for initial inclinations of 6 • (top-left), 63 • (top-right), 98 • (bottom-left), and 116.6 • (bottom-right), and where the initial mean anomaly, perigee, and node angles were all set to zero. The colorbar represents the norm RMS of the difference between the recovered and simulated positions over five orbital periods, according to each formulation. The colorbar limit of each map was set to the maximum error found in the Milankovitch scheme (cf. Fig. 6) to provide a better contrast. Grid points leading to unphysical errors or to values that exceed this limit are represented in white. Clockwise starting from the top left: Brouwer-Lyddane recorded (maximum, mean) errors (in km) of (58.5237, 1.2662), (59.5538, 1.2013), (59.6549, 1.2698), and (82198.6187, 218.0410);and Gim-Alfriend recorded errors of (58.4994, 1.0281), (59.8837, 1.6868), (59.6433, 1.0878), and (48982.3446, 196.3659) Fig. 10 Error maps in the inclination-eccentricity plane using the Brouwer-Lyddane (top) and Gim-Alfriend (bottom) transformations, respectively. Each panel samples an equidistant grid of initial osculating (i, a) values, for initial semi-major axes of R + 800 km (left) and a GEO (right), and where the initial mean anomaly, perigee, and node angles were all set to zero. The colorbar represents the norm RMS of the difference between the recovered and simulated positions over five orbital periods, according to each formulation. The colorbar limit of each map was set to the maximum error found in the Milankovitch scheme (cf. Fig. 7) to provide a better contrast. Grid points leading to unphysical errors or to values that exceed this limit are represented in white. From left to right, the (maximum, mean) errors (in km) were (1.2992, 0.2601) and (94.1968, 1.9108) for Brouwer-Lyddane, and (2.5030, 0.2355) and (80.6873, 1.8132) for Gim-Alfriend Fig. 11 Error maps in the mean anomaly-periapsis angle plane (left) and periapsis-node angles plane (right) using the Brouwer-Lyddane (top) and Gim-Alfriend (bottom) transformations, respectively. Each panel samples an equidistant grid of initial osculating (M, ω) or (ω, Ω) values, for an initial semi-major axis of R + 800 km, eccentricity of 0.1, and inclination of 98 • . The remaining orbital element needed to form the full state vector was set to zero in each case. The colorbar represents the norm RMS of the difference between the recovered and simulated positions over five orbital periods, according to each formulation. The colorbar limit was set to the maximum error found across all schemes as there were no singular or outlier cases for these maps. From left to right, the (maximum, mean) errors (in km) were (1.3241, 0.4100) and (1.0332, 0.4704) for Brouwer-Lyddane, and (1.0173, 0.2060) and (1.3206, 0.4133) for Gim-Alfriend 123

Fig. 1

 1 Fig.1Schematic showing a comparison between the numerical and mean solutions for an arbitrary orbital element, when starting from the same initial condition. The mean solution changes approximately linearly over one orbit and is referred to as the averaged variation. The short-period oscillations are the fluctuations that happen per orbit in the real motion

Fig. 3

 3 Fig.3Radial, along-track, and cross-track errors between the recovered and simulated osculating trajectories, using the Brouwer-Lyddane (orange, dash-dot), Milankovitch (blue, dashed), and FFT (gray, solid) transformations, respectively. The initial osculating orbits (a, e, i) = (R + 800 km, 0.001, 98 • ) (top) and (26562 km, 0.75, 63 • ) (bottom), each with (ω, Ω) = (90 • , 180 • ) and M = 0 (left) or M = 45 • (right), were converted to their corresponding mean states using each respective transformation and subsequently propagated following Eq. 30. At every time step, the simulated mean trajectory was converted to osculating according to each transformation and compared against the simulated dynamics of Eq. 29. The norm RMS of the difference between the recovered and simulated positions (in km) over five orbital periods varied between 0.0556 and 20.9714 for Brouwer-Lyddane, 0.0666 and 0.3114 for Milankovitch, and 0.0533 and 20.9796 for the FFT transformation

Fig. 4 Fig. 5 Fig. 6

 456 Fig. 4 Radial, along-track, and cross-track errors between the recovered and simulated osculating trajectories, using the full Brouwer-Lyddane (orange, dash-dot), BL without long-period terms (blue, dashed), and STELA (black, solid) transformations, respectively. The initial osculating orbits (a, e, i) = (R + 800 km, 0.001, 98 • ) (left) and (26562 km, 0.75, 63 • ) (right), each with (ω, Ω) = (90 • , 180 • ) and M = 45 (left) or M = 0 • (right), were converted to their corresponding mean states using each respective transformation and subsequently propagated. At every time step, the simulated mean trajectory was converted to osculating according to each transformation and compared against the simulated osculating dynamics. The norm RMS of the difference between the recovered and simulated positions (in km) over five orbital periods was 0.0556 and 20.9714 for BL, 0.0555 and 20.9785 for BL w/out LP, and 0.0543 and 22.4028 for the STELA transformation

  ê⊥ + ê⊥ êê ⊥ + ê⊥ ê⊥ ê) + sin 3 f ê⊥ ê⊥ ê⊥ d f = μ H 3 I I I 111 êêê + I I I 112 (êêê ⊥ + êê ⊥ ê + ê⊥ êê) + I I I 122 (êê ⊥ ê⊥ + ê⊥ êê ⊥ + ê⊥ ê⊥ ê) + I I I 222 ê⊥ ê⊥ ê⊥ . (57)

  p) 2 I I I 112 ê + I I I 122 ê⊥ + 2(ê • p)(ê ⊥ • p) (I I I 122 -I I I 111 ) ê + (I I I 222 -I I I 112 ) ê⊥

  p) 2 I I I 112 ê + I I I 122 ê⊥ + 2(ê • p)(ê ⊥ • p) I I I 122 -I I I 111 ê + I I I 222 -I I I 112 ê⊥ + 1 2

  p) 2 I I I 112 ê + I I I 122 ê⊥ + 2(ê • p)(ê ⊥ • p) (I I I 122 -I I I 111 ) ê + 2(ê • p)(ê ⊥ • p) I I I 222 -I I I 112 ê⊥ + 1 2 M 1 -5( ĥ • p) 2 ĥ + 2( ĥ • p) p • e , = I I 12 -I I 12 , (66b) I I I 112 = I I I 112 -I I I 112 , (66c) I I I 222 = I I I 222 -I I I 222 , (66d)

  p) 2 ( I I I 112 -I I I 112 M)ê + I I I 122 ê⊥ + 2(ê • p)(ê ⊥ • p) I I I 122 -I I I 111 ê + I I I 222 -I I I 222 M -I I I 112 + I I I 112 M ê⊥

  p) 2 ( I I I 112 -I I I 112 π)ê + I I I 122 ê⊥ + 2(ê • p)(ê ⊥ • p) I I I 122 -I I I 111 ê + I I I 222 -I I I 222 π -I I I 112 + I I I 112 π ê⊥

  p) 2 ( I I I 112 -I I I 112 -I I I 112 (M -π))ê + ( I I I 122 -I I I 122 )ê ⊥ + 2(ê • p)(ê ⊥ • p) I I I 122 -I I I 122 -I I I 111 + I I I 111 ê + I I I 222 -I I I 222 -I I I 222 (M -π) -I I I 112 + I I I 112 + I I I 112 (M -π)

  ) I I 12 = I I 12 -I I 12 -I I 12 (M -π) 112 = I I I 112 -I I I 112 -II I 112 (M -π) = I V 112 -I V 112 -I V 112 (M -π) + e 2 )C 0,1 k + 60eC 0,2 k + 5(4 + 3e 2 )C 0,3 k + 15eC 0,4 k + 3e 2 C 0,5 k sin k M, (78d)I V 222 = I V 222 -I V 222 -I V222 (M -π) II = I I 11 -I I 11 -I I 22 + I I 22 IV = I 1 -I 1 + I I I 111 -I I I 111 -I I I 122 + I I I 122 -5( I V 122 -I V 122 ) e 2 )S 0,1 k + 36eS 0,2 k + (28 + 9e 2 )S 0,3 k + 18eS 0,4 k + 3e 2 S 0,5 k cos k M .

The open-source and extensively validated tool, Semi-analytic Tool for End of Life Analysis (STELA), developed by CNES, is also based on the equinoctial parameters (Le Fèvre et al.

2014).

An expression ab, formed by the juxtaposition of two vectors in R

constitutes what is called a dyad. For notational convenience, we denote the cross-product dyadic as a so that the tilde operator turns a vector into a skew-symmetric dyadic, allowing to write the cross product as a linear vector function. Equivalent expressions of the cross product are:a × b = a • b = a • b =b • a = -b • a.

Retaining the long-periodic perturbations in the transformation from osculating to mean elements is technically a misuse of Brouwer's theory[START_REF] Breakwell | On error bounds and initialization in satellite orbit theories[END_REF]); yet, it is, perhaps inadvertently, done in the literature[START_REF] Schaub | J 2 invariant relative orbits for spacecraft formations[END_REF] and we show herein that it does not have much bearing on the solutions outside the critical inclination.

 4 We have verified that the expressions for a and i do indeed agree, but have not done so for the remaining elements, which is a more painstaking and arduous task.

Truncating to order 12 in the Hansen coefficients was found to yield accurate results, in general, with higher orders being needed with increasing eccentricity.

Acknowledgements A deep debt of gratitude is extended to Claudio Bombardelli for help in clarifying a

Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.