Efficient formulation of the periodic corrections in Brouwer's gravity solution

The periodic terms of Brouwer's gravity solution are reconstructed in a nonsingular set of variables which are derived from the well-known polar-nodal variables. This change does not affect the essence of the solution, which still keeps all the benefits of the action-angle variables approach, and yields two major improvements. Namely, the periodic corrections of Brouwer's solution are now valid for any eccentricity below one and any inclination except the critical inclination, and, besides, are significantly simpler than the nonsingular corrections in Lydanne's reformulation of Brouwer's theory.


Introduction
Usual orbit propagators rely on Brouwer's (1959) analytical solution to the artificial satellite theory (AST). Brouwer limited to the main terms of the gravity potential and proceeded by averaging the mean anomaly and the argument of the perigee to decompose the solution into secular, long-period and short-period terms, and obtained a perturbation solution that takes into account the secular effects of the satellite problem up to the second order of J 2 -the coefficient of the zonal harmonic of degree 2 in the geopotential expansion-whereas the periodic effects are achieved up to the first order of J 2 .
In his pioneering approach, instead of resorting to the usual expansion of the gravity potential as a Fourier series in the mean anomaly with coefficients that are truncated expansions in the eccentricity, Brouwer is able to find a closedform solution by averaging implicit functions of the mean anomaly, in this way extending the validity of the theory to the higher-eccentricity orbits and making it more efficient due to the resulting shorter series. However, in spite of being a significant achievement, Brouwer's approach suffers from several drawbacks.
First of all, because Brouwer's solution is computed in Delaunay canonical variables, which are singular for circular and equatorial orbits, it introduces errors of the first order in the computation of periodic effects when the eccentricity approaches to zero. This drawback is commonly fixed by reformulating Brouwer's solution in non-singular variables as, for instance, Poincaré's canonical variables (Lyddane, 1963) -in which case the modified solution is sometimes renamed as Brouwer-Lyddane theory.
In second place, as a result of the double averaging both the Delaunay action and the polar component of the angular momentum are converted into formal integrals of the satellite problem, a fact that restricts the validity of the solution to a limited class of orbits. Indeed, after the first averaging over the mean anomaly the semi-major axis remains constant except for short-period effects, meaning that Brouwer's solution does not apply to orbits whose semi-major axis is affected of long-period or secular effects. As far as the gravity model assumes an axisymmetric earth semi-major axis resonances are not of concern in Brouwer's theory. On the other hand, the second averaging over the argument of the perigee 1 makes that both the eccentricity and the inclination remain constant except for periodic effects, thus excluding resonant-inclination orbits from the range of validity of the solution. In particular the critical inclination of 63.4 deg, which corresponds to the 1:1-resonance (Lara, 2014a), 2 remains out of the scope of Brouwer's solution. Since resonances are part of the dynamics this "issue" is not expected to be "fixed", but some tricks are customarily introduced in the code in order to avoid the problem of small divisors in software implementations of Brouwer's or other perturbation solutions by double averaging (Coffey et al., 1996;Hoots et al., 2004). Alternatively, analytical solutions based on intermediary orbits of AST do 1 Commonly, it is said that the averaging over the argument of the perigee is valid only for orbits with rotating periapsis, in this way missing the important class of frozen-perigee orbits (Cutting et al., 1978;Coffey et al., 1994;Lara, 1999Lara, , 2008. 2 Note that the critical inclination resonance and the so-called Lidov-Kozai resonance caused by third-body perturbations (Lidov, 1962;Kozai, 1962b) are of a similar nature. Still, the small value of the J 2 perturbation typical of AST, on one hand, and the double bifurcation phenomenon, on the other, make the characteristics of AST solutions radically different from the case of third-body perturbations. Indeed, while the eccentricity of librating-periapsis orbits may undergo important secular (or very-long period) variations with concomitant changes in inclination due to third-body perturbations, variations in eccentricity are of the order of J 2 in AST, and hence the libratingperiapsis orbits remain with almost constant "critical" inclination.
not suffer the small divisors problem in the critical inclination, but they are generally constrained to the contribution of first order effects of J 2 , although they may include some of the second-order effects of J 2 in specific degenerate cases. For a review on the topic of main problem intermediaries, see (Deprit, 1981); see also (Lara, 2014b).
Finally, the evaluation of the long Fourier series that provide the short-period corrections when expressed in Delaunay variables requires a non-negligible computational effort, yet notable savings are achieved with the efficient factorings found by Brouwer (1959). This fact motivated attempts in improving the efficiency of the Brouwer-Lyddane theory by reformulating it in new sets of nonsingular variables (Hoots, 1981), or by introducing different simplifications in software implementations of the Brouwer-Lyddane solution in order to speed computations, as is the case of SGP4 (Hoots and Roehrich, 1980;Vallado et al., 2006). The efficient evaluation of Fourier series is even of more concern when second-order periodic corrections are taken into account, a problem to which Kozai (1962a) found partial remedy by deriving the short-period corrections for the radial distance and the argument of the latitude rather than the mean anomaly and the argument of the perigee. Since Kozai found the relations among the different variables by means of differential formulas, he was compelled to compute first the short-period corrections in Delaunay variables in the usual way. However, as pointed out by Izsak (1963), the short-period corrections for the radial distance and the argument of the latitude can be derived directly by the simple expedient of furnishing the determining function of the canonical transformation of the averaging in polar-nodal variables, in this way greatly simplifying their computation.
The simplifications provided by polar-nodal variables in the computation of short-period corrections are easily justified looking at the main problem disturbing function, where the Legendre polynomial of degree 2 comprises the square of the cosine of the argument of the latitude (the true anomaly plus the argument of the perigee), on the one hand, while the cube of the radial distance divides the whole disturbing term, on the other. Since the inverse of the radial distance has a summand that is proportional to the cosine of the true anomaly, it happens that the disturbing function of the main problem will comprise powers from 2 to 5 of the cosine of the true anomaly, which is reminded to be an implicit function of the mean anomaly. Thus, the disturbing function of the main problem involves sine and cosine functions of different multiples m of the true anomaly with m ranging from 1 to 5, and this structure is preserved by the short-period transformation equations of the mean anomaly averaging. On the contrary, the disturbing function is dramatically simplified when using polar-nodal canonical variables instead of Delaunay variables. Now, the radial distance is a variable per se, and hence the number of trigonometric terms is limited to just one: the cosine of two times the argument of the latitude.
The advantages provided by polar-nodal variables motivated efforts by Cid and Lahulla (1969) in developing a polar-nodal variables alternative to Brouwer's solution with improved computational performance. The approach of Cid an Lahulla was succesful in what respects to first-order terms, where a single averaging over the argument of the latitude is enough for the reduction of the main problem to an integrable Hamiltonian. However, while the transformation equations of the averaging over the argument of the latitude are realized in polar-nodal variables by a simple and compact set of formulas of straightforward evaluation, the analytical solution depends on elliptic integrals whose evaluation is much more cumbersome than the evaluation of Brouwer's secular terms. Moreover, Cid an Lahulla did not find the same success when trying to extend their approach to higher orders, a case in which the transformation equations of the averaging are adversely affected of mixed secular-periodic terms (Cid and Lahulla, 1971). Remark that the appearance of mixed terms does not invalidate the procedure, but it does constrain the validity of the analytical solution to short time intervals.
Because polar-nodal variables are not action-angle variables of the Kepler problem (the unperturbed case of the satellite problem), the appearance of mixed secular-periodic terms in Cid and Lahulla's second-order solution is not unexpected. Hence, Kozai's (1962a) approach of using Delaunay variables in the computation of the secular terms while keeping the polar-nodal variables just for computing the short-period transformation equations seems to be the correct way. Still, Kozai computes the long-period corrections in Delaunay variables and, therefore, finds the same problems as Brouwer for dealing with low eccentricities and inclinations. Kozai's suggestion of replacing the mean anomaly, the eccentricity and the argument of the perigee by the mean argument of latitude and the projections of the eccentricity vector in the nodal frame avoids the singularities in the case of small eccentricities (Kozai, 1961), but cannot deal efficiently with the problem of small inclinations, which does not cause trouble in a perturbation solution of the main problem but it definitely does as soon as the effects of odd zonal harmonics are incorporated into the disturbing function.
To amend these issues, the suggestion of Izsak (1963) of taking advantage of the invariance of Poisson brackets with respect to canonical transformations in order to compute the periodic corrections directly in polar-nodal variables is carried out here to its full extent. Namely, while the averaging is performed in Delaunay variables, both generating functions of the short-and long-period elim-ination as well as the corresponding transformation equations are reformulated in polar-nodal variables, in this way avoiding the problem of low eccentricities, on one hand, and providing a compact formulation of the periodic corrections, on the other, which is of straightforward evaluation. This departure from Lyddane's (1963) approach of computing the satellite ephemeris through some set of nonsingular variables -not necessarily constrained to the case of canonical variables (Deprit and Rom, 1970)-whose flow is written as a function of the Delaunay secular terms, does not prevent completely the appearance of small divisors, which still may remain in the case of low inclinations. But this case is trivially reformulated and thoroughly simplified by using a new set of non-singular variables which, not surprisingly, are quite similar to the alternate set of variables used by Hoots (1981) in his attempts to accelerate evaluation of the Brouwer-Lyddane solution.
As another option, one might be tempted to try the alternative approach of computing the periodic corrections in the canonical set of Cartesian variables, which may well completely avoid the problem of small divisors. However, this non-singular transformation can only be achieved at the cost of a considerable increase in the size of the transformation equations, who still would need to be carefully inspected in order to trace terms affected of non-essential singularities related to equatorial orbits. These remaining critical terms should then be recast in such a form that the appearance of (virtual) small divisors could be fully prevented.

Short-period elimination
Since this research deals only with perturbations of gravitational origing, the problem of disturbed Keplerian motion can take benefit from Hamiltonian formulation. Thus, the motion of a massless particle in the gravitational field of the earth is derived from the Hamiltonian where H 0 represents the integrable Keplerian Hamiltonian and D is the disturbing function, which comprises the non-centralities of the geopotential. From the usual solution of Laplace's equation in spherical coordinates, the forces model is further limited to the zonal harmonics case, in which the disturbing function is written where µ is the earth's gravitational parameter, α is the earths' equatorial radius, r is the radial distance from the earth's center of mass, ϕ is latitude, P m,0 are Legendre polynomials of degree m, and C m,0 = −J m are corresponding zonal harmonic coefficients.
The problem of small inclinations in Brouwer's solution is related to the effects of odd zonal harmonics, so to illustrate this case it is enough to consider the impact of C 3,0 , in addition to the main problem, and hence the zonal gravitational potential in Eq. (2) is further truncated to the degree m = 3. Besides, because of the different orders of the harmonic coefficients, where J 3 = O(J 2 ) 2 , it is found convenient to make the Hamiltonian perturbative arrangement in which where the relation sin ϕ = sin I sin( f + ω) has been used, with I the orbit inclination, ω the argument of the perigee, and f the true anomaly, s and c are abbreviations for the sine and cosine of the inclination, respectively, a is the semi-major axis, and, in consequence with the Hamiltonian formulation, all the symbols, that is to say: a, r, ω, f , and I, are assumed to be functions of some set of canonical variables. In particular, since Brouwer's perturbation solution is computed by averaging, it relies in the action-angle variables of the Kepler problem, the so-called Delaunay variables. Namely, the mean anomaly ℓ and its conjugate momentum L = √ µ a (the Delaunay action), the argument of the perigee g = ω and its conjugate momentum G = L √ 1 − e 2 (the total angular momentum), where e is the orbital eccentricity, and the argument of the node h and its conjugate momentum H = G cos I (the polar component of the angular momentum). By using this canonical set it is simple to see that h is cyclic in Eq. (3) and, therefore, H is an integral.
The Hamiltonian reduction of Eq. (3) by perturbation methods is thoroughly documented in the literature, and hence results are provided without giving details in the method. In particular, the computations carried out are based on the implementation of the Lie transforms method known as Deprit's triangle algorithm, which is nowadays considered standard for Hamiltonian perturbations. Readers interested in the Lie transforms method can find all the required details in the original papers of Hori (1966) and Deprit (1969), as well as in modern celestial mechanics textbooks like (Meyer and Hall, 1992;Boccaletti and Pucacco, 2002), or other specialized books as (Ferraz-Mello, 2007).
At the first order of the perturbation approach, Deprit's triangle gives where {P, Q} notes the Poisson bracket of two functions P and Q of the canonical variables, which in this case are the Delaunay variables. Brouwer chooses the new Hamiltonian term H 0,1 as the averaging of H 1,0 over the mean anomaly 3 where, η is the eccentricity function and, for the sake of abbreviating notation, the function ǫ 2 ≡ ǫ 2 (G; µ) has been introduced, which is given by where p = G 2 /µ is the semilatus rectum. The corresponding term U 1 of the generating function is computed from Eq. (7) by quadrature where n = µ 2 /L 3 is the mean motion. In view of the differential relation a 2 η dℓ = r 2 d f , Eq. (9) can be integrated in closed form of the eccentricity to give +3es 2 sin( f + 2g) + 3s 2 sin(2 f + 2g) + es 2 sin(3 f + 2g) , (10) (cf. Eq. (15) of Brouwer, 1959, keeping where, here, ρ ∈ (ℓ, g, h, L, G, H) and W 1 = U 1 . Corresponding transformation equations in Delaunay variables can be expressed as Fourier series which involve sine and cosine functions of 10 different arguments of the form β = k f + 2mg with k = 0, . . . , 5 and m = −1, 0, 1 (cf. the first order terms in Eqs. (3.12) and (3.13) of Kozai, 1962a, for instance). However, important simplifications can be achieved by using the function instead of wholy expanding the transformation equations as Fourier series. In this way, the number of trigonometric arguments is reduced to just four: f , f + 2g, 2 f + 2g, and 3 f + 2g (cf. Eqs. (20) and (21) of Brouwer, 1959).
Alternatively, as pointed out by Izsak (1963), the generating function U 1 can be expressed in polar-nodal variables (r, θ, ν, R, Θ, N), which stand for the radial distance, the argument of latitude, the argument of the node, the radial velocity, the total angular momentum, and the polar component of the angular momentum. Rewriting Eq. (10) in polar-nodal variables as V 1 ≡ U 1 (r, θ, ν, R, Θ, N) is straightforward, leading to 4 where the functions are the projections of the eccentricity vector in the orbital frame when written in polar-nodal variables, which are trivially derived from Eq. (12) and its time derivative Then, the first order transformation equations are obtained in polar-nodal variables again from Eq. (11), where now ρ ∈ (r, θ, ν, R, Θ, N) and W 1 = V 1 . In this case, the equation of the center is φ ≡ φ(r, R, Θ), and, in particular, Hence, it is easily obtained which must be evaluated in prime variables for direct corrections ∆ρ = ρ − ρ ′ and in original variables from the inverse corrections ∆ρ ′ = ρ ′ − ρ. Remarkably, now the evaluation of the corrections only requires dealing with sine and cosine functions of the single argument 2θ. Note that the evaluation of the equation of the center φ = f − ℓ = f − u + e sin u is required in Eqs. (17) and (18), where u = 2 arctan 1 − e 1 + e tan f 2 , e = √ κ 2 + σ 2 , and f is unambiguously computed from cos f = κ/e and sin f = σ/e.
In the new notation, the first order of Deprit's triangle in Eq. (7) is rewritten Because K 1,0 in Eq. (24) does not depend on g, the new first-order Hamiltonian term is chosen K 0,1 = K 1,0 , and hence {K 0,0 , X 1 } = 0 from Eq. (25). However, this does not mean to make null the first order term X 1 of the long-period generating function. Quite on the contrary, since the generating function of the long-period averaging does not depend on ℓ then {K 0,0 , X 1 } necessarily vanishes in Eq. (25). Therefore, the term X 1 can only be determined at the next order of the perturbation algorithm. The Poisson bracket {K 0,0 , X 2 } vanishes likewise, and the second order of Deprit's triangle in Eq. (22) is simplified in this case to where the term K 0,2 is chosen as the average of K 2,0 over the argument of the perigee. That is, It follows the computation of X 1 from Eq. (26) by a quadrature: which is trivially solved to give where, for conciseness, the notation has been introduced. Next, the long-period generating function Y 1 = −ǫ 2 Θ s 2 14 − 15s 2 8 4 − 5s 2 κ 2 − σ 2 sin 2θ − 2κ σ cos 2θ +ǫ 3 Θ s (κ cos θ + σ sin θ).

The case of low inclinations
Due to the contribution of the odd zonal harmonic C 3,0 , it happens that δθ and δν are singular for equatorial orbits. However, as the simple inspection of Eqs. (32) and (33) may suggest, the trouble in the case of low inclinations is easily remedied by computing the long-period corrections to the non-singular elements Indeed, terms of the order of s 2 and higher only produce higher order effects in the lower inclination orbits, and hence can be neglected. Therefore, the corrections are straightforwardly derived from Eqs. (31)-(36). The transformation from non-singular to Cartesian variables is obtained by means of the usual rotations applied to the projections of the position and velocity vectors in the orbital frame. Thus, where R 1 , R 3 , are the usual rotation matrices After replacing ν = ψ − θ and sin θ = ξ/s, cos θ = χ/s, in Eq. (45), the transformation from non-singular to Cartesian variables can be obtained from the sequence x = r (t cos ψ + q sin ψ), (46) y = r (t sin ψ − q cos ψ), X = R (t cos ψ + q sin ψ) − Θ r (q cos ψ + τ sin ψ), Finally, in spite of a state (x, y, 0, X, Y, 0) corresponding to an exactly (instantaneous) equatorial orbit would be rarely obtained when working in real arithmetic, it deserves to mention that in the space (original or double prime) that it might happen ξ = χ = 0 and it does not make sense to speak of the node or the argument of latitude. However, periodic corrections still exist for ξ and χ. Indeed, while short-period corrections ∆ξ and ∆χ vanish for equatorial orbits, as derived from Eqs (61)-(62), corresponding long-period corrections do not, and Eqs. (39) and (40) result in δξ = −ǫ 3 σ, δχ = ǫ 3 κ.

Conclusions
Soon after Brouwer's solution was announced, the reformulation of the shortperiod corrections in polar-nodal variables was suggested as a way of simplifying their evaluation. This reformulation is extended here to the case of the long-period corrections to show that, as odd as it may seem to introduce short-period terms in the computation of long-period corrections, this artifact prevents the usual deterioration of the corrections in the case of low-eccentricity orbits. The new formalism yields significantly less computational effort than Lyddane's non-singular variables approach, yet the case of low-inclination orbits must be treated separately. However, the elementary inspection of the long-period corrections in polar-nodal variables reveals a simple set of non-singular elements that may be used for dealing properly with that case. This lack of universality in the new formalism is easily overcome by introducing a flag in the code, and is largely compensated by the notable simplifications that the periodic corrections accept in the case of almost-equatorial orbits.