Analytical theory of an artificial satellite of the moon

See discussions, stats, and author profiles for this publication at: https://www.researchgate.net/publication/8488512 Analytical Theory of an Artificial Satellite of the Moon Article in Annals of the New York Academy of Sciences · June 2004 DOI: 10.1196/annals.1311.025 · Source: PubMed CITATIONS READS 3 45 2 authors, including: Bernard De Saedeleer Université Catholique de Louvain 16 PUBLICATIONS 94 CITATIONS SEE PROFILE All content following this page was uploaded by Bernard De Saedeleer on 31 December 2016. The user has requested enhancement of the downloaded file. All in-text references underlined in blue are added to the original document and are linked to publications on ResearchGate, letting you access and read them immediately. Analytical Theory of an Artificial Satellite of the Moon BERNARD DE SAEDELEER AND JACQUES HENRARD Département de Mathématique, University of Namur, Namur, Belgium ABSTRACT: The dynamics of an artificial satellite of the moon is quite different from the dynamics of an artificial satellite of the Earth. Indeed, the C22 term is only 1/10 of the J2 term, and the effect of the Earth on the lunar satellite is much larger than the effect of the Moon on a terrestrial satellite. The method used here is the Lie method for averaging the Hamiltonian of the problem, in canonical variables. The solution is developed in powers of the small factors linked to J2 and C22. Short period terms (linked to l, the mean anomaly) are eliminated first, and then the long period terms (linked to g, the mean argument of the periaster, and to h, the longitude of the ascending node), which finally gives the secular motion. The results are obtained in a closed form, without any series developments in eccentricity or inclination. Thus, the solution applies for a wide range of values, except for few isolated critical values. The results are only very preliminary. As a side result, we were able to check the solution given by Kozai for the effect of the J2 term on an artificial satellite. KEYWORDS: lunar artificial satellite; non zonal; Lie; Hamiltonian; C22 INTRODUCTION The lunar gravity field is far from being central. The first lunar mascons were discovered in 1968,1 and the non-centrality of the gravity field was confirmed later by more precise measurements and models of increasing order2–4 in spherical harmonics (see, e.g., Ref. 5 for a summary). The most recent model nowadays is that due to Konopliv et al.4 Normalized coefficients are considered, to account for their relative importance. The order of magnitudes of the second order coefficients for the Earth (see Ref.6 p.155) and the Moon (see Ref. 2, p.1018) are given in TABLE 1. The Moon is much less flattened than the Earth, which makes the value of the C22 coefficient closer to that of J2 (one order of magnitude instead of three in the case of the Earth); our aim is to introduce the effect of C22 in addition to that of J2. THE MAIN PROBLEM OF A LUNAR ARTIFICIAL SATELLITE The case of a satellite around the Moon is quite different from that of a satellite around the Earth because of the following facts: • the C22 term is only 1/10 of the J2 term, so effects should be considered; Address for correspondence: Bernard De Saedeleer, Département de Mathématique, University of Namur, Rempart de la Vierge, 8, B-5000 Namur, Belgium. [email protected] Ann. N.Y. Acad. Sci. 1017: 434–449 (2004). ©2004 New York Academy of Sciences. doi: 10.1196/annals.1311.025 DE SAEDELEER & HENRARD: AN ARTIFICIAL SATELLITE OF THE MOON 435 TABLE 1. Orders of magnitude for J2 and C22 C20 ≡ J2 C22 Earth ⊕ − 10−3 2 × 10−6 Moon (| − 2 × 10−4 2 × 10−5 • the effect of the Earth on the lunar satellite is much larger than the effect of the Moon on a terrestrial satellite, so the former effect is mixed effects of the shape of the lunar gravity field; • the moon is a slowly rotating body; and • the moon has no dense atmosphere. In this preliminary study, we make the following assumptions: • the orbit of the Moon is circular (we neglect the eccentricity of 0.055); • the motion of the Moon is uniform (librations are neglected); • the lunar equator lies in the ecliptic (we neglect the inclination of the lunar equator to the ecliptic of 1.5° and the inclination of the lunar orbit to the ecliptic of about 5°); • the perturbation of the Sun is negligible; and • the longitude of the lunar longest meridian λ22 is equal to the longitude of the Earth λ⊕ (librations are neglected). The last assumption is one of the well known Cassini's laws,7 stating that the Moon is always presenting the same face to us (synchronous rotation). This means that the Moon rotates about its axis perpendicular to the plane of its orbit at an angular velocity γ(| that is equal to the mean angular velocity in the orbit n(|; that is to say, γ(| = n(| = 2π/T (| , with the sidereal rotation period of the Moon T (| = 27.3216615solar days. In this paper, we focus on the combined effect of the perturbations J2 and C22. The effect of the Earth will be introduced in a forthcoming paper. PARTIAL PERTURBATIVE HAMILTONIANS IN J2 AND C22 Many efficient perturbation methods are available within the frame of the Hamiltonian formalism; we use here the method of the Lie transform8 (or the Lie triangle). The Hamiltonian must be written in powers of a perturbative parameter ε: H(0) = ei (0) ---- H i ; i ≥ 0 i! ∑ the unperturbed potential being H 0( 0 ) . In the case of the spherical harmonic description of the gravitational potential (see Ref. 6 for the mathematical foundations), we have only first powers in the perturbative coefficients Cmn and Smn, thus we should write: H ( 0 ) = H 0( 0 ) + (0). ∑ εmn Hmn m, n 436 ANNALS NEW YORK ACADEMY OF SCIENCES In the present description, we only use J2 and C22, so that we can even write: H ( 0 ) = H 0( 0 ) + εH 1( 0 ) + δHB 1( 0 ) if we define ε = J2R2 and δ = − C22R2. To retain the Hamiltonian formalism, it is required to work in canonical variables; we choose the classical Delaunay variables (qi, pi) = (l, g, h, L, G, H) defined by l = E – e sin E L = µa g = ω G = µa ( 1 – e 2 ) h = Ω H = µa ( 1 – e 2 ) cos I . In these variables, the central term is written − µ2/2L2 and the Hamiltonian equations of motion are simply: ∂H ∂H L̇ = – --------l˙ = --------∂l ∂L ∂H ∂H Ġ = – --------ġ = --------(1) ∂g ∂G ∂H ∂H Ḣ = – --------- . ḣ = --------∂h ∂H Now we need to write the partial perturbative Hamiltonians H 1( 0 ) and HB 1( 0 ) in these variables and in an inertial frame; that is to say, with respect to a constant direction in space. We define our inertial frame (x, y, z) as follows (see FIGURE 1): the origin is taken at the center of the Moon; the x direction is toward the first point of Aries ϒ, the y is the direction normal to x and contained in the lunar equatorial plane containing x, and the z direction is the right-handed normal to (x, y). To be able to use the expressions of the spherical harmonics for the potential, we first need to define spherical coordinates (r, λ′, φ), so that the longitude of the satellite λ′ starts from the x-axis in the equatorial plane, the latitude φ being defined as the deviation from the (x, y) plane. Within the inertial frame, the Hamiltonian is written H = T + V = T + V00 + V20 − V22, where µ V 00 = – --r 1 P 20(x) = --- ( 3x 2 – 1 ) µ R 2 2 V 20 = --- ⎛ ---⎞ J 2 P 20 ( sin φ ) and r ⎝ r⎠ P 22(x) = 3 ( 1 – x 2 ). µ R 2 V 22 = --- ⎛ ---⎞ C 22 P 22 ( sin φ ) cos ( 2 ( λ′ – λ 22 ) ) r ⎝ r⎠ The Hamiltonian includes the central field term V00 and the perturbative potentials V20 and V22. Use has been made of the Legendre associated functions P20 and P22. We can translate their argument (sin φ) into Delaunay variables by the way of spherical trigonometry (see FIGURE 1, where the plane of the orbit is at an inclination I), sin φ = sin I sin(f + g). Then, µ V 20 = εH 1( 0 ) = – ( J 2 R 2 ) -------- ( 1 – 3 sin2 I sin2 ( f + g ) ), 2r 3 so that DE SAEDELEER & HENRARD: AN ARTIFICIAL SATELLITE OF THE MOON 437 µ H 1( 0 ) = – -------- ( 1 – 3 sin2 I sin2 ( f + g ) ). 2r 3 Expanding to the first power only in the trigonometric functions finally gives µ H 1( 0 ) = -------- ( 1 – 3c 2 – 3s 2 cos ( 2 f + 2g ) ). (2) 4r 3 However, the coefficient C22 makes the longitude λ′ appear in addition to the latitude φ. Since the spherical harmonics are defined with respect to the main axis of inertia of the attracting body, we need to define λ22 as the longitude of the lunar longest meridian (minimum inertia). This angle causes the Hamiltonian to be timedependent, since λ22 = λ⊕ travels at the rate of the synchronous rotation, λ̇ ⊕ = n(|. To eliminate this dependency, we work in a rotating system whose x′ axis now passes through the Earth; we then define new longitudes with respect to λ⊕, λ = λ′ − λ⊕, and we also redefine h = Ω − λ⊕ (the angles always appear in this combination). The Hamiltonian H needs to be redefined in order to take this rotation into account, H ( 0 ) = H 0( 0 ) + εH 1( 0 ) + δHB 1( 0 ) – n (| H . (3) The equations of motion are the same as in (1), we only changed slightly the content of H. Now HB 1( 0 ) can be written µ HB 1( 0 ) = ----- P 22( sin φ) cos ( 2λ ). (4) r3 FIGURE 1. Simplified selenocentric sphere. 438 ANNALS NEW YORK ACADEMY OF SCIENCES Once again, the factor cos(2λ) can be translated into Delaunay variables in the same way, using spherical trigonometry (see FIG. 1). It is now useful to define the shortcuts (s, c) = (sin I, cos I) and (X, Y) = (cos(f + g), sin(f + g)), so that we have the following spherical relationships: sin φ = sY X = cos φ cos ( λ – h ) cos φ = cY sin ( λ – h ) + X cos ( λ – h ) . (5) Solving (5) for cos(2λ), [ 1 + Y 2 ( s 2 – 2 ) ] cos ( 2h ) – 2cXY sin ( 2h ) cos ( 2λ ) = ---------------------------------------------------------------------------------------------------- . cos2 φ (6) It turns out that the denominator cos2 φ appearing in (6) cancels when multiplied by P22(sin φ) = 3 cos2 φ, giving the following expression for HB 1( 0 ) , after gathering the trigonometric arguments and reducing the trigonometric terms to first power only, 3µ HB 1( 0 ) = -------- { 2s 2 cos ( 2h ) + ( c + 1 ) 2 cos ( 2 f + 2g + 2h ) 4r 3 + (c – 1 ) 2 cos ( 2 f (7) + 2g – 2h ) }. At this stage, there remains in (2) and (7) only r to be expressed as a function of (l, g, h) in order to be able to apply a perturbation method. It turns out that the functions r = r(l, g, h) and f = f(l, g, h) cannot be expressed in a closed form and that one usually falls back at this point into series development in the eccentricity. We wish to avoid this, at least for the following reasons: the results would be much less compact and, hence, lack ease of interpretation; moreover they would no longer be valid for higher values of the eccentricity. Instead, we prefer to use the following set of auxiliary variables (ξ, f, g, h, a, n, e, η, s, c), which is closed, and allow for high eccentricities 1 + e cos f a 1 ξ = --- = ------------------------ = -----------------------2 r 1 – e cos E 1–e f L2 a = -----µ µ2 n = -----L3 e = G 2 1 – ⎛ ----⎞ ⎝ L⎠ s = sin I = g H 2 1 – ⎛ ----⎞ ⎝ G⎠ η = G 1 – e 2 = ---L (8) H c = cos I = ---G h The only drawback of this set is that it is redundant and we need to perform partial derivatives with respect to the canonical variables; but this not too heavy a task. We have, for example, ξ 2 e sin f ⎞ ∂A 2 dA ∂A ∂A ∂ξ ∂A ∂ f ∂A ∂A ⎛ –---------------------- + ( ξ η ). + + = = + ⎝ ⎠ ∂f η dl ∂l ∂ξ ∂l ∂ f ∂l ∂l ∂ξ (9) DE SAEDELEER & HENRARD: AN ARTIFICIAL SATELLITE OF THE MOON 439 The partial derivatives of the auxiliary variables (8) with respect to the canonical variables (l, g, h, L, G, H) are summarized in TABLE 2. Computation of the partial derivatives themselves requires some caution, but is not too complicated. Use has to be made of the Kepler equation (l = E − e sinE) that links the anomalies. One might think that the complexity introduced by the partial derivatives, such as (9), is very penalizing, since there are many products. It turns out that it only increases the number of terms to be treated (an order of magnitude in the size of the computations is given in a later section), but that it does not introduce additional problems for the integration itself. As a matter of fact, in all the derivatives, ξ only appears in the numerator, so that the only problem that could be critical in being able to integrate is to have ξ1 or ξ0 finally (since one has to divide out the terms by ξ2 to switch the integration from l into f). In this new set of variables, (8), the factor µ/r3 appearing in (2) and (7) is written 3 ξ n2, which allows us to express the perturbations in the following final compact form: ξ3n2 H 1( 0 ) = ----------- ( 1 – 3c 2 – 3s 2 cos ( 2 f + 2g ) ) 4 ξ3n2 HB 1( 0 ) = -----------3 { 2s 2 cos ( 2h ) + ( c + 1 ) 2 cos ( 2 f + 2g + 2h ) 4 (10) + ( c – 1 ) 2 cos ( 2 f + 2g – 2h ) }. TABLE 2. Table of partial derivatives ∂/∂L ∂/∂G ∂/∂H ∂/∂l ξ ξ2 η2 ------------ cos f na 2 e ξ2 η – ------------ cos f na 2 e 0 ξ2 e – -------- sin f η a 2 -----an 0 0 0 n 3 – ----a2 0 0 0 s 0 c2 ---------------na 2 ηs c – ---------------na 2 ηs 0 c 0 c – ------------na 2 η 1 ------------na 2 η 0 e η2 -----------na 2 e η – -----------na 2 e 0 0 η η – --------na 2 1 --------na 2 0 0 f 1 + ξη 2 ------------------- sin f na 2 e 1 + ξη 2 – ------------------- sin f ηna 2 e 0 ξ2 η 440 ANNALS NEW YORK ACADEMY OF SCIENCES PERTURBATION METHOD USING SEVERAL PARAMETERS As perturbation technique, we use the Lie transform.8 The essence of the Lie transform is to build a canonical transformation analytic in ε at ε = 0 so as to achieve a specific requirement in the new transformed Hamiltonian; that is, we choose to suppress the fastest angle. The main advantages are that it is done explicitly (whereas the von Zeipel method is implicit), in a practical recursive form, and that the inverse is constructed in the same way. The Hamiltonians are written in powers of a perturbative parameter ε. For the initial Hamiltonian (input), we write εi (0) ---- H i ; i ≥ 0 i! H(0) = ∑ and for the transformed Hamiltonian (output), where the fast angle has been removed, we write H0 = εi ---- H 0( 0 ) . i ≥ 0 i! ∑ This transformation is symbolized by the Lie triangle, where the first column represents the input Hamiltonian, and the diagonal the output Hamiltonian, H 0( 0 ) H 1( 0 ) H 0( 1 ) H 2( 0 ) H 1( 1 ) H 0( 2 ) ⋮ ⋮ ⋮ ⋱ The triangle is filled by the way of the following recursive formula, H i( j ) = H i( +j –11 ) + i ∑ C ik ( H i( –j –k 1 ) ;W k + 1 ) , k=0 where the Wk are generating functions. For example, we have, for the first order, H 0( 1 ) = H 1( 0 ) + ( H 0( 0 ) ;W 1 ) and for the second order H 0( 2 ) = H 1( 1 ) + ( H 0( 1 ) ;W 1 ) and H 1( 1 ) = H 2( 0 ) + ( H 1( 0 ) ;W 1 ) + ( H 0( 0 ) ;W 2 ). In the case of spherical harmonics, which provides a linear description in Cmn, all H k( 0≥) 2 are zero, and we can write more precisely H 0( 2 ) = ( H 1( 0 ) + H 0( 1 ) ;W 1 ) + ( H 0( 0 ) ;W 2 ) . (11) In these recurrences, use has to be made of the Poisson parenthesis, defined in our canonical variables by ∂A ∂B ∂B ∂A ∂A ∂B ∂B ∂A ∂A ∂B ∂B ∂A (12) – – + + – . ∂l ∂L ∂l ∂L ∂g ∂G ∂g ∂G ∂h ∂H ∂h ∂H As noted previously, when working in the new compact set of variables (ξ, f, g, h, a, n, e, η, s, c) we need to compute their partial derivatives with respect to (l, g, h, L, G, H) in order to compute the Poisson parenthesis; these are given in the TABLE 2. In the present description, the perturbative parameter is not only ε, but also δ (see (3)). Thus, we must rewrite the recursive formula of the first order as ( A ;B ) = DE SAEDELEER & HENRARD: AN ARTIFICIAL SATELLITE OF THE MOON εH 0( 1 ) + δHB 0( 1 ) = εH 1( 0 ) + δHB 1( 0 ) + ( H 0( 0 ) ;W 1* ) . 441 (13) The Lie algorithm consists of imposing the form of the new Hamiltonian. Here, we choose it to be independent of the fastest angle l, by making the following choice: 1 2π εH 0( 1 ) + δHB 0( 1 ) = ------ ∫ ( εH 1( 0 ) + δHB 1( 0 ) ) dl 2π 0 (14) = … = ( εH 0( 1 ) + δHB 0( 1 ) ). We write H 0( 1 ) as H 0( 1 ) in order to remember that l has been eliminated. Then we use (13) to compute the generator W 1* . H 0( 0 ) being a function of L only, the Poisson parenthesis reduces to one term, giving 1 l W 1* = --- ∫ ( εH 1( 0 ) + δHB 1( 0 ) – εH 0( 1 ) – δHB 0( 1 ) ) dl n 0 (15) = … = εW 1 + δWB 1 . Note that, for the case δ = 0, we should retrieve the results for the effect of J2 alone; this is of course so, as is shown in the next section. For the second order, we must also rewrite (11) as H 0( 2 ) = ( εH 1( 0 ) + δHB 1( 0 ) + εH 0( 1 ) + δHB 0( 1 ) ;W 1* ) + ( H 0( 0 ) ;W 2* ), so that H 0( 2 ) and W 2* contain second order terms; that is, terms in ε2, (16) εδ, and δ2. COMPUTATIONS, TOOL, AND RESULTS The number of computations is quite large, thus, hand calculations are avoided. For a first attempt, we used symbolic software, but on one hand the results are not so readable (because the results are often displayed in several pages of large terms for a single expression), and on the other hand the computational speed is quite large. Thus, we preferred to use a specific Fortran code that was clearer and faster; the typical gain in CPU time is two orders of magnitude. This tool is called the MM, standing for Moon’s series manipulator. It was developed at our University. In this tool, each expression is given by a series of linear trigonometric functions, with polynomial coefficients. That is why we preferred to reduce the perturbation to the linear forms (10). An example of such a series is given in TABLE 3. When multiplying such series, more and more terms are produced and they are always reduced to linear form by the following well-known relationships: TABLE 3. The series defining ␰ = (1 + e cos f)/(1 – e2) f g h ξ a n e η c s (f − l) nδ cos 0 0 0 0 0 0 0 −2 0 0 0 0 0.1000000000000000D +01 cos 1 0 0 0 0 0 1 −2 0 0 0 0 0.1000000000000000D +01 coefficient 442 ANNALS NEW YORK ACADEMY OF SCIENCES TABLE 4. The initial perturbations εH 1( 0 ) + δHB 1( 0 ) f g h ξ a n e η c s (f − l) nδ cos 0 0 0 3 0 2 0 0 0 0 0 0 cos 0 0 0 3 0 2 0 0 2 0 0 0 − 0.7500000000000000D +00 cos 2 2 0 3 0 2 0 0 0 2 0 0 − 0.7500000000000000D +00 cos 0 0 2 3 0 2 0 0 0 0 0 1 cos 0 0 2 3 0 2 0 0 2 0 0 1 − 0.1500000000000000D +01 cos 2 2 2 3 0 2 0 0 0 0 0 1 0.7500000000000000D +00 cos 2 2 2 3 0 2 0 0 1 0 0 1 0.1500000000000000D +01 cos 2 2 2 3 0 2 0 0 2 0 0 1 0.7500000000000000D +00 cos 2 2 −2 3 0 2 0 0 0 0 0 1 0.7500000000000000D +00 cos 2 2 −2 3 0 2 0 0 1 0 0 1 − 0.1500000000000000D +01 cos 2 2 −2 3 0 2 0 0 2 0 0 1 coefficient 0.2500000000000000D +00 0.1500000000000000D +01 0.7500000000000000D +00 sin ( A + B ) + sin ( A – B ) sin A cos B = ----------------------------------------------------------2 sin ( A + B ) – sin ( A – B ) cos A sin B = ----------------------------------------------------------2 cos ( A + B ) + cos ( A – B ) cos A cos B = ------------------------------------------------------------2 cos ( A – B ) – cos ( A + B ) sin A sin B = ------------------------------------------------------------- . 2 The concern to keep a linear expression is motivated by the fact that we want to retain ease in the integrations. The series of the perturbation in the parameters J2 and C22 written in (10) is given in TABLE 4. Note that the variable nδ represent the exponent in the perturbative parameter δ. For an expression of order p, we then have terms containing factors like ε p – nδ δ nδ . Note also that we retain the variable ξ as such for as long as possible in order to have a compact form for the series. When integrating the expression with respect to l, we simply divide by ξ2η (since dl = df/ξ2η) in order to integrate with respect to f, and at this stage we only replace ξ by its own series ξ(f), defined in TABLE 3. TABLE 5. The first-order averaged Hamiltonian εH 0( 1 ) + δHB 0( 1 ) f g h ξ a n e η c s (f − l) nδ cos 0 0 0 0 0 2 0 −3 0 0 0 0 cos 0 0 0 0 0 2 0 −3 2 0 0 0 − 0.7500000000000000D +00 cos 0 0 2 0 0 2 0 −3 0 0 0 1 cos 0 0 2 0 0 2 0 −3 2 0 0 1 − 0.1500000000000000D +01 coefficient 0.2500000000000000D +00 0.1500000000000000D +01 DE SAEDELEER & HENRARD: AN ARTIFICIAL SATELLITE OF THE MOON 443 The result of the computation (14) of the first-order averaged Hamiltonian is given in TABLE 5, and the corresponding first order generator W 1* computed from (15) is given in Table 6. TABLE 6. The first order generator W 1* = εW 1 + δWB 1 f g h ξ a n e η c s (f − l) nδ coefficient cos 0 0 0 0 0 1 0 −3 0 0 1 0 0.2500000000000000D +00 cos 0 0 0 0 0 1 0 −3 2 0 1 0 − 0.7500000000000000D +00 sin 1 0 0 0 0 1 1 −3 0 0 0 0 sin 1 0 0 0 0 1 1 −3 2 0 0 0 − 0.7500000000000000D +00 0.2500000000000000D +00 sin 1 2 0 0 0 1 1 −3 0 2 0 0 − 0.3750000000000000D +00 sin 2 2 0 0 0 1 0 −3 0 2 0 0 − 0.3750000000000000D +00 sin 3 2 0 0 0 1 1 −3 0 2 0 0 − 0.1250000000000000D +00 cos 0 0 2 0 0 1 0 −3 0 0 1 1 cos 0 0 2 0 0 1 0 −3 2 0 1 1 − 0.1500000000000000D +01 0.1500000000000000D +01 sin 1 0 2 0 0 1 1 −3 0 0 0 1 sin 1 0 2 0 0 1 1 −3 2 0 0 1 − 0.7500000000000000D +00 sin 1 0 −2 0 0 1 1 −3 0 0 0 1 0.7500000000000000D +00 0.7500000000000000D +00 sin 1 0 −2 0 0 1 1 −3 2 0 0 1 − 0.7500000000000000D +00 sin 1 2 2 0 0 1 1 −3 0 0 0 1 0.3750000000000000D +00 sin 1 2 2 0 0 1 1 −3 1 0 0 1 0.7500000000000000D +00 sin 1 2 2 0 0 1 1 −3 2 0 0 1 0.3750000000000000D +00 sin 1 2 −2 0 0 1 1 −3 0 0 0 1 0.3750000000000000D +00 sin 1 2 −2 0 0 1 1 −3 1 0 0 1 − 0.7500000000000000D +00 sin 1 2 −2 0 0 1 1 −3 2 0 0 1 0.3750000000000000D +00 sin 2 2 0 0 1 0 −3 0 0 0 1 0.3750000000000000D +00 2 sin 2 2 2 0 0 1 0 −3 1 0 0 1 0.7500000000000000D +00 sin 2 2 2 0 0 1 0 −3 2 0 0 1 0.3750000000000000D +00 sin 2 2 −2 0 0 1 0 −3 0 0 0 1 0.3750000000000000D +00 sin 2 2 −2 0 0 1 0 −3 1 0 0 1 − 0.7500000000000000D +00 sin 2 2 −2 0 0 1 0 −3 2 0 0 1 0.3750000000000000D +00 sin 3 2 0 0 1 1 −3 0 0 0 1 0.1250000000000000D +00 2 sin 3 2 2 0 0 1 1 −3 1 0 0 1 0.2500000000000000D +00 sin 3 2 2 0 0 1 1 −3 2 0 0 1 0.1250000000000000D +00 0.1250000000000000D +00 sin 3 2 −2 0 0 1 1 −3 0 0 0 1 sin 3 2 −2 0 0 1 1 −3 1 0 0 1 − 0.2500000000000000D +00 sin 3 2 −2 0 0 1 1 −3 2 0 0 1 0.1250000000000000D +00 444 ANNALS NEW YORK ACADEMY OF SCIENCES After having performed the computation, the first-order averaged Hamiltonian can be rewritten in full, 3n 2 n2 εH 0( 1 ) + δHB 0( 1 ) = ε --------- ( 1 – 3c 2 ) + δ --------- ( 2s 2 cos ( 2h ) ) . 4η 3 4η 3 (17) This form was expected: looking at (10) and considering that the constant part of ξ3, after dividing by ξ2η, gives finally η−3. Of course, in the case δ = 0 we retrieve the known results for the effect of J2 alone. The expression for εH 0( 1 ) is the same as that given by ( – F 1* ) defined in Equation 13 of Reference 9, when setting 2k2 = ε. Moreover, when expressing the averaged Hamiltonian as H = H 0( 0 ) ( L ) + εH 0( 1 ) ( L, G, H ) + O(ε 2), (18) we can deduce the first order averaged equations of motion ∂H 0( 1 ) ∂H 0( 0 ) l˙ = --------------- + ε --------------∂L ∂L L˙ = 0 ∂H 0( 1 ) ġ = 0 + ε --------------∂G G˙ = 0 ∂H 0( 1 ) h˙ = 0 + ε --------------∂H H˙ = 0 (19) which can then be written in the usual Keplerian elements as ˙ n 3 l = n + ε – --- --------------------------4 a2( 1 – e2 )2 1 – e 2 ( 1 – 3 cos2 I ) n 3 ġ = ε – --- --------------------------- ( 1 – 5 cos2 I ) 4 a2( 1 – e2 )2 (20) n 3 h˙ = ε – --- --------------------------- ( 2 cos2 I ). 4 a2( 1 – e2 )2 The two formulæ giving the effect of J2 on g and h are well known.10–12 The peculiar value of the inclination that causes ġ to vanish, known as the critical inclination Ic = 63°26′, is also well known10 and even led to some controversy in the past.12 These two formulæ can also be derived from the famous Lagrange planetary equations.6,7 The first order generator, whose series is given in TABLE 6, can also be rewritten in full, W 1* = εW 1 + δWB 1 , (21) with n W 1 = --------- [ 2 ( 1 – 3c 2 ) ( f – l ) + 2e ( 1 – 3c 2 ) sin f 8η 3 – and 3es 2 sin ( f + 2g ) – 3s 2 sin ( 2 f + 2g ) – (22) es 2 sin ( 3 f + 2g ) ] DE SAEDELEER & HENRARD: AN ARTIFICIAL SATELLITE OF THE MOON 445 TABLE 7. The second-order averaged Hamiltonian (multiplied by F) cos f g h ξ a n e η c s (f − l) nδ 0 0 0 0 0 0 0 0 2 0 0 0 coefficient 0.8000000000000000D +01 cos 0 0 0 0 0 0 0 0 4 0 0 0 − 0.4000000000000000D +02 cos 0 0 0 0 0 0 0 1 0 0 0 0 − 0.4000000000000000D +01 0.2400000000000000D +02 cos 0 0 0 0 0 0 0 1 2 0 0 0 cos 0 0 0 0 0 0 0 1 4 0 0 0 − 0.3600000000000000D +02 cos 0 0 0 0 0 0 2 0 0 0 0 0 0.5000000000000000D +01 cos 0 0 0 0 0 0 2 0 2 0 0 0 − 0.1800000000000000D +02 cos 0 0 0 0 0 0 2 0 4 0 0 0 cos 0 2 0 0 0 0 2 0 0 2 0 0 − 0.2000000000000000D +01 cos 0 2 0 0 0 0 2 0 2 2 0 0 cos 0 0 2 0 0 0 0 0 0 2 0 1 − 0.1280000000000000D +03 0.5000000000000000D +01 0.3000000000000000D +02 0.1600000000000000D +03 cos 0 0 2 0 0 0 0 0 2 2 0 1 cos 0 0 2 0 0 0 0 1 0 2 0 1 − 0.4800000000000000D +02 cos 0 0 2 0 0 0 0 1 2 2 0 1 0.1440000000000000D +03 cos 0 0 2 0 0 0 2 0 0 2 0 1 − 0.5200000000000000D +02 cos 0 0 2 0 0 0 2 0 2 2 0 1 − 0.2000000000000000D +02 cos 0 2 2 0 0 0 2 0 0 0 0 1 − 0.4000000000000000D +01 cos 0 2 2 0 0 0 2 0 1 0 0 1 0.5200000000000000D +02 cos 0 2 2 0 0 0 2 0 2 0 0 1 0.5600000000000000D +02 cos 0 2 2 0 0 0 2 0 3 0 0 1 − 0.6000000000000000D +02 cos 0 2 2 0 0 0 2 0 4 0 0 1 − 0.6000000000000000D +02 cos 0 2 −2 0 0 0 2 0 0 0 0 1 − 0.4000000000000000D +01 cos 0 2 −2 0 0 0 2 0 1 0 0 1 − 0.5200000000000000D +02 cos 0 2 −2 0 0 0 2 0 2 0 0 1 0.5600000000000000D +02 cos 0 2 −2 0 0 0 2 0 3 0 0 1 0.6000000000000000D +02 cos 0 2 −2 0 0 0 2 0 4 0 0 1 − 0.6000000000000000D +02 cos 0 0 0 0 0 0 0 0 0 0 0 2 − 0.1760000000000000D +03 0.3840000000000000D +03 cos 0 0 0 0 0 0 0 0 2 0 0 2 cos 0 0 0 0 0 0 0 0 4 0 0 2 − 0.8000000000000000D +02 cos 0 0 0 0 0 0 0 1 0 4 0 2 − 0.7200000000000000D +02 cos 0 0 0 0 0 0 2 0 0 0 0 2 − 0.5400000000000000D +02 cos 0 0 0 0 0 0 2 0 2 0 0 2 0.1560000000000000D +03 cos 0 0 0 0 0 0 2 0 4 0 0 2 0.1000000000000000D +02 cos 0 0 4 0 0 0 0 0 0 4 0 2 − 0.8000000000000000D +02 cos 0 0 4 0 0 0 0 1 0 4 0 2 − 0.7200000000000000D +02 446 ANNALS NEW YORK ACADEMY OF SCIENCES TABLE 7/continued. f g h ξ a n e η c s (f − l) nδ coefficient 0.1000000000000000D +02 cos 0 0 4 0 0 0 2 0 0 4 0 2 cos 0 2 0 0 0 0 2 0 0 2 0 2 − 0.3600000000000000D +02 cos 0 2 0 0 0 0 2 0 2 2 0 2 0.6000000000000000D +02 cos 0 2 4 0 0 0 2 0 0 2 0 2 0.3000000000000000D +02 cos 0 2 4 0 0 0 2 0 1 2 0 2 0.6000000000000000D +02 cos 0 2 4 0 0 0 2 0 2 2 0 2 0.3000000000000000D +02 cos 0 2 −4 0 0 0 2 0 0 2 0 2 0.3000000000000000D +02 cos 0 2 −4 0 0 0 2 0 1 2 0 2 − 0.6000000000000000D +02 cos 0 2 −4 0 0 0 2 0 2 2 0 2 0.3000000000000000D +02 n WB 1 = --------- [ 12s 2 ( f – l ) cos ( 2h ) + 6es 2 sin ( f + 2h ) 8η 3 + 6es 2 sin ( f – 2h ) + 3e ( 1 + c ) 2 sin ( f + 2g + 2h ) + 3e ( 1 – c ) 2 sin ( f + 2g – 2h ) + 3 ( 1 + c ) 2 sin ( 2 f + 2g + 2h ) (23) + 3 ( 1 – c ) 2 sin ( 2 f + 2g – 2h ) + e ( 1 + c ) 2 sin ( 3 f + 2g + 2h ) + e ( 1 – c ) 2 sin ( 3 f + 2g – 2h ) ]. Many symmetries between the pairs ( H 1( 0 ) and W1) and ( HB 1( 0 ) and WB1) give confidence in the new first order results in δ. Once more, for the case δ = 0, we retrieve, of course, the known results for the isolated effect of J2: the expression for εW1 is the same as that given by (− S1), defined in Equation (15) of Reference 9. The result of computing the second-order averaged Hamiltonian is given in TABLE 7, where the series has been multiplied by a factor F = (3n2/64a2η7)−1 in order to make the result as readable as possible. In fact, the original series contained about 341 terms, that reduce to 44 terms after simplification (see next section). Again, we can further rewrite these results in full, by defining the second-order averaged Hamiltonian as ε 2 H 0( 2 ) + εδHM 0( 2 ) + δ 2 HB 0( 2 ) , (24) H 0( 2 ) × F = [ 5 ( s 4 – 8c 4 ) – 4η ( 1 – 3c 2 ) 2 – η 2 ( 5s 4 – 8c 2 ) – 2 ( 1 – 15c 2 ) cos ( 2g ) ], (25) with and HM 0( 2 ) × F = 32s 2 ( 5c 2 – 4 ) cos ( 2h ) + 48ηs 2 ( 3c 2 – 1 ) cos ( 2h ) – 4e 2 s 2 ( 5c 2 + 13 ) cos ( 2h ) – 4e 2 ( cos ( 2g + 2h ) + cos ( 2g – 2h ) ) + 52e 2 c ( cos ( 2g + 2h ) – cos ( 2g – 2h ) ) + 56e 2 c 2 ( cos ( 2g + 2h ) + cos ( 2g – 2h ) ) – 60e 2 c 3 ( cos ( 2g + 2h ) – cos ( 2g – 2h ) ) – 60e 2 c 4 ( cos ( 2g + 2h ) + cos ( 2g – 2h ) ), (26) DE SAEDELEER & HENRARD: AN ARTIFICIAL SATELLITE OF THE MOON 447 and also HB 0( 2 ) × F = 8 ( – 22 + 48c 2 – 10c 4 ) + 2e 2 ( – 27 + 78c 2 + 5c 4 ) – 72ηs 4 ( cos ( 4h ) + 1 ) + 10 ( e 2 – 8 )s 4 cos ( 4h ) – 12e 2 s 2 ( 3 – 5c 2 ) cos ( 2g ) (27) + 30e 2 s 2 ( 1 + c ) 2 cos ( 2g + 4h ) + 30e 2 s 2 ( 1 – c ) 2 cos ( 2g – 4h ). Again, symmetries (especially by the way of the tracer η) give high confidence in the new second-order averaged results in δ. Obviously, in the case δ = 0 we retrieve again the known results for the secondorder effect of J2 alone. The expression for ε2 H 0( 2 ) is the same as that given by ( – 2F 2* ) defined in Equation (29) of Reference 9. From the second-order averaged Hamiltonian, we could now derive the second-order equations of motion, but this is not our purpose here. We do not display the results for the second-order generator W 2* , because it is quite long (1,395 terms at the moment, coming from 2,711 terms originally). Further simplification needs to be done. For the part in ε2, containing only 241 terms, it looks very similar to other results,13 but the exact correspondence needs to be checked more closely in future work. Note that the generators may sometimes differ by a constant factor, independent of l. SIZE OF THE COMPUTATIONS AND TECHNICAL DETAILS A critical issue in making such computations is the number of terms appearing in the series (typically thousands of terms here). Consider, for example, computation of the first Poisson parenthesis of (16), namely, ( εH 1( 0 ) + δHB 1( 0 ) + εH 0( 1 ) + δHB 0( 1 ) ; εW 1 + δWB 1 ). (28) Expression (28) contains 2,302 terms that need to be integrated. The highest degree in f is (ξ6,7f), which finally gives up to 10f when expanded. This order of magnitude was expected when looking at (10) where we see factors (ξ3, 2f), at (22) where we see factors (ξ0, 3f), and at TABLE 2 where we see factors (ξ2, f). The last column of TABLE 8 gives the number of terms generated by the expansion of each power of ξ in the computation of (28). The other two columns give the same quantity, but for computations considering only ε or δ. The nonlinear increase in the number of terms stems from the coupling between ε and δ. Of course, we have to be careful to keep the number of terms at a reasonable level. Therefore, as said previously, we are concerned to simplify the expressions as far as possible. The redundancy of the set (8) implies, indeed, that we have several equivalent forms for expressions, with two main consequences: that some simplifications do not always appear at first glance and that comparison of results is not always direct. To avoid these drawbacks, we can choose to express the results in terms of the preferred variables (e, c) rather than (η, s), by replacing each occurrence of η2 by (1 − e2) and similarly s2 by (1 − c2), performing a substitute operation. For example, 448 ANNALS NEW YORK ACADEMY OF SCIENCES TABLE 8. Number of terms appearing in the second-order computation of (28) C22 J2 J2 and C22 ξ6 141 747 1,563 ξ5 177 1,060 2,142 ξ4 185 963 2,039 ξ3 78 518 1,051 ξ2 12 45 108 Total less than 1,000 less than 4,000 less than 7,000 in the computation of (24), the simple fact of multiplying the initial rough series by η9 (the common denominator) and then replacing η2 by (1 − e2) reduces the number of terms by a factor of about three (it decreases from 315 to 107). Among these simplifications are series of the kind e2 e2 e4 ------ – ------ + -----η7 η9 η9 that, in fact, cancel to zero. After these simplifications are made, we can obtain more compact forms by going in the opposite direction, using a collect operation; this reduces the number of terms from 107 to 44, which seems to be the most compact and simplified form. Recall that a further simplification needs to be done for the second-order generator, W 2* , because it is still quite long (1,395 terms for the time being, from 2,711 originally). Now a few words about how to integrate each term of the series. In fact, we do the integration in successive steps, written as separate subroutines. Terms with powers at least equal to two in ξ are integrated first. Among these, those containing a factor (f − l) receive a special treatment. The factor (f − l) is indeed known to play an important role in the problem of an artificial satellite.14 Moreover, it is worth performing the sequence of steps in the correct order; for example, in the computation of W 2* , some terms need not be integrated if we have previously integrated other terms by parts, since they then cancel each other.15 The lowest degree terms to integrate while computing (28) are in ξ0 (they are not listed in TABLE 8), and they present a special integration case (e.g., ∫2π ξ 0 dl ), since they 0 cannot be divided by ξ2η as simply as before. Furthermore, numerical accuracy considerations are also required; for example, one of the 1,395 terms appearing in the computation of W 2* is of order 10−14 (from the fact that we are working in double precision 10−16 with factors of order up to 102), and need to be canceled out. CONCLUSIONS We have second-order results for the combined effect of J2 and C22. One future task is to validate these results theoretically by numeric integration of the effect of J2 and C22. An experimental validation could also follow, by a comparison with data for a real lunar satellite. DE SAEDELEER & HENRARD: AN ARTIFICIAL SATELLITE OF THE MOON 449 We aim also to finish our validation of the computation of the second-order generator W 2* (by comparison with other results13), then to eliminate the next fastest angle (g or h) in order to obtain the secular motion. Subsequently, we would be ready to add other perturbations (Earth, Sun, other lunar harmonics like C31, librations, etc.). For this we need to improve the software MM (e.g., in order to increase the allowed size for the number of variables and probably number of terms). REFERENCES 1. MULLER, P.M. & W.L. S JOGREN. 1968. Mascons: lunar mass concentrations. Science 161: 680–684. 2. BILLS, B.G. & A.J. F ERRARI. 1980. A harmonic analysis of lunar gravity. J. Geophys. Res. 85: 1013–1025. 3. LEMOINE, F.G., D.E. S MITH, M.T. ZUBER, et al. 1997. A 70th degree lunar gravity model (GLGM-2) from Clementine and other tracking data. J. Geophys. Res. 102: 16,339–16,359. 4. KONOPLIV, A.S., S.W. A SMAR, E. CARRANZA, et al. 2001. Recent gravity models as a result of the lunar prospector mission. Icarus 150: 1–18. 5. FLOBERGHAGEN, R., R. NOOMEN, P.N.A.M V ISSER & G.D. RACCA. 1996. Global lunar gravity recovery from satellite-to-satellite tracking. Planet Space Sci. 44(10): 1081– 1097. 6. KAULA, W.M. 1966. Theory of Satellite Geodesy. Blaisdell Publishing Company. 7. COOK, A. 1988. The Motion of the Moon. IOP Publishing Ltd. 8. DEPRIT, A. 1969. Canonical transformations depending on a small parameter. Celest. Mech. Dyn. Astr. 1: 12–30. 9. BROUWER, D. 1959. Solution of the problem of artificial satellite theory without air drag. Astron. J. 64: 378–397. 10. SZEBEHELY, V.G. 1989. Adventures in Celestial Mechanics. University of Texas Press. 11. ROY, A.E. 1968. The theory of the motion of an artificial lunar satellite I. Development of the disturbing function. Icarus 9: 82–132. 12. JUPP, A.H. 1988. The critical inclination problem: 30 years of progress. Celest. Mech. Dyn. Astr. 43: 127–138. 13. KOZAI, Y. 1962. Second-order solution of artificial satellite theory without air drag. Astron. J. 67: 446–461. 14. METRIS, G. 1991. Mean values of particular functions in the elliptic motion. Celest. Mech. Dyn. Astr. 52: 79–84. 15. AKSNES, K. 1971. A Note on “the main problem of satellite theory for small eccentricities, by A. Deprit and A. Rom, 1970”. Celest. Mech. Dyn. Astr. 4: 119–121. View publication stats