BBrouwer’s satellite solution redux *Martin Lara †‡ November 3, 2020
Abstract
Brouwer’s solution to the artificial satellite problem is revisited to showthat the complete Hamiltonian reduction is rather achieved in the plain Poin-caré’s style, through a single canonical transformation, than using a sequenceof partial reductions based on von Zeipel’s alternative for dealing with per-turbed degenerate Hamiltonian systems.
Brouwer’s analytical solution to the artificial satellite problem [7] based on vonZeipel’s partial reduction method for dealing with perturbed degenerate Hamilto-nians [42] fiercely resists obsolescence sixty years after publication. Indeed, inspite of the spectacular increase of the computational power, widespread softwarepackages for approximate ephemeris prediction still rely on Brouwer’s seminal re-sults [18, 41]. Furthermore, the success of Brouwer’s closed-form solution amongpractitioners as well as the reputation gained among theorists by Brouwer’s step-wise normalization, make that these days some authors designate the method as“the von Zeipel-Brouwer theory” [15].Since then, the merits of Brouwer’s decomposition of the solution of perturbedKeplerian motion into secular, long-, and short-period effects, seem not to havebeen questioned. Moreover, after the invention of Hamiltonian simplification meth-ods [10], it has been suggested that carrying out additional decompositions, thusincreasing the number of canonical transformations, could be the proper way tosuccess in the search for separable perturbation Hamiltonians of celestial mechan-ics problems [11]. Conversely, it has been recently pointed out that the use of * Submitted to Celestial Mechanics and Dynamical Astronomy † GRUCACI, University of La Rioja, and Space Dynamics Group – UPM ‡ [email protected] a r X i v : . [ m a t h . D S ] N ov amiltonian simplification procedures could be merely optional in the construc-tion of higher-order analytical approximations to the satellite problem [30]. Then,it emerges the question of which is the real value of splitting a normalization pro-cedure, either partial or complete, into several different stages, a topic that maywell deserve additional study.We walk a step in that direction and, in the classical style of Poincaré’s “newmethod” [37], undertake the construction of Brouwer’s second-order completelyreduced Hamiltonian of the main problem in the artificial satellite theory (the J problem) by means of a single canonical transformation. The difficulties stemmingfrom the degeneracy of the Kepler Hamiltonian, who has two null frequencies, areeasily overcome with the addition of suitable integration constants to the generatingfunction of the transformation that yields the complete Hamiltonian reduction.The use of arbitrary functions in the construction of perturbation solutions isnot new at all [35]. In fact, it can be traced back to Poincaré’s efforts in approach-ing degenerate perturbation problems [37, Chap. XI]. They also play a fundamentalrole in a reverse partial normalization process in which the angular momentum isnormalized in first place [3, 28, 38, 30]. On the other hand, in spite of average per-turbation Hamiltonians do not exist in general [14], the use of arbitrary constants toguarantee the purely periodic nature of the generating function became customaryin attempts to bring the mean elements dynamics as close as possible to the trueaverage dynamics [33, 39, 27].To the first order, the construction of Brouwer’s closed-form solution by meansa single transformation amounts to the sum of the two transformations computedby Brouwer for the short- and long-period elimination. This is due to the lineariza-tion provided by the first order truncation of a perturbation theory. In view of nodifferences arise between Brouwer’s and the current approach when the periodiccorrections are constrained to first order effects, we feel compelled to supplementBrouwer’s analytical solution with second-order periodic corrections, yet limitedto the J contributions. We compare our results with corresponding correctionsobtained in the traditional way, in which the normalization is split into the elimina-tion of the parallax, the elimination of the perigee, and the Delaunay normalization[3, 26]. At the second order the single transformation is no longer the addition ofthe different canonical transformations. As expected, the (single) periodic correc-tions are now much more involved than those corresponding to each of the partialreductions or simplifications, and are also more intricate than the composition ofall of them. However, the length of the series defining the solution is only a partof the whole picture, and we found clear computational advantages in the evalua-tion of the single transformation. The improvements stem from the fact that manyinclination polynomials pertaining to the periodic corrections admit factorization.Because common factors repeat many times throughout the corrections, the com-2iler is able to perform a higher optimization of the code in the case of the singletransformation than in the case in which the transformation is split in differentstages.In the process of computing the second-order corrections, we will recognizehow artificial the controversy created about the integration of the equation of thecenter was. Indeed, difficulties confronted by researchers involved in the automa-tization of celestial mechanics computations were, in fact, derived from their ownprogramming strategies [12, 21]. On the contrary, the trouble had been easilysidestepped by celestial mechanicians relying on traditional hand computations[23, 1]. We will show that standard integration by parts reduces the equation ofthe center issue to the well known integration of cosine functions in elliptic mo-tion [22, 40]. We hasten to say that the controversy was in no way futile sinceit provoked the appearance of Hamiltonian simplification methods and led to thedevelopment of sophisticated computational strategies [16]. In order to fully determine the second-order term of the generating functionof Brouwer’s theory, the third-order term of the completely reduced Hamiltonianneeds to be previously specified. The use of higher-order secular terms should im-prove further the long-term performance of Brouwer’s solution. However, to beeffective in the propagation of an initial state vector, the initialization constantsof the analytical solution, and, in particular, the secular frequencies, must be com-puted within comparable accuracy to that of the secular terms. Rather than carryingout the long and tedious computations required in the determination of the third-order generating function, we take the clever shortcut proposed by Breakwell andVagners [6]. That is, we limit the computation of third-order corrections to thecase of the secular mean motion, which, besides, is directly obtained from the sec-ular Hamiltonian. With this effortless procedure the addition of third-order secularterms clearly improves the performance of Brouwer’s solution.
Constraining the dynamical model of the artificial satellite problem to the J per-turbation (main problem), Brouwer’s gravitational Hamiltonian takes the form [7] H = − µ a + µr R ⊕ r C , (cid:20) − s + 32 s cos(2 f + 2 ω ) (cid:21) , (1)where the Earth’s gravitational field is materialized by the physical constants µ ,the gravitational parameter, R ⊕ , the equatorial radius, and C , = − J , the non- A brief review of the history of Hamiltonian simplification methods can be found in the Intro-duction of [29]. The symbols a , r , f and ω , stand for semi-major axis, radius from the Earth’s center of mass, true anomaly, and argument ofthe perigee, respectively, whereas s ≡ sin I abbreviates the sine of the inclination I . Since we are dealing with Hamiltonian mechanics, these symbols must be un-derstood as functions of some set of canonical variables. In particular we assume,with Brouwer, that the Hamiltonian is written in terms of the Delaunay coordinates (cid:96) , the mean anomaly, g = ω , and h , the longitude of the ascending node, and theirconjugate momenta L = √ µa , G = L (1 − e ) / , with e denoting the eccentricity,and H = G cos I , standing for the Delaunay action, the total angular momentum,and the projection of the angular momentum vector along the Earth’s rotation axis,respectively. That H is an integral of Eq. (1) becomes evident from the cyclic char-acter of h . Besides, the Hamiltonian itself is constant because the time does notappear explicitly on it.The small value of the Earth’s J coefficient identifies Eq. (1) like a case of per-turbed Keplerian motion, which, therefore, can be reduced to a separable Hamil-tonian by perturbation methods. This is achieved by finding a canonical transfor-mation T : ( (cid:96), g, h, L, G, H, (cid:15) ) (cid:55)→ ( (cid:96) (cid:48) , g (cid:48) , h (cid:48) , L (cid:48) , G (cid:48) , H (cid:48) ) , from osculating to meanvariables, depending on the small parameter (cid:15) ∼ J , such that the transformedHamiltonian in mean (prime) variables becomes a function of only the momenta,namely H ◦ T = K ( − , − , − , L (cid:48) , G (cid:48) , H (cid:48) ; (cid:15) ) . The transformation T , we learnedfrom Poincaré [37], is derived from a determining function that is solved in theform of a Taylor series up to some truncation order of the small parameter (cid:15) .Brouwer, for his part, after introducing the method of solution, seems to refuseapproching the direct computation of the transformation T since the beginning, bysimply declaring that“[. . . ] it is more convenient to choose a determining function in sucha manner that the mean anomaly is not present in the transformedHamiltonian while the argument of the perigee is permitted to appear.”Next, after invoking von Zeipel, he proceeded stepwise by partial reduction, firstcomputing a canonical transformation that only removes the short-period termsfrom the Hamiltonian, and then carrying out a second canonical transformationthat removes the long-period terms. In this way Brouwer outstandingly achievesthe complete Hamiltonian reduction in closed form.Conversely, we ignore the presumed convenience of Brouwer’s procedure andapproach the perturbation problem searching for a single determining function inthe original style of Poincaré, yet we better rely on the equivalent but more func-tional method of Lie transforms [19, 9, 13]. Thus, we write Eq. (1) in the usual Note that k = − C , R ⊕ in Brouwer’s notation. H = (cid:80) m ≥ ( (cid:15) m /m !) H m, , with H , = − µ a , H , = − µr R ⊕ r (cid:20) − s + 32 s cos(2 f + 2 ω ) (cid:21) , H m, = 0 for m ≥ , and (cid:15) ≡ J = − C , . Recall, that all the symbols arefunctions of the Delaunay canonical variables, in which the Lie operator L = { ; H , } , where the curly brackets stand for the Poisson brackets operator, takesthe simple form L = n∂/∂(cid:96) , where n = µ /L is the mean motion. This al-lows us to compute the determining function W = (cid:80) m ≥ ( (cid:15) m /m !) W m +1 fromthe sequence given by the homological equation W m = 1 n (cid:90) ( (cid:101) H ,m − H ,m )d (cid:96) + C m . (2)At each step m , terms (cid:101) H ,m in Eq. (2) are known, coming either from the origi-nal Hamiltonian or stemming from intermediate computations at previous orders.Terms H ,m are selected in such a way that they cancel those terms of (cid:101) H ,m per-taining to the kernel of the Lie operator. Finally, the integration “constants” C m —arbitrary functions of the Delaunay variables fulfilling the condition ∂ C m /∂(cid:96) =0 — will be chosen like such trigonometric functions of g that they prevent the ap-pearance of purely long-period terms at the next order of the perturbation approach,in this way making feasible the complete normalization at once. The method isstandard these days, and the required details can be found in textbooks as, for in-stance, [34, 4].In preparation of the solution, the equivalence r j = 1 r (cid:18) e cos fp (cid:19) j − , (3)where p = aη is the orbit parameter and η = (1 − e ) / is the so-called eccen-tricity function, is applied to the instances j > in (cid:101) H , = H , , which is thenwritten in the convenient form (cid:101) H , = H , R ⊕ r η (cid:88) i =0 B i ( s ) i +1 (cid:88) j = i (2 − j (cid:63) ) i e | j − i | cos( jf + 2 ig ) , (4)where B = 1 − s , B = s , and we abbreviate j (cid:63) ≡ j mod 2 . On accountof j ≥ k in Eq. (4), we immediately verify that (cid:101) H , is not affected of purelylong-period terms. Then, the complete reduction is achieved at the first order bychoosing the new Hamiltonian term H , like the average of H , over the meananomaly. 5he average is obtained in closed form with the help of the Keplerian differen-tial relation between the true and mean anomalies ηa d (cid:96) = r d f . It is equivalentto removing all the terms with j > from Eq. (4) after multiplied by the factor r / ( a η ) . We trivially obtain the usual result H , = (cid:104) (cid:101) H , (cid:105) (cid:96) ≡ H , R ⊕ p η (cid:18) − s (cid:19) , which is, of course, the same expression obtained by Brouwer. Then, Eq. (2) isrearranged in the form W = 1 n (cid:34) H , φ + (cid:90) (cid:32) (cid:101) H , r a η − H , (cid:33) d f (cid:35) + C , (5)where φ = f − (cid:96) denotes the equation of the center, and the integrand in Eq. (5)only embraces periodic functions of f . We obtain W = − G R ⊕ p (cid:20) B φ + (cid:88) i =0 B i i +1 (cid:88) j =max( i, (2 − j (cid:63) ) i j e | j − i | sin( jf + 2 ig ) (cid:21) + C , (6)where the first term of the right hand member is the same as Brouwer’s first or-der determining function of the short-period elimination, and C is an integrationconstant that is left undetermined by the time being.On account of H , ≡ , the known terms at the second order of the Lie trans-forms approach are (cid:101) H , = {H , ; W } + {H , ; W } , from which the terms per-taining to the kernel of the Lie operator must be cancelled by the adequate selectionof H , . The usual choice is H , = (cid:104) (cid:101) H , (cid:105) (cid:96) , yet additional terms could be left inthe new Hamiltonian in particular cases [10, 29]. However, this process wouldleave purely long-period terms in the new Hamiltonian in addition to the secularterms, both certainly pertaining to the kernel of the Lie operator. Since this isagainst the total normalization criterion, purely long-period terms should vanishidentically in (cid:101) H , , a requirement that is achieved with the proper selection of C ,whose partial derivatives with respect to g , G , and L , appear formally in (cid:101) H , .We attack the computation of the second-order of the perturbation theory byparts. To that effect, we make (cid:101) H , = (cid:101) H (cid:48) , + (cid:101) H ∗ , with (cid:101) H (cid:48) , = {H , + H , ; V } and (cid:101) H ∗ , = {H , + H , ; C } . Straightforward evaluation of the Poisson brackets,followed by the use of Eq. (3) and standard trigonometric reduction, yields (cid:101) H (cid:48) , = H , α p a r η η (cid:88) i =0 s i i +4 (cid:88) j =( − i i −| i − j | (cid:88) k =0 B i,j,k η k e | j − i | cos( jf + 2 ig ) H , α p η (cid:20) η (3 s − + 3(5 s − s (cid:88) j =1 − j (cid:63) j e j (cid:63) cos( jf + 2 g ) (cid:21) + H , α p (cid:16) s − (cid:17) s p r φη (cid:88) j =1 (2 − j (cid:63) ) e j (cid:63) sin( jf + 2 g ) , (7)where the needed coefficients B i,j,k ( s ) are listed in Table 1.Table 1: Inclination polynomials B i,j,k in Eq. (7). B , , = B , , B , , = 3 (cid:0) s − (cid:1) B , , = B , , B , , = 10 (cid:0) s − s + 8 (cid:1) B , , = 12 (cid:0) s − (cid:1) B , , = 5 B , , = B , , B , , = B , , B , , = B , , B , , = 2 (cid:0) s + 8 s − (cid:1) B , , = B , , B , , = 6 B , , = 2 (cid:0) s − s + 60 (cid:1) B , , = − (cid:0) s − (cid:1) B , , = B , , = 0 B , , = 4 (cid:0) s − s + 24 (cid:1) B , , = B , , B , , = B , , B , , = B , , B , , = − (cid:0) s − (cid:1) B , , = − B , , = 2 (cid:0) s − s + 24 (cid:1) B , , = 150 − s B , , = B , , B , , = − (cid:0) s + 8 s − (cid:1) B , , = − (cid:0) s − (cid:1) B , , = − B , , = 2 (cid:0) s − (cid:1) B , , = 2 − s B , , = B , , B , − , = B , , B , , = − (cid:0) s − (cid:1) B , , = − B , , = B , , B , , = − (cid:0) s − (cid:1) B , , = B , , B , , = 4 (cid:0) s − (cid:1) B , , = − (cid:0) s − (cid:1) B , , = − We intentionally split (cid:101) H (cid:48) , into three different blocks. Namely, all the termson the first row of Eq. (7) are free from the equation of the center and factoredby a /r , hence being of trivial integration in the true anomaly. Terms of thesecond row are free from both φ and r ; the integration of terms of this type reducesto the well-known case of the integration of cosine functions in elliptic motion[7, 22, 40]. Terms on the third row are of the form ( p/r ) φ sin( mf + α ) , with m integer, and are readily integrated by parts. That is, on account of d(cos mf ) / d (cid:96) = − ( m/η )( p/r ) sin mf , mη (cid:90) p r φ sin mf d (cid:96) = − φ cos mf + sin mfm − (cid:90) cos mf d (cid:96), (8)in this way leading to the previous case of integration of cosine functions. Partic-ularization for definite integration follows from the fundamental theorem of calcu-lus. It is worth noting that if, on the contrary, we replace r by the conic equation in7he third row of Eq. (7), in order to arrange this row like the Fourier series H , α p
98 (5 s − s φ (cid:88) j = − q | j − | e | j − | sin( jf + 2 g ) , with q = 3 e +2 , q = ( e +4) , q = , and q = , then the equation of the cen-ter shows as an isolated function of the mean anomaly when the summation indextakes the value j = 0 . This arrangement brings no problem in the computation ofdefinite integrals, which can be carried out using the general rules for computing (cid:104) φ sin mf (cid:105) (cid:96) and (cid:104) φ cos mf (cid:105) (cid:96) provided in [32]. On the contrary, while indefiniteintegration is still possible, it requires the sophisticated use of special functions,which could make notably difficult to progress in the perturbation approach [36].On the other hand, the evaluation of the Poisson brackets involving the integra-tion constant C yields (cid:101) H ∗ , = H , R ⊕ p L ∂ C ∂g (cid:20) − s + a ηr (cid:88) i =0 2 i +3 (cid:88) j = − i j (cid:63) (cid:88) k =0 b i,j,k η k e j (cid:48) cos( jf + 2 ig ) (cid:21) −H , R ⊕ p ∂ C ∂G ηs a ηr (cid:88) j =1 [1 + ( j + 1) (cid:63) ] e j (cid:63) sin( jf + 2 g )+ H , R ⊕ p ∂ C ∂L η a ηr (cid:88) i =0 b i i +3 (cid:88) j = − i q i,j e j (cid:63) sin( jf + 2 ig ) , (9)where j (cid:48) = ( | j − i | − j (cid:63) , ( j + 1) (cid:63) ≡ ( j + 1) mod 2 , and the inclination poly-nomials b i,j,k ( s ) and the eccentricity polynomials q i,j ( e ) are provided in Table 2.Table 2: Non-null inclination b i,j,k and eccentricity polynomials q i,j in Eq. (9). b , , = 4 − s b , − , = s q , = e + 4 b , , = − (29 s − b , , = s − q , = 4 e b , , = (17 s − b , , = 1 − s q , = q , − = e b , , = 2 − s b , , = 5 s − q , = − q , b , , = − (3 s − b , , = s − q , = − e + 2) b , , = 1 − s q , = − q , b , , = 3 s q , = − e b , , = s q , = − e In the same way as we did in the first order, we choose H , = (cid:104) (cid:101) H , (cid:105) (cid:96) = (cid:104) (cid:101) H (cid:48) , (cid:105) (cid:96) + (cid:104) (cid:101) H ∗ , (cid:105) (cid:96) to guarantee that it cancels all the terms of (cid:101) H , perteining to thekernel of the Lie derivative. Firstly, we compute (cid:104) (cid:101) H (cid:48) , (cid:105) (cid:96) as follows. To average8he first row of Eq. (7) over the mean anomaly, it is first multiplied by the factor r / ( a η ) to carry out the integration in the true anomaly, and then those terms thatare free from f , which are those with j = 0 , are selected. The term free from f in the second row averages to itself while the remaining terms in this row areaveraged using the rule π (cid:90) π cos( mf + α ) d (cid:96) = (cid:18) − e η (cid:19) m (1 + mη ) cos α, (10)cf. [22]. Finally, the terms on the third row of Eq. (7) are averaged by parts withthe help of Eqs. (8) and (10). We finally obtain, (cid:104) (cid:101) H (cid:48) , (cid:105) (cid:96) = H , R ⊕ p η (cid:104) s − s + 8) + η (6 s − + η (5 s + 8 s − (cid:105) + H , R ⊕ p η (15 s − s e cos 2 g, (11)which is precisely Brouwer’s second-order Hamiltonian after the elimination ofshort-period terms. The average of Eq. (9) is readily obtained with analogous pro-cedures, to obtain (cid:104) (cid:101) H ∗ , (cid:105) (cid:96) = −H , R ⊕ p s −
4) 1
L ∂ C ∂g . (12)Visual inspection of Eqs. (11) and (12), immediately shows that if we complete thecomputation of the first-order term of the generating function in Eq. (6) choosing C = G R ⊕ p s − s − s e sin 2 g, (13)then Eq. (12) turns into the opposite of the term in the last row of Eq. (11), the onlyone that depends on g , thus mutually canceling out. Hence H , = H , R ⊕ p η (cid:104) s − s + 8) + η (6 s − + η (5 s + 8 s − (cid:105) , (14)which is the same as the second-order term of Brouwer’s Hamiltonian after theelimination of long-period terms. In this way we have achieved Brouwer’s totalHamiltonian reduction of the main problem at once, with a single canonical trans-formation.It is not a big surprise that C is the same integration constant obtained inAlfriend and Coffee’s elimination of the perigee [3, 28] or in the author’s reverse9ormalization of the angular momentum [30], since the motion in the orbital planeis decoupled from the rotation of that plane in each case.The computation of first-order periodic corrections is now straightforward fromthe simple evaluation of Poisson brackets, namely ξ − ξ (cid:48) = J ∆ ξ , where ∆ ξ ≡{ ξ, W } and ξ denotes either a canonical variable or some wanted function of thecanonical variables [9]. For instance, for the first-order periodic corrections to thesemi-major axis we obtain ∆ a = a R ⊕ p η (cid:88) i =0 B i ( s ) i (cid:88) j = − i A i,j ( η ) e | j − i | cos( jf + 2 ig ) , (15)where A , = 10 − η − η , A , = A , = A , = 15 − η , A , = A , = A , = 6 , A , = A , − = A , = 1 , and the coefficients B i are the same asthose in Eq. (4). Recall that Eq. (15) must be evaluated in mean (prime) variablesin the direct transformation from mean elements to osculating ones, and in orig-inal (unprimed) variables in the inverse transformation from osculating to meanvariables. The second-order term of the generating function is now computed making m = 2 in Eq. (2). Namely W = V + C , with V = 1 n (cid:90) ( (cid:101) H , − H , )d (cid:96). The needed integrals in the computation of V are either trivial, solved with thehelp of Eq. (8) for those terms involving the equation of the center, or using thedifferential relation between the mean and true anomalies for those other that arefree from φ but depend on trigonometric functions of f . Straightforward computa-tions yield V = G R ⊕ p φ (cid:20) − η (5 s + 8 s − − s − s + 8) − (15 s − e s cos 2 g + 12 s (5 s − (cid:88) j =1 − j (cid:63) j e j (cid:63) cos( jg + 2 g ) (cid:21) + G R ⊕ p (cid:88) i =0 j max (cid:88) j = j min (cid:88) k =0 β i,j,k ( s ) η k s i e j (cid:63) sin( jf + 2 ig )(5 s − − i (cid:63) (1 + η ) (cid:98) (3 − i ) (cid:99) , (16)10able 3: Non-null inclination polynomials β i,j,k in Eq. (16). , , : − s − s − s + 2400 s − , , : 12( − s + 16 s + 4) , , : − s − s − s + 8928 s − , , : 2(1855 s − s + 972) , , : 3( − s + 3030 s − s + 2368 s − , , : 2(1045 s − s + 540) , , : 3 s (975 s − s + 1728 s − , , : − s − s − , , : − β , , , , : − s − s − , , : − β , , , , : − β , , , , : 6(1925 s − s + 7452 s − s + 768) , , : − β , , , , : 6(125 s − s + 1660 s − s + 256) , , : − s − s − , , : − β , , , , : − s − s − , , : − β , , , , : − β , , , , : 2625 s − s + 7408 s − s + 512 , , : − s − s − , , : s (825 s − s + 1616 s − , , : − β , , , − , : − β , − , , , : 3(225 s − s + 208) , − , : − β , − , , , : − β , , , − , : 6(135 s − s + 100) , , : 60(50 s − s + 38) , − , : 6(7 s − s − , , : − s − s + 122) , , : − s − s + 364) , , : 8(75 s − s + 61) , , : − s − s + 656) , , : − s − s − , , : 48(5 s − , , : 12(5 s − s − , , : − s − s − , , : 3(5 s − s − , , : 12( − s + 240 s − , , : 3(5 s − s − , , : 12( − s + 240 s − , , : − β , , , , : β , , , , : − s − where j min = 2( i + 1) (cid:63) − , j max = 4 + i + (cid:98) ( i − (cid:99) , and the inclinationpolynomials β i,j,k are listed in Table 3.As before, the integration constant C will be determined by imposing to theknown terms of the next order (cid:101) H , = {H , + H , , W } + {H , + 2 H , , W } , where H , = H , + {H , , W } , the condition of being free from pure long-periodic terms. Again, the known terms are split into terms directly computableand those depending on the arbitrary function C . That is, (cid:101) H , = (cid:101) H (cid:48) , + (cid:101) H ∗ , ,where (cid:101) H (cid:48) , = {H , + H , , W } + {H , + 2 H , , V } , (cid:101) H ∗ , = {H , + 2 H , , C } . It follows the customary computation of H , so that it cancels the terms of (cid:101) H , H , = (cid:104) (cid:101) H , (cid:105) (cid:96) = (cid:104) (cid:101) H (cid:48) , (cid:105) (cid:96) + (cid:104) (cid:101) H ∗ , (cid:105) (cid:96) . After straightforward evaluation of the Poisson brackets, we obtain (cid:104) (cid:101) H (cid:48) , (cid:105) (cid:96) = H , R ⊕ p η (cid:88) i =0 4 − i + i (cid:63) (cid:88) k =0 β i,k ( s ) η k s i e i (5 s − − i (cid:63) (1 + η ) i (cid:63) cos 2 ig, (17)where i (cid:63) = i mod 2 and the inclination polynomials β i,k are given in Table 4.Analogously, (cid:104) (cid:101) H ∗ , (cid:105) (cid:96) = −H , R ⊕ p
92 (5 s −
4) 1
L ∂ C ∂g . (18)In this process, we only found integrals of the same type as we did at the secondorder, and hence there were no special difficulties in solving them, yet in this casewe needed to deal with notably longer series than in previous orders.Table 4: Inclination polynomials β i,k in Eq. (17). β , = − (cid:0) s − s + 158960 s − s + 45152 s − (cid:1) β , = − (cid:0) s − (cid:1) (cid:0) s − (cid:1) (cid:0) s − s + 8 (cid:1) β , = 2 (cid:0) s − s + 130852 s − s + 30176 s − (cid:1) β , = − (cid:0) s − (cid:1) (cid:0) s − (cid:1) (cid:0) s + 8 s − (cid:1) β , = s (cid:0) s − (cid:1) (cid:0) s − s + 590 s − (cid:1) β , = 525 s − s + 5632 s − β , = 5925 s − s + 14848 s − β , = (cid:0) − s (cid:1) (cid:0) s − s + 120 (cid:1) β , = (cid:0) s − (cid:1) (cid:0) s + 36 s − (cid:1) β , = (cid:0) s − (cid:1) (cid:0) s − (cid:1) Once more, the simple inspection of Eqs. (17) and (18) shows that if we nowchoose C = G R ⊕ p (cid:88) i =1 4 − i + i (cid:63) (cid:88) k =0 β i,k ( s ) η k s i e i (5 s − i +1 (1 + η ) i (cid:63) sin 2 ig i , (19)then Eq. (18) cancels the terms of Eq. (17) depending on the argument of theperigee out, to yield H , = H , R ⊕ p η (5 s − (cid:88) k =0 β ,k ( s ) η k , (20)12hich is completely reduced as desired.Beyond the first order, direct and inverse transformations are no longer oppo-site. At the second order the direct transformation is given by ξ = ξ (cid:48) + J ∆ ξ + J δ (cid:48) ξ , where δ (cid:48) ξ = { ∆ ξ, W } + { ξ, W } is evaluated in prime variables. Theinverse transformation is ξ (cid:48) = ξ − J ∆ ξ + J δξ , where δξ = { ∆ ξ, W } + { ξ, −W } , is evaluated in the original variables. For instance, replacing ξ by a weobtain the inverse second-order periodic correction to the semi-major axis δa = a R ⊕ p η (cid:20) η (5 s + 8 s −
8) + 48 η (15 − s e cos 2 g + (cid:88) i =0 6+2 i (cid:88) j = − i − i (cid:63) −| j − i | (cid:88) k =0 (3 s − i (cid:63) s i A i,j,k ( s ) η k e | j − i | cos( jf + 2 ig ) (cid:21) , (21)where the inclination coefficients A i,j,k are provided in Table 21.Table 5: Non-null inclination coefficients A i,j,k in Eq. (21). , − , : 9 , , : 1980 , , : − , , : 462 A , , , − , : 108 , , : − , , : 15120 , , : − A , , , , : 594 , , : 594 , , : 480 , , : 10 A , , , , : − , , : − , , : − , , : 210 A , , , , : 1980 , , : 108 , , : − s − s − , , : 24(71 s − s + 56) , , : − , , : 9 , , : 240 , , : − s − s + 8) , , : 4455 , − , : − , , : − , , : 792 A , , , , : − , − , : − , , : 8640 , , : − A , , , , : 135 , − , : − , , : 360 , , : 15 A , , , , : 7128 , − , : 72 , , : − , , : 120 A , , , , : − , − , : − , , : − s − s − , , : − A , , , , : 1080 , − , : 720 , , : − , , : 495 A , , , , : 8316 , − , : 24 , , : 3240 , , : − A , , , , : − , , : − , , : 144 , , : 6 A , , , , : 3780 , , : 3240 , , : − , , : 15 A , , , , : − , , : 144 , , : − , , : 220 A , , , , : 7128 , , : − , , : 720 , , : − A , , , , : − , , : − , , : 24 , , : − s − , , : 1080 , , : 8640 , , : − , , : 66 A , , , , : 4455 , , : 360 , , : 72 , , : − A , , , , : − , , : − , , : − , , : 12 A , , , , : 135 , , : − s − s − , , : − , , : 2(27 s − s + 8) Initialization of the secular constants and performancetests
Soon after Brouwer’s solution appeared in print, different reports pointed out anapparent contradiction between the accuracy expected from the series truncationorder and the comparatively large in-track errors obtained in a variety of testsagainst numerical integrations [5]. The issue, however, did not happen when fittingBrouwer’s solution to observational data. Hence, the apparent discrepancy waseasily identified with an inconsistency in the use of Brouwer’s theory. Indeed, toget the expected accuracy provided by the secular terms, the initialization of theconstants of Brouwer’s solution should be done with analogous accuracy. How-ever, Brouwer only provided the periodic corrections up to the first order of J ,and hence the direct initialization of the secular mean motion for given initial con-ditions yields analogous accuracy. The trouble is, of course, solved if the inverseperiodic corrections are computed up to the same order as the secular terms.On the other hand, since the trouble arises from an inaccurate computation ofthe secular mean motion, the theory can be patched by supplementing Brouwer’sfirst order corrections only with the inverse second-order correction to the semi-major axis, either using Eq. (21) or in the different much shorter formulation givenby [31]. Alternatively, the errors in the initialization procedure can be palliatedby fitting the secular frequencies to data obtained from a preliminary numericalintegration over several revolutions, or by a calibration of the secular mean motion n (cid:48) = µ /L (cid:48) from the energy equation [6].The latter approach is particularly appealing because it totally avoids the needof carrying out additional computations to those already carried out by Brouwer.Thus, for given initial conditions ( (cid:96) , g , h , L , G , H ) , the initial Hamiltonianin osculating elements evaluates to H ( (cid:96) , g , − , L , G , H ) = E . On the otherhand, after the complete Hamiltonian reduction E = − µ L (cid:48) + k (cid:88) m =1 J m m ! H ,m ( L (cid:48) , G (cid:48) , H ) + O ( J k +12 ) . (22)However, the constants L (cid:48) and G (cid:48) are computed from the osculating initial condi-tions through the inverse periodic corrections only up to O ( J k − ) . While this factdoes not compromise the accuracy of Eq. (22) in what respects to the terms H ,m ( m ≥ ) because they are multiplied by corresponding factors J m , it certainly doesin the case of the Keplerian term. What Breakwell and Vagners [6] propose is thento replace L (cid:48) by the calibrated value (cid:98) L = µ √ (cid:104) − E + (cid:80) km =1 ( J m /m !) H ,m ( L (cid:48) , G (cid:48) , H ) (cid:105) / , (23)14btained by solving the Keplerian term from Eq. (22). If now L (cid:48) is replaced inEq. (22) by the calibrated value (cid:98) L then the energy equation will remain certainlyaccurate to O ( J k +12 ) . Therefore, the initialization of the secular frequencies isnotably improved using the values n (cid:96) = µ (cid:98) L + k (cid:88) m =1 J m m ! ∂ H ,m ∂L (cid:48) , n g = k (cid:88) m =1 J m m ! ∂ H ,m ∂G (cid:48) , n h = k (cid:88) m =1 J m m ! ∂ H ,m ∂H . Obviously, the use of Breakwell and Vagners’ calibration procedure is not con-strained to the second-order of Brouwer’s theory, and also applies to any truncationorder. In our particular case, in which we had already computed the second-orderdirect and inverse periodic corrections of Brouwer’s theory, the calibration of the(mean) Delaunay action allowed us to improve the accuracy of Brouwer’s secu-lar terms to the third order of J without need of computing the long series thatcomprise the third-order term of the generating function.We checked that the new extended third-order Brouwer’s solution, in whichthe second-order corrections consist of a single transformation, enjoys the sameaccuracy as a third-order solution computed in the traditional way of splitting theHamiltonian reduction into the sequence provided by the elimination of the paral-lax, followed by the elimination of the perigee, and ending with a Delaunay nor-malization [3]. At the third order, the initialization of the constants of the later wasalso calibrated in Breakwell and Vagners’ style. In both cases, the required directand inverse transformations were computed in polar variables to avoid singularitiesin the case of circular orbits, but also for efficiency reasons [20, 2, 24, 25].An example of the accuracy obtained with the new single-transformation ap-proach is presented in Fig. 1 for a Topex-type orbit, with a = 7707 . km, e = 0 . , I = 66 . , Ω = 180 . ◦ , ω = 270 ◦ , and (cid:96) = 180 ◦ , yet a varietyof cases have been tested for different types of orbits with analogous results. Thefigure depicts the root sum square (RSS) of the position errors of the analytical so-lution for different truncations of the extended Brouwer’s solution when comparedwith the “true” solution along one month propagation. To guarantee the accuracyof the latter, the true solution was obtained from the numerical integration of thedifferential equations of the main problem in Cartesian coordinates using extendedprecision.Labels {I:S:D} in the plot denote the truncation order of the I nverse correc-tions, S ecular terms, and D irect corrections of the orbital theory. The notation{I+:S:D} means that the inverse corrections are improved in Breakwell and Vagn-ers’ style. Thus, the label {1:2:1} denotes the original Brouwer’s solution (with asingle transformation), which at the end of one month accumulates a RSS error ofabout 2.5 km. The simple calibration of the secular mean motion using Eq. (23),15 }{ + :2:1 } { } { + :3:2 } R SS po s i t i one rr o r ( m ) Figure 1: Root-sum-square error of the Cartesian coordinates provided by differentversions of Brouwer’s solution. Abscissas are days.case {1+:2:1}, clearly bends the RSS errors curve towards the meter level with anegligible increase of the computational burden, reaching less than 20 m at the endof the propagation period. Figures are further improved when the orbit is prop-agated with the full second-order theory, in which the single second-order trans-formation is used both in the initialization of the constants of Brouwer’s secularterms (inverse periodic corrections) and in the ephemeris computation (direct peri-odic corrections). Now, the RSS errors fall below the meter level at day 30, curvelabeled {2:2:2}, yet the computational burden increases by about one third due toevaluation of second-order direct corrections, contrary to the lighter first-order cor-rections. Finally, supplementing Brouwer’s theory with third-order secular termsand the consequent calibration of the secular mean motion, case {2+:2:2}, keepsthe RSS errors to the level of just a few cm along the whole propagation intervalwith negligible increase of the computational burden with respect to the {2:2:2}case.For efficiency in the evaluation of perturbation solutions, arrangement of theseries that comprise the periodic corrections for optimal evaluation is an impor-tant consideration [8, 17]. In this task, we limited our efforts to minor arrange-ments of the code, like the factorization of the inclination polynomials involvedin the different summations and the following use of Horner’s algorithm, and leftthe code optimization job to the compiler. Because we did the same for both an-alytical solutions (Brouwer’s with single periodic corrections, and the traditional16arallax-perigee-Delaunay solution), even if optimal evaluation is not achieved, thecomparisons are not expected to be biased towards a particular theory.After repeated evaluation of the periodic corrections for a variety of initialconditions, we found that the evaluation of the periodic corrections of the tra-ditional analytical solution spends roughly twice the time needed by the single-transformation approach in the evaluation of the periodic corrections. This resultwas a priori unexpected because the series comprising the corrections of the newapproach, which only involves a single transformation, are clearly longer than thecomposition of those involved in the three transformations needed in the traditionalapproach. The improved evaluation efficiency is then attributed to the fact that thecompiler is able to carry out a better optimization of the code in the case of thesingle-transformation approach. This fact may be understood when taking into ac-count that, for instance, the coefficient (5 s − appears about 300 times in ourarrangement of the periodic corrections of the single transformation, but only 73times in the classical parallax-perigee-Delaunay transformation arranged with thesame factorization criterion, where, in particular, this factor only appears in thecorrections related to the elimination of the perigee. Thus, cancelling this commonfactor by the compiler is roughly four times more efficient in the first case. Unde-niably, making a smarter organization of the code before sending it to the compilermight modify these figures. However, the balance is so radical on the side of thesingle transformation that these presumed improvements due to an additional pre-processing of the code are not expected to be relevant enough to revert the figures. Experience gained through the use of Hamiltonian simplification methods promptedus to question Brouwer’s splitting normalization strategy. The convenience of di-viding a normalization process into different stages has been taken for grantedsince the initial efforts in fully automatizing the computation of perturbation the-ories. Needless to say that we agree in which this way of proceeding may easethe construction of the perturbation solution. However, what is not so obvious isthat the evaluation of the solution constructed this way must necessarily yield theless computational burden. On the contrary, results in this paper seem to point inthe direction that the claimed benefits of partial normalization as well as Hamilto-nian simplification procedures can be counterbalanced by other type of considera-tions, at least for the lower orders of normalization that suffice in many practicalcases. Prospective application of the strategy proposed here to other instances ofperturbed Keplerian motion, or to the computation of higher orders of the mainproblem of artificial satellite theory, should contribute to make clear the issue.17rouwer’s closed-form approach and full automatization of the computationof perturbation theories seem two legitimate aims in this epoch of computationalplenty. However, as demonstrated by the equation of the center controversy, ratherthan running perturbation algorithms in batch processes, one should not disregardthe power of modern hand computations carried out with the help of existing soft-ware tools. Indeed, as far as mathematical simplification remains in the categoryof arts, inspection of intermediate expressions turns into a convenient practice thatmay eventually lead the user to straightforward simplifications that make feasibleor just simpler the next step of a partially automated procedure. Like chess players,celestial mechanicians are rarely able to anticipate more than a few moves in theoutcome of a perturbation approach. On the contrary, they need to wait for theopponent’s reaction in order to implement a winning strategy, which, in addition,is most times settled on an empirical basis. It was, in particular, the case of thecurrent research, in which the help provided by the computer algebra system con-verted into a simple task the critical inspection of the seminal solution obtained byBrouwer.
Acknowledgements
Partial support by the Spanish State Research Agency and the European RegionalDevelopment Fund under Projects ESP2016-76585-R and ESP2017-87271-P (AEI/ERDF, EU) is recognized. The author acknowledges with pleasure the help ofSylvio Ferraz-Mello in finding particular passages of volume 2 of Poincaré’s
Méth-odes Nouvelles . References [1] K. Aksnes. A note on ‘The main problem of satellite theory for small eccentricities,by A. Deprit and A. Rom, 1970’.
Celestial Mechanics , 4(1):119–121, September1971.[2] K. Aksnes. On the Use of the Hill Variables in Artificial Satellite Theory.
Astronomyand Astrophysics , 17(1):70–75, February 1972.[3] K. T. Alfriend and S. L. Coffey. Elimination of the perigee in the satellite problem.
Celestial Mechanics , 32(2):163–172, February 1984.[4] D. Boccaletti and G. Pucacco.
Theory of orbits. Volume 2: Perturbative and ge-ometrical methods . Astronomy and Astrophysics Library. Springer-Verlag, BerlinHeidelberg New York, 1st edition, 2002.[5] N. L. Bonavito, S. Watson, and H. Walden. An Accuracy and Speed Comparisonof the Vinti and Brouwer Orbit Prediction Methods. Technical Report NASA TND-5203, Goddard Space Flight Center, Greenbelt, Maryland, May 1969.
6] J. V. Breakwell and J. Vagners. On Error Bounds and Initialization in Satellite OrbitTheories.
Celestial Mechanics , 2:253–264, June 1970.[7] D. Brouwer. Solution of the problem of artificial satellite theory without drag.
TheAstronomical Journal , 64:378–397, November 1959.[8] S. Coffey and A. Deprit. Fast evaluation of Fourier series.
Astronomy and Astro-physics , 81:310–315, January 1980.[9] A. Deprit. Canonical transformations depending on a small parameter.
CelestialMechanics , 1(1):12–30, 1969.[10] A. Deprit. The elimination of the parallax in satellite theory.
Celestial Mechanics ,24(2):111–153, 1981.[11] A. Deprit and B. Miller. Simplify or Perish.
Celestial Mechanics , 45:189–200, 1989.[12] A. Deprit and A. Rom. The Main Problem of Artificial Satellite Theory for Smalland Moderate Eccentricities.
Celestial Mechanics , 2(2):166–206, June 1970.[13] E. Deprit and A. Deprit. Poincaré’s méthode nouvelle by skew composition.
CelestialMechanics and Dynamical Astronomy , 74(3):175–197, July 1999.[14] S. Ferraz-Mello. Do Average Hamiltonians Exist?
Celestial Mechanics and Dynam-ical Astronomy , 73:243–248, January 1999.[15] S. Ferraz-Mello.
Canonical Perturbation Theories - Degenerate Systems and Reso-nance , volume 345 of
Astrophysics and Space Science Library . Springer, New York,January 2007.[16] L. M. Healy. The Main Problem in Satellite Theory Revisited.
Celestial Mechanicsand Dynamical Astronomy , 76(2):79–120, 2000.[17] L. M. Healy and J. J. Travisano. Automatic rendering of astrodynamics expressionsfor efficient evaluation.
Journal of the Astronautical Sciences , 46(1):65–81, 1998.[18] F. R. Hoots and R. L. Roehrich. Models for Propagation of the NORAD ElementSets. Project SPACETRACK, Rept. 3, U.S. Air Force Aerospace Defense Command,Colorado Springs, CO, December 1980.[19] G.-i. Hori. Theory of General Perturbation with Unspecified Canonical Variables.
Publications of the Astronomical Society of Japan , 18(4):287–296, 1966.[20] I. G. Izsak. A note on perturbation theory.
The Astronomical Journal , 68(8):559–561,October 1963.[21] W. H. Jefferys. Automated, Closed Form Integration of Formulas in Elliptic Motion.
Celestial Mechanics , 3:390–394, September 1971.[22] Y. Kozai. Mean values of cosine functions in elliptic motion.
The AstronomicalJournal , 67:311, June 1962.[23] Y. Kozai. Second-order solution of artificial satellite theory without air drag.
TheAstronomical Journal , 67:446–461, September 1962.
24] M. Lara. Efficient Formulation of the Periodic Corrections in Brouwer’s GravitySolution.
Mathematical Problems in Engineering , vol. 2015(Article ID 980652):1–9, 2015.[25] M. Lara. LEO intermediary propagation as a feasible alternative to Brouwer’s gravitysolution.
Advances in Space Research , 56(3):367–376, August 2015.[26] M. Lara. Review of analytical solutions for low earth orbit propagation and study ofthe precision improvement in the conversion of osculating to mean elements. Techni-cal Report CM 2019/SER0023, Universidad de La Rioja, Logroño, La Rioja, Septem-ber 2019.[27] M. Lara, J. F. San-Juan, Z. J. Folcik, and P. Cefola. Deep Resonant GPS-DynamicsDue to the Geopotential.
Journal of the Astronautical Sciences , 58(4):661–676, Oc-tober 2011.[28] M. Lara, J. F. San-Juan, and L. M. López-Ochoa. Delaunay variables approach tothe elimination of the perigee in Artificial Satellite Theory.
Celestial Mechanics andDynamical Astronomy , 120(1):39–56, September 2014.[29] M. Lara. A new radial, natural, higher-order intermediary of the main problem fourdecades after the elimination of the parallax.
Celestial Mechanics and DynamicalAstronomy , 131(9):20, September 2019.[30] M. Lara. Solution to the main problem of the artificial satellite by reverse normaliza-tion.
Nonlinear Dynamics , 101(2):1501–1524, July 2020.[31] R. H. Lyddane and C. J. Cohen. Numerical comparison between Brouwer’s theoryand solution by Cowell’s method for the orbit of an artificial satellite.
AstronomicalJournal , 67:176–177, April 1962.[32] G. Metris. Mean values of particular functions in the elliptic motion.
Celestial Me-chanics and Dynamical Astronomy , 52:79–84, March 1991.[33] G. Métris and P. Exertier. Semi-analytical theory of the mean orbital motion.
Astron-omy and Astrophysics , 294:278–286, February 1995.[34] K. R. Meyer, G. R. Hall, and D. Offin.
Introduction to Hamiltonian DynamicalSystems and the N-Body Problem , volume 90 of
Applied Mathematical Sciences .Springer, New York, 2nd edition, 2009.[35] J.A. Morrison. Generalized method of averaging and the von zeipel method. InRaynor L. Duncombe and Victor G. Szebehely, editors,
Methods in Astrodynamicsand Celestial Mechanics , volume 17 of
Progress in Astronautics and Rocketry , pages117 – 138. Elsevier, 1966.[36] C. Osácar and J. F. Palacián. Decomposition of functions for elliptical orbits.
Celes-tial Mechanics and Dynamical Astronomy , 60(2):207–223, October 1994.[37] Henri Poincaré.
Les méthodes nouvelles de la mécanique céleste. Tome 2 . Gauthier-Villars et fils (Paris), 1893.
38] J. F. San-Juan, D. Ortigosa, L. M. López-Ochoa, and R. López. Deprit’s Eliminationof the Parallax Revisited.
Journal of the Astronautical Sciences , 60:137–148, June2013.[39] D. Steichen. An averaging method to study the motion of lunar artificial satellitesII: Averaging and Applications.
Celestial Mechanics and Dynamical Astronomy ,68(3):225–247, July 1998.[40] F. Tisserand.
Traité de mécanique céleste. Tome I: Perturbations des planètes d’aprésla méthode de la variation des constantes arbitraries . Gauthier-Villars et fils, Quaides Grands-Augustins, 55, Paris, 1889.[41] D. A. Vallado, P. Crawford, R. Hujsak, and T. S. Kelso. Revisiting Spacetrack Re-port
AIAA/AAS Astrodynamics Specialist Conference andExhibit , pages 1–88, USA, August 2006. Guidance, Navigation, and Control and Co-located Conferences, American Institute of Aeronautics and Astronautics.[42] H. von Zeipel. Research on the Motion of Minor Planets. NASA TT F-9445, 1965.(NASA Translation of: Recherches sur le mouvement des petites planètes, Arkiv förmatematik, astronomi och fysik, vol. 11, 1916, vol. 12, 1917, vol. 13, 1918)., pages 1–88, USA, August 2006. Guidance, Navigation, and Control and Co-located Conferences, American Institute of Aeronautics and Astronautics.[42] H. von Zeipel. Research on the Motion of Minor Planets. NASA TT F-9445, 1965.(NASA Translation of: Recherches sur le mouvement des petites planètes, Arkiv förmatematik, astronomi och fysik, vol. 11, 1916, vol. 12, 1917, vol. 13, 1918).