A&A 471, 1069-1075 (2007)
DOI: 10.1051/0004-6361:20077568

Long-term harmonic development of lunar ephemeris

S. M. Kudryavtsev

Sternberg Astronomical Institute of Moscow State University, Universitetsky Pr. 13, Moscow, 119992, Russia

Received 29 March 2007 / Accepted 28 May 2007

Abstract
Aims. The aim of this study is to develop new analytical series representing lunar coordinates to accuracy compatible with the accuracy of the modern numerical ephemeris of the Moon.
Methods. An improved method of spectral analysis of tabulated function is used to make harmonic development of the latest long-term numerical ephemeris of the Moon LE-406 which covers a six thousand-year interval. A feature of the method is that the development is made directly to Poisson series where both amplitudes and arguments of the series' terms are high-degree polynomials of time.
Results. The new analytical development includes 42 270 Poisson series' terms of minimal amplitude equivalent to 1 cm and is valid over 1500-2500. A simplified version of the development includes 7952 series' terms of minimal amplitude equivalent to 1 m and is valid over 3000BC-3000AD. Over 1500-2500 the maximum difference between lunar coordinates calculated by means of the new analytical series and numerical ephemeris LE-406 is 3.2 m in geocentric distance, 0 $.\!\!^{\prime\prime}$0056 in ecliptic longitude, and 0 $.\!\!^{\prime\prime}$0018 in ecliptic latitude. This is 9-70 times better than the accuracy of the latest analytical theory of lunar motion ELP/MPP02, and the number of terms in the new development is less than that in ELP/MPP02.

Key words: celestial mechanics - methods: analytical - Moon

1 Introduction

Lunar motion theory is a classical task of celestial mechanics. It was a research topic dealt with by Newton, Clairaut, d'Alembert, Euler, Lagrange, Laplace, Hansen, Adams, Delaunay, Andoyer, Newcomb and others. Also important is the famous analytical theory by Hill & Brown (Hill 1878a,b,c; Brown 1897, 1899, 1905, 1908); for a long time the relevant tables (Brown 1919) were the basis for calculation of lunar coordinates. Recent progress in computer technology and new requirements for accuracy of lunar ephemeris stimulated further studies to improve the Hill-Brown method and series (Eckert et al. 1954; Schmidt 1980) as well as on the development of new analytical theories of lunar motion (Deprit et al. 1971a,b; Henrard 1979, 1980, 1981; Chapront-Touzé & Chapront 1983, 1988; Chapront & Chapront-Touzé 1997; Bidart 2001; Chapront et al. 2002; Chapront & Francou 2003).

At the same time, numerical ephemerides of the Moon have been successfully developed; among them the most recent and accurate are the ephemerides of DE/LE-series done at the Jet Propulsion Laboratory (JPL), USA (Standish & Williams 1981; Standish et al. 1995; Standish 1998, 2003a,b, 2006), the ephemerides of the EPM-series done at Russian Institute of Applied Astronomy (Krasinsky 2002; Pitjeva 2001, 2003, 2005) and the ephemerides of INPOP-series done at IMCCE-Observatoire de Paris, France (Fienga et al. 2006). Currently, the accuracy of the above-mentioned numerical ephemerides of the Moon is better than that of available analytical theories of lunar motion, but an important advantage of the latter is their compactness and computer platform independence which is important in many practical applications. In particular, JPL-based numerical ephemerides of the Moon and planets have been recently replaced with analytical motion theories of these bodies within all key elements of ground software systems used for Hubble Space Telescope (HST) mission support (McCuttcheon 2003). For this purpose the HST personnel chose the analytical theory of lunar motion ELP2000-82B (Chapront-Touzé & Chapront 1983) and the "planetary series 1996'' (Chapront & Francou 1996) representing planetary ephemerides by Poisson series. According to HST tests, the selected analytical motion theories reproduce the modern numerical planetary/lunar ephemerides DE/LE-405/406 (Standish 1998) to an accuracy of at least 0 $.\!\!^{\prime\prime}$008 for all major planets (which meets all requirements for the mission support), and to an accuracy of 0 $.\!\!^{\prime\prime}$5 for the Moon. The main reason of the relatively low accuracy obtained for the Moon is that ELP2000-82B theory is adjusted to an old version of numerical lunar ephemeris, LE-200 (Standish & Williams 1981). (Recently Chapront & Francou (2003) developed a new analytical theory of lunar motion, ELP/MPP02, which is adjusted to LE-405/406, but coefficients of the theory are not published yet).

Along with construction of purely analytical (or semi-analytical) theories of lunar motion and development of purely numerical lunar ephemeris a combined approach can be used. This is a spectral analysis of values for lunar coordinates pre-calculated with a small sampling step on the basis of the latest long-term numerical ephemeris of the Moon. The form of the resulting series is very similar to that given by the modern analytical theories of lunar motion and keeps all the advantages of the latter, and the series accuracy proves to be compatible with the accuracy of the source numerical ephemeris. Note that a similar approach was used by Chapront (1995, 2000) in his improvement of analytical planetary theories. However, a typical disadvantage of any spectral method is the well-known problem of "close frequencies'' which can only be resolved by increasing the time interval over which the analysis is done. We thus evolved a new modification of the spectral analysis method which allowed us to perform harmonic development of the latest long-term numerical ephemeris of the Moon LE-406 over six thousand years (the complete time interval covered by LE-406). The feature of the method is that the development is directly made to Poisson series where both amplitudes and arguments of the series' terms are high-degree polynomials of time as opposed to the classical Fourier analysis where the terms' amplitudes and frequencies are constants. This approach leads to an essential improvement in accuracy of harmonic development of a function tabulated over a long-term interval and to reduction of the series' length. The details of the method are given in the next section (see also Kudryavtsev 2004).

2 Spectral analysis to Poisson series

Let f(t) be an arbitrary function tabulated by its numerical values over an interval of time [-T,T] with a small sampling step.

Over the same interval we will build an analytical representation of the function by a finite h-order Poisson series of the following form

 

$$f(t)\approx\sum_{k=1}^{N}\Bigl[\left(A^c_{k0}+A^c_{k1}t+\cdots+A^c_{kh}t^h % \right)\cos{\omega_k(t)}$$

$$+\left(A^s_{k0}+A^s_{k1}t+\cdots+A^s_{kh}t^h \right)\sin{\omega_k(t)}\Bigr]$$ (1)

where $A^c_{k0},\dots,A^s_{kh}$ are constants and $\omega_k(t)$ are some pre-defined arguments which are assumed to be q-degree polynomials of time t

$$\omega_k(t)=\nu_kt+\nu_{k2}t^2+\cdots+\nu_{kq}t^q.$$ (2)

For this we find the projections of f(t) on a basis generated by functions

 

$$\vec{c}_{kl}(t)\equiv{t}^l\cos{\omega_k(t)}, \quad \vec{s}_{kl}(t)\equiv{t}^l\sin{\omega_k(t)} \quad$$ (3)

$$(k=1,2,\cdots,N; \: l=0,1,\cdots,h)$$

through numerical computation of the following scalar products

$$A^c_{kl}=\langle f,\vec{c}_{kl}\rangle \equiv \frac{1}{2T}\int_{-T}^Tf(t)t^l \cos{\omega_k(t)}\chi(t)~{\rm d}t,$$ (4)

$$A^s_{kl}=\langle f,\vec{s}_{kl}\rangle \equiv \frac{1}{2T}\int_{-T}^Tf(t)t^l \sin{\omega_k(t)}\chi(t)~{\rm d}t$$ (5)

by using the definition

$$\langle f,g\rangle \equiv \frac{1}{2T}\int_{-T}^Tf(t) \bar g(t)\chi(t)~{\rm d}t$$ (6)

where $\bar g$ is the complex conjugate to the g function and $\chi(t)=1+\cos\frac{\pi}{T}t$ is the Hanning filter chosen as the weight function following suggestions by Laskar et al. (1992), Laskar (1993). The proper choice of arguments $\omega_k(t)$ depends on the specific task (e.g. they can be multipliers of Delaunay variables and/or planetary mean longitudes, etc.).

However, the basis functions $\vec{c}_{k_1l_1},~\vec{s}_{k_1l_1},~ \vec{c}_{k_2l_2},~\vec{s}_{k_2l_2},\dots~$ are not usually orthogonal. We thus have to perform an orthogonalization process over the expansion coefficients to improve the quality of representation (1) and avoid superfluous terms. For this procedure we used the algorithm developed by Sidlichovský & Nesvorný (1997). Expressions (7)-(14) present a minimal algorithm which we have generalized as indicated below.

Let f(t) be a tabulated complex function and let $\{\vec{e}_i \}_{i=1,2,\cdots,M}$ be a set of M basis functions [in our study equal to the complete set of $\vec{c}_{kl}(t),~\vec{s}_{kl}(t)$so that $M=2\times N \times{(h+1)}$]. The function f(t) is developed in the basis $\{\vec{e}\}$ as

$$f(t)=\sum^M_{i=1}A_i^{(M)}\vec{e}_i+f_M(t)$$ (7)

where A_i^(M) is a coefficient at $\vec{e}_i$ after expanding f(t) over M basis functions, and f_M(t), the difference between the original function and its approximation by M terms, proves to be minimal. Let us define the projections $F_i\equiv \langle f_{i-1},\vec{e}_i\rangle$ and $Q_{ij}\equiv \langle \vec{e}_i,\vec{e}_j\rangle$. The original algorithm by Sidlichovský & Nesvorný (1997) employs a certain normalized basis $\{\vec{e}\}$ where the latter scalar product $\left(Q_{ij}\right)$is always a real-valued function. We expand their result to the case of an arbitrary non-normalized basis $\{\vec{e}\}$ where Q_ij can take complex values as well.

Thus coefficients A_i^(M) are iteratively calculated as follows. At the first step

$$\alpha_{11}=\frac{1}{\sqrt{Q_{11}}}, \quad A_1^{(1)}=\alpha^2_{11}F_1, \quad f_1(t)=f_0(t)-A_1^{(1)}\vec{e}_1$$ (8)

where $f_0(t)\equiv f(t)$ and $\alpha_{ij}$ are hereafter some calculated complex constants.

At the $m{\rm th}$ step, for every $j=1,2,\cdots,m-1$ we compute the following complex coefficients:

$$B_j^{(m)}=-\sum^j_{s=1}\bar \alpha_{js}Q_{ms},$$ (9)

$$\alpha_{mm}=\left( {Q_{mm}-\sum^{m-1}_{s=1}\bar B_s^{(m)}B_s^{(m)}}\right) ^{-\frac{1}{2}}$$ (10)

(by construction the coefficient $\alpha_{mm}$ takes on a real value for every m),

$$\alpha_{mj}=\alpha_{mm}\sum^{m-1}_{s=j}B_s^{(m)}\alpha_{sj},$$ (11)

$$A_m^{(m)}=\alpha^2_{mm}F_m, \quad$$ (12)

$$A_j^{(m)}=A_j^{(m-1)}+\alpha_{mm}\alpha_{mj}F_m,$$ (13)

$$f_m(t)=f_{m-1}(t)-\alpha_{mm}F_m\sum^m_{i=1}\alpha_{mi}\vec{e}_i$$ (14)

where $\bar \alpha_{js}$ and $\bar B_s^{(m)}$ are complex conjugate values of the relevant quantities.

For the selected basis (3) the projections $F_i\equiv~ \langle f_{i-1},\vec{e}_i\rangle$ are numerically calculated according to expressions (4)-(5). The values for scalar products of the basis functions $Q_{ij}\equiv \langle \vec{e}_i,\vec{e}_j\rangle$ can be found analytically through the following steps.

Step 1. As far as trigonometric functions can be represented in exponential form, we shall further deal with definite integrals of the form

$$I_n(\nu)\equiv \frac{1}{2T}\int_{-T}^Tt^n {\rm e}^{{\rm i}\n... ...frac{\pi}{T}t}\right){\rm d}t=I_n^a(\nu)+I_n^b(\nu)+I_n^c(\nu)$$ (15)

where ${\rm i}\equiv \sqrt{-1}$ and

$$I_n^a(\nu)\equiv \frac{1}{2T}\int_{-T}^Tt^n{\rm e}^{{\rm i}\nu t}{\rm d}t,$$ (16)

$$I_n^b(\nu)\equiv \frac{1}{4T}\int_{-T}^Tt^n{\rm e}^{{\rm i}(\nu+\frac{\pi}{T})t}{\rm d}t,$$ (17)

$$I_n^c(\nu)\equiv \frac{1}{4T}\int_{-T}^Tt^n{\rm e}^{{\rm i}(\nu-\frac{\pi}{T})t}{\rm d}t.$$ (18)

It is easy to find

$$I_n^a(0)=2I_n^b\left(-\frac{\pi}{T}\right)=2I_n^c\left(\frac{... ...+1} & \mbox{ if}&n \quad \mbox{is even.}\\ \end{array}\right.$$ (19)

Otherwise, if n=0

$$I_0^a(\nu)=\frac{\sin(\nu T)}{\nu T},$$ (20)

$$I_0^b(\nu)=-\frac{\sin(\nu T)}{2(\nu T+\pi)},$$ (21)

$$I_0^c(\nu)=-\frac{\sin(\nu T)}{2(\nu T-\pi)}\cdot$$ (22)

If $n\ge 1$ we calculate the integrals iteratively

$$I_n^a(\nu)=\frac{\rm i}{\nu T}\left(TnI^a_{n-1}-T^n\psi(\nu)\right),$$ (23)

$$I_n^b(\nu)=\frac{\rm i}{\nu T+\pi}\left(TnI^b_{n-1}+\frac{1}{2}T^n\psi(\nu)\right),$$ (24)

$$I_n^c(\nu)=\frac{\rm i}{\nu T-\pi}\left(TnI^c_{n-1}+\frac{1}{2}T^n\psi(\nu)\right)$$ (25)

where

$$\psi(\nu)= \left\{\begin{array}{rcl} \cos{\nu T} & \mbox{ if}... ... T} & \mbox{ if}&n \quad \mbox{is even.}\\ \end{array}\right.$$ (26)

Step 2. We partially expand the exponential function of the argument (2) to a power series of t by assuming smallness of the second and further items in the right-hand side of expression (2) with respect to the first term. (In particular, it is true for Delaunay variables, planetary mean longitudes and many other expansion functions of celestial mechanics.) This task can easily be performed by means of a computer algebra system (we have used MAPLE V program package by Waterloo Maple Software). The result is as follows (where the maximal degree of the polynomial argument, q, has been restricted to a value of 4, although similar expansions can be obtained for any value of degree q)

$${\rm e}^{{\rm i}(\nu t+\nu_2t^2\!+\!\nu_3t^3\!+\!\nu_4t^4)}\... ...gl({\rm i}\nu_4-\frac{1}{2}\nu_2^2\Bigr)t^4 +\cdots\Bigr]\cdot$$ (27)

In our study we obtain the expansion (27) up to the terms proportional to t²⁴. We then define a new integral function of a polynomial argument $\omega (t)$ as

$$J_n(\omega)\equiv \frac{1}{2T} \int_{-T}^{T}t^n {\rm e}^{{\r... ...nu_3t^3+\nu_4t^4)}\left(1+\cos{\frac{\pi}{T}t}\right){\rm d}t.$$ (28)

By combining expressions (15) and (27) we can write

$$J_n(\omega) = I_n(\nu) + {\rm i}\nu_2I_{n+2}(\nu) + {\rm i}\nu_3I_{n+3}(\nu)$$

$$+ \Bigl({\rm i}\nu_4-\frac{1}{2}\nu_2^2\Bigr)I_{n+4}(\nu) + \cdots$$ (29)

Step 3. We expand trigonometric functions to exponential ones, and then the required scalar products can be calculated as

$$\langle \vec{c}_{k_1l_1},\vec{s}_{k_2l_2}\rangle \frac{1}{2T}\int_{-T}^{T}t^{l_1+l_2} \Bigl(\frac{{\rm e}^{{\rm i}... ...l(\frac{{\rm e}^{{\rm i}\omega k_2}-{\rm e}^{-{\rm i}\omega k_2}}{2\rm i}\Bigr)$$ =

$$\times \left(1+\cos{\frac{\pi}{T}t}\right){\rm d}t$$

$$\frac{1}{4\rm i}\bigl[-J_{l_1+l_2}( \omega_{k_1}-\omega_{k_2}) -J_{l_1+l_2}(-\omega_{k_1}-\omega_{k_2})$$ =

$$+J_{l_1+l_2}( \omega_{k_1}+\omega_{k_2}) +J_{l_1+l_2}(-\omega_{k_1}+\omega_{k_2})\bigr],$$ (30)

and analogously

$$\langle \vec{s}_{k_1l_1},\vec{c}_{k_2l_2}\rangle = \frac{1}{4\rm ... ..._{l_1+l_2}( \omega_{k_1}-\omega_{k_2}) -J_{l_1+l_2}(-\omega_{k_1}-\omega_{k_2})$$

$$+J_{l_1+l_2}( \omega_{k_1}+\omega_{k_2}) -J_{l_1+l_2}(-\omega_{k_1}+\omega_{k_2})\bigr],$$ (31)

$$\langle \vec{s}_{k_1l_1},\vec{s}_{k_2l_2}\rangle= \frac{1}{4}\bigl[J_{l_1+l_2}( \omega_{k_1}-\omega_{k_2}) -J_{l_1+l_2}(-\omega_{k_1}-\omega_{k_2})$$

$$-J_{l_1+l_2}( \omega_{k_1}+\omega_{k_2}) +J_{l_1+l_2}(-\omega_{k_1}+\omega_{k_2})\bigr],$$ (32)

 

$$\langle \vec{c}_{k_1l_1},\vec{c}_{k_2l_2}\rangle= \frac{1}{4}\bigl[J_{l_1+l_2}( \omega_{k_1}-\omega_{k_2}) +J_{l_1+l_2}(-\omega_{k_1}-\omega_{k_2})$$

$$+J_{l_1+l_2}( \omega_{k_1}+\omega_{k_2}) +J_{l_1+l_2}(-\omega_{k_1}+\omega_{k_2})\bigr].$$ (33)

This is a modification of the spectral analysis method allowing expansion of a tabulated function to Poisson series with the arguments being high-degree polynomials of time.

3 Expansion form of lunar coordinates

The formalism described in the previous section has been applied to accurate harmonic development of the numerical lunar ephemeris LE-405/406. As a set of variables describing the position of the Moon in space we chose spherical coordinates of its centre: r (geocentric distance), V (ecliptic longitude reckoned along the moving ecliptic from the mean equinox of epoch) and U (ecliptic latitude reckoned from the moving ecliptic). The same variables are used in all modern ELP-theories of lunar motion (although the ELP-theories use a different origin of the ecliptic longitude V).

On the base of LE-406 numerical lunar ephemeris, we calculated daily values for spherical coordinates of the Moon r, V, U over 3000BC-3000AD (the complete time interval covered by LE-406). We then performed a spectral analysis of the tabulated values using the method presented in the previous section. As a result, the geocentric spherical coordinates of the Moon r, V, U are represented by Poisson series of the following form

 

$$\sum_{k=1}^{N_r} \bigl[ A_{k0}^{(r)}\cos\bigl(\omega_k^{(r)}(t)+\... ...^{(r)}\bigr) +A_{k1}^{(r)}t\cos\bigl(\omega_k^{(r)}(t)+\varphi_{k1}^{(r)}\bigr)$$                         r(t) =

$$+A_{k2}^{(r)}t^2\cos\bigl(\omega_k^{(r)}(t)+\varphi_{k2}^{(r)}\bigr)\bigr],$$ (34)

 

$$\bar V(t)+\sum_{k=1}^{N_V} \bigl[ A_{k0}^{(V)}\sin\bigl(\omega_k^{(V)}(t)+\varphi_{k0}^{(V)}\bigr)$$                         V(t) =

$$+A_{k1}^{(V)}t\sin\bigl(\omega_k^{(V)}(t)+\varphi_{k1}^{(V)}\bigr)$$

$$\quad +A_{k2}^{(V)}t^2\sin\bigl(\omega_k^{(V)}(t)+\varphi_{k2}^{(V)}\bigr)\bigr],$$ (35)

 

$$U(t)=\sum_{k=1}^{N_U} \bigl[ A_{k0}^{(U)}\sin\bigl(\omega_k^{(U)}(t)+\varphi_{k0}^{(U)}\bigr)$$

$$+A_{k1}^{(U)}t\sin\bigl(\omega_k^{(U)}(t)+\varphi_{k1}^{(U)}\bigr)$$

$$+A_{k2}^{(U)}t^2\sin\bigl(\omega_k^{(U)}(t)+\varphi_{k2}^{(U)}\bigr)\bigr]$$ (36)

where

 

$$\bar V(t) 218\hbox{$.\!\!^\circ$ }31664563+17325643723\hbox{$.\!\!^{\prime\prime}$ }0470t-527\hbox{$.\!\!^{\prime\prime}$ }90t^2$$ =

$$+6\hbox{$.\!\!^{\prime\prime}$ }665t^3-0\hbox{$.\!\!^{\prime\prime}$ }5522t^4$$ (37)

is the mean longitude of the Moon referred to the moving ecliptic and mean equinox of epoch (Simon et al. 1994); t is hereafter TDB measured in thousands of Julian years (365250d) from J2000.0 (JD2451545.0); $\omega_k(t)$ are time polynomial arguments, and A_ki, $\varphi_{ki}$ are, respectively, the terms' amplitudes and some additional constants (i=0,1,2, i.e. the maximum order of the Poisson series in our development is chosen equal to 2). The algorithm of spectral analysis of a tabulated function to Poisson series suggested in the previous section works for any assigned order of the series, but a limitation here is the time interval over which the function is tabulated (numerical tests prove that expanding a function tabulated over a relatively small interval of time to high-order Poisson series leads to strong correlation between amplitudes of different Poisson terms of the same frequency and does not improve the quality of the development). It is difficult to estimate the right order of the series in advance, so in the present study we restricted the maximum order of Poisson series for the lunar coordinates to 2 following the similar form of the ELP-solutions. Indeed, a further analysis of differences between the lunar coordinates given by the final expansions (34)-(37) and by the complete LE-405/406 numerical ephemeris did not reveal any systematic residuals which would indicate a higher order of Poisson series is necessary for representing lunar coordinates over the considered six thousand-year interval.

The arguments $\omega_k(t)$ in expansions (34)-(36) are defined as follows. First, from Simon et al. (1994) we took fourth-degree time polynomial expressions for the mean longitude of the ascending node of the Moon $\Omega$ (referred to the mean ecliptic and equinox of J2000.0), for Delaunay variables D, l', l, F (mean elongation of the Moon from the Sun, mean anomaly of the Sun, mean anomaly of the Moon, and mean longitude of the Moon subtracted by $\Omega$, respectively), for mean longitudes of the eight major planets (referred to the mean ecliptic and equinox of J2000.0), and for the general precession in longitude p_A (based on Williams et al. 1991)

 

$$\Omega 125\hbox{$.\!\!^\circ$ }04455501 - 69679193\hbox{$.\!\!^{\prime\prime}$ }631t + 636\hbox{$.\!\!^{\prime\prime}$ }02t^2$$ =

$$+ 7\hbox{$.\!\!^{\prime\prime}$ }625t^3 - 0\hbox{$.\!\!^{\prime\prime}$ }3586t^4,$$ (38)

$$297\hbox{$.\!\!^\circ$ }85019547 + 16029616012\hbox{$.\!\!^{\prime\prime}$ }090t - 637\hbox{$.\!\!^{\prime\prime}$ }06t^2$$                             D =

$$+ 6\hbox{$.\!\!^{\prime\prime}$ }593t^3 - 0\hbox{$.\!\!^{\prime\prime}$ }3169t^4,$$ (39)

$$357\hbox{$.\!\!^\circ$ }52910918 + 1295965810\hbox{$.\!\!^{\prime\prime}$ }481t - 55\hbox{$.\!\!^{\prime\prime}$ }32t^2$$                             l' =

$$+ 0\hbox{$.\!\!^{\prime\prime}$ }136t^3 - 0\hbox{$.\!\!^{\prime\prime}$ }1149t^4,$$ (40)

$$134\hbox{$.\!\!^\circ$ }96340251 + 17179159232\hbox{$.\!\!^{\prime\prime}$ }178t + 3187\hbox{$.\!\!^{\prime\prime}$ }92t^2$$                             l =

$$+ 51\hbox{$.\!\!^{\prime\prime}$ }635t^3 - 2\hbox{$.\!\!^{\prime\prime}$ }4470t^4,$$ (41)

$$93\hbox{$.\!\!^\circ$ }27209062 + 17395272628\hbox{$.\!\!^{\prime\prime}$ }478t - 1275\hbox{$.\!\!^{\prime\prime}$ }12t^2$$                             F =

$$- 1\hbox{$.\!\!^{\prime\prime}$ }037t^3 + 0\hbox{$.\!\!^{\prime\prime}$ }0417t^4,$$ (42)

$$\lambda_{Me} 252\hbox{$.\!\!^\circ$ }25090552 + 5381016286\hbox{$.\!\!^{\prime\prime}$ }88982t$$ =

$$- 1\hbox{$.\!\!^{\prime\prime}$ }92789t^2+ 0\hbox{$.\!\!^{\prime\prime}$ }00639t^3,$$ (43)

$$\lambda_{Ve} 181\hbox{$.\!\!^\circ$ }97980085 + 2106641364\hbox{$.\!\!^{\prime\prime}$ }33548t$$ =

$$+ 0\hbox{$.\!\!^{\prime\prime}$ }59381t^2 - 0\hbox{$.\!\!^{\prime\prime}$ }00627t^3,$$ (44)

$$\lambda_{Ea} 100\hbox{$.\!\!^\circ$ }46645683 + 1295977422\hbox{$.\!\!^{\prime\prime}$ }83429t$$ =

$$- 2\hbox{$.\!\!^{\prime\prime}$ }04411t^2- 0\hbox{$.\!\!^{\prime\prime}$ }00523t^3,$$ (45)

$$\lambda_{Ma} 355\hbox{$.\!\!^\circ$ }43299958 + 689050774\hbox{$.\!\!^{\prime\prime}$ }93988t$$ =

$$+ 0\hbox{$.\!\!^{\prime\prime}$ }94264t^2- 0\hbox{$.\!\!^{\prime\prime}$ }01043t^3,$$ (46)

$$\lambda_{Ju} 34\hbox{$.\!\!^\circ$ }35151874 + 109256603\hbox{$.\!\!^{\prime\prime}$ }77991t - 30\hbox{$.\!\!^{\prime\prime}$ }60378t^2$$ =

$$+ 0\hbox{$.\!\!^{\prime\prime}$ }05706t^3 + 0\hbox{$.\!\!^{\prime\prime}$ }04667t^4,$$ (47)

$$\lambda_{Sa} 50\hbox{$.\!\!^\circ$ }07744430 + 43996098\hbox{$.\!\!^{\prime\prime}$ }55732t + 75\hbox{$.\!\!^{\prime\prime}$ }61614t^2$$ =

$$- 0\hbox{$.\!\!^{\prime\prime}$ }16618t^3 - 0\hbox{$.\!\!^{\prime\prime}$ }11484t^4,$$ (48)

$$\lambda_{Ur} 314\hbox{$.\!\!^\circ$ }05500511 + 15424811\hbox{$.\!\!^{\prime\prime}$ }93933t - 1\hbox{$.\!\!^{\prime\prime}$ }75083t^2$$ =

$$+ 0\hbox{$.\!\!^{\prime\prime}$ }02156t^3,$$ (49)

$$\lambda_{Ne} 304\hbox{$.\!\!^\circ$ }34866548 + 7865503\hbox{$.\!\!^{\prime\prime}$ }20744t + 0\hbox{$.\!\!^{\prime\prime}$ }21103t^2$$ =

$$- 0\hbox{$.\!\!^{\prime\prime}$ }00895t^3,$$ (50)

 

$$p_A=50288\hbox{$.\!\!^{\prime\prime}$ }200t + 111\hbox{$.\!\!^{\p... ...\hbox{$.\!\!^{\prime\prime}$ }0773t^3 - 0\hbox{$.\!\!^{\prime\prime}$ }2353t^4.$$ (51)

Every argument $\omega_k(t)$ in expansions (34)-(36) is a linear combination of integer multipliers of expressions (38)-(51). Given that Delaunay variables are only defined by fourth-degree time polynomials, we have truncated the original expressions for $\lambda_{Me}, \cdots, \lambda_{Ne}$ and p_A to the same level of accuracy in our study. To select the set of arguments we preliminarily evaluated the spectrum of the tabulated numerical values of r, V, U at numerous combinations of integer multipliers of basic frequencies of expressions (38)-(51) by using the classical Fourier analysis (at this stage all the frequencies were cut to linear functions of time). Here we considered all combinations of the multipliers where the latter run within certain limits. The wider the limits are the more combinations can be taken into consideration, but based on our computer efficiency we ultimately chose the following maximum ranges:

• from -6 to 6 - for multipliers of Delaunay variables and of the mean longitude of the lunar ascending node;
• from -20 to 20 - for multipliers of the planetary mean longitudes;
• from -4 to 4 - for multipliers of the general precession in longitude.

To accelerate the evaluation process the Fast Fourier Transformation (FFT) of the data arrays at the frequencies specially defined by the FFT was made first. The approximate amplitude of a spectrum's wave at every considered combination of the basic frequencies was then found through interpolating the results of FFT. This method is inspired by studies by Laskar et al. (1992), Laskar (1993), where the authors used interpolating of the FFT results to find the fundamental frequencies of a dynamical system. Subsequently, we made improved spectral analysis of the original data arrays at all combinations of frequencies for which the preliminary (FFT) amplitude exceeded or was equal to a pre-set minimal level (equivalent to 1 cm in our study). This was done by using expressions (4)-(33) to account for the fourth-degree polynomial form of arguments $\omega_k(t)$. We emphasize that when making the improved spectral analysis we only use the time-dependent part of arguments $\omega_k(t)$ of form (2), but when calculating the lunar coordinates by means of expansions (34)-(37) the complete expressions (38)-(51) (incl. the free terms) should be used to obtain the arguments  $\omega_k(t)$. The constants  $\varphi_{ki}$ in expressions (34)-(36) are calculated so as to off-set the effect of those free terms as well as account for the additional phases resulting from transformation of the original series of form (1) to series with a single trigonometric function.

Transformation of lunar rectangular coordinates from the reference frame of LE-405/406 (defined by the mean geoequator and equinox of J2000.0) to the reference frame defined by the moving ecliptic and mean equinox of epoch has been done with use of the following precession quantities (Simon et al. 1994)

 

$$\theta_A 20042\hbox{$.\!\!^{\prime\prime}$ }0207t - 42\hbox{$.\!\!^{\prime\prime}$ }6566t^2 - 41\hbox{$.\!\!^{\prime\prime}$ }8238t^3$$ =

$$- 0\hbox{$.\!\!^{\prime\prime}$ }0731t^4 - 0\hbox{$.\!\!^{\prime\prime}$ }0127t^5 + 0\hbox{$.\!\!^{\prime\prime}$ }0004t^6,$$ (52)

$$\zeta_A 23060\hbox{$.\!\!^{\prime\prime}$ }9097t + 30\hbox{$.\!\!^{\prime\prime}$ }2226t^2 + 18\hbox{$.\!\!^{\prime\prime}$ }0183t^3$$ =

$$- 0\hbox{$.\!\!^{\prime\prime}$ }0583t^4 - 0\hbox{$.\!\!^{\prime\prime}$ }0285t^5 - 0\hbox{$.\!\!^{\prime\prime}$ }0002t^6,$$ (53)

$$23060\hbox{$.\!\!^{\prime\prime}$ }9097t + 109\hbox{$.\!\!^{\prime\prime}$ }5270t^2 + 18\hbox{$.\!\!^{\prime\prime}$ }2667t^3$$                            z_A =

$$- 0\hbox{$.\!\!^{\prime\prime}$ }2821t^4 - 0\hbox{$.\!\!^{\prime\prime}$ }0301t^5 - 0\hbox{$.\!\!^{\prime\prime}$ }0001t^6,$$ (54)

 

$$\varepsilon_A 23\hbox{$^\circ$ }26\hbox{$^\prime$ }21\hbox{$.\!\!^{\prime\prime... ...468\hbox{$.\!\!^{\prime\prime}$ }0927t - 0\hbox{$.\!\!^{\prime\prime}$ }0152t^2$$ =

$$+ 1\hbox{$.\!\!^{\prime\prime}$ }9989t^3 - 0\hbox{$.\!\!^{\prime\prime}$ }0051t^4 - 0\hbox{$.\!\!^{\prime\prime}$ }0025t^5.$$ (55)

Although such a transformation does not completely meet the present recommendations on precession quantities (Capitaine et al. 2003; McCarthy & Petit 2003) we prefer the formulae (52)-(55) because they are highly consistent with the other expressions (37)-(51) taken from the same source (Simon et al. 1994). In fact, the use of the precession quantities (52)-(55) does not decrease the accuracy of calculations, because one needs only to employ the same formulae (52)-(55) to transform lunar coordinates referred to the moving ecliptic and mean equinox of epoch (obtained by means of analytical expansions (34)-(37)) to the LE-405/406 reference plane.

The characteristics of the new harmonic development of lunar coordinates r, V, U are given in the following section.

4 Harmonic development of lunar ephemeris LE-405/406

There are two versions of the new Poisson series representing the geocentric ecliptic spherical coordinates of the Moon r, V, U. The complete solution, LEA-406a, includes 42270 terms of minimal amplitude equivalent to 1 cm and is valid over 1500-2500. The simplified solution, LEA-406b, includes 7952 terms of minimal amplitude equivalent to 1 m and is valid over 3000BC-3000AD.

Tables 1, 2 present characteristics of both solutions. For every coordinate we give the total number of non-zero trigonometric and Poisson terms in the relevant expression ((34), (35) or (36)). The value of minimum amplitude for Poisson terms corresponds to the maximum range of the terms' amplitude at the ends of the interval of time where the solution is valid. The value of maximum error stands for the maximum deviation between the lunar coordinates calculated by means of the new analytical development LEA-406 and lunar coordinates provided by the numerical ephemeris LE-405/406 at every 0.1 days within the considered interval of time. We made a comparison with the ephemeris LE-405 over the time intervals 1900-2100 and 1600-2200 and a comparison with the ephemeris LE-406 over the time intervals 1500-2500 and 3000BC-3000AD. (The ephemeris LE-405 is formally more accurate than LE-406, but covers a shorter time interval, 1600-2200; the ephemeris LE-406 is an extension of LE-405 over the interval of time 3000BC-3000AD.)

Table 1: Harmonic development of lunar ephemeris LE-405/406 over 1500-2500 (complete solution LEA-406a).

Table 2: Harmonic development of lunar ephemeris LE-405/406 over 3000BC-3000AD (simplified solution LEA-406b).

	Figure 1: Differences between lunar coordinates given by numerical ephemeris LE-406 and by analytical series LEA-406a over 1500-2500.
Open with DEXTER

	Figure 2: Dependence of LEA-406a accuracy on the number of series' terms.
Open with DEXTER

Figure 1 shows the differences between the lunar coordinates r, V, U calculated by means of the latest long-term numerical ephemeris LE-406 and by the complete harmonic development LEA-406a over 1500-2500.

Figure 2 demonstrates how the accuracy of calculation of lunar coordinates r, V, U over 1500-2500 depends on the number of terms taken from LEA-406a series.

The main terms of LEA-406a series are presented in Tables 3-5. The complete LEA-406 harmonic development of the lunar ephemeris LE-405/406 is given in Tables 6-11 only available in electronic form at the CDS.

5 Comparison of LEA-406 harmonic development with ELP/MPP02 analytical theory

LEA-406, the new harmonic development of lunar ephemeris, and ELP/MPP02 (Chapront & Francou 2003), the modern analytical theory of lunar motion, are both adjusted to the latest long-term numerical ephemeris of the Moon LE-405/406 (Standish 1998) which is recommended for use by the current IERS Conventions (McCarthy & Petit 2003). Both solutions are expansions of the geocentric ecliptic spherical coordinates of the Moon r, V, U to Poisson series where amplitudes of the series' terms are second-order polynomials of time, and their arguments are polynomials of the fourth degree of time.

For quantitative comparison of LEA-406 harmonic development with ELP/MPP02 lunar motion theory we give characteristics of the latter in Table 12 (following Chapront & Francou 2003).

When comparing the characteristics of the two analytical representation of lunar ephemeris given in Tables 1, 2 and in Table 12 one sees that the accuracy of the new harmonic development LEA-406 is better than that of ELP/MPP02 analytical theory over both short-term and long-term intervals. Over a one thousand-year interval centered at J2000.0 (1500-2500) the gain in accuracy is from a factor of 9 to a factor of 70 (depending on the coordinate). The total number of terms included to the new harmonic development of lunar ephemeris is 42 270 in its complete version (LEA-406a) and 7952 in its simplified version (LEA-406b) vs. 45 053 terms composing the analytical lunar motion theory ELP/MPP02.

Table 3: The main terms of LEA-406a harmonic development of lunar geocentric distance r.

Table 4: The main terms of LEA-406a harmonic development of lunar ecliptic longitude V.

Table 5: The main terms of LEA-406a harmonic development of lunar ecliptic latitude U.

Table 12: Expansion of lunar coordinates in ELP/MPP02 theory.

Acknowledgements

The author is grateful to Dr. Jean Chapront for his valuable comments on the manuscript. Research supported in part by the Russian Foundation for Basic Research under grant no. 05-02-16436.

References

• Bidart, P. 2001, A&A, 366, 351 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Bretagnon, P., & Francou, G. 1988, A&A, 202, 309 [NASA ADS]
• Brown, E. W. 1897, Memoirs R.A.S., 53 (In the text)
• Brown, E. W. 1899, Memoirs R.A.S., 54 (In the text)
• Brown, E. W. 1905, Memoirs R.A.S., 57 (In the text)
• Brown, E. W. 1908, Memoirs R.A.S., 59 (In the text)
• Brown, E. W. 1919, Tables of the motion of the Moon (New Haven: Yale University Press) (In the text)
• Capitane, N., Wallace, P. T., & Chapront, J. 2003, A&A, 412, 567 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Chapront, J. 1997, A&A, 109, 181 (In the text)
• Chapront, J. 2000, Celest. Mech. Dyn. Astr., 78, 75 [NASA ADS] [CrossRef] (In the text)
• Chapront, J., & Francou, G. 1996, http://cdsweb.u-strasbg.fr/cgi-bin/Cat?VI/87 (In the text)
• Chapront, J., & Chapront-Touzé, M. 1997, Celest. Mech. Dyn. Astr., 66, 31 [NASA ADS] [CrossRef] (In the text)
• Chapront-Touzé, M., & Chapront, J. 1983, A&A, 124, 50 [NASA ADS] (In the text)
• Chapront-Touzé, M., & Chapront, J. 1988, A&A, 190, 342 [NASA ADS] (In the text)
• Chapront, J., & Francou, G. 2003, A&A, 404, 735 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Chapront, J., Chapront-Touzé, M., & Francou, G. 2002, A&A, 387, 700 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Deprit, A., Henrard, J., & Rom, A. 1971a, AJ, 76, 269 [NASA ADS] [CrossRef] (In the text)
• Deprit, A., Henrard, J., & Rom, A. 1971b, AJ, 76, 723 [NASA ADS]
• Eckert, W. J., Jones, R., & Clark, H. K. 1954, Improved lunar ephemeris 1952-1959 (Washington D.C.: US Government Printing Office) (In the text)
• Fienga, A., Manche, H., Laskar, J., & Gastineau, M. 2006, A&A, in press (In the text)
• Henrard, J. 1979, Celest. Mech., 19, 337 [NASA ADS] (In the text)
• Henrard, J. 1980, Celest. Mech., 22, 335 [NASA ADS] (In the text)
• Henrard, J. 1981, Celest. Mech., 25, 417 [NASA ADS] (In the text)
• Hill, G. W. 1878a, Am. J. Math, 1, 5 [CrossRef] (In the text)
• Hill, G. W. 1878b, Am. J. Math, 1, 129 [CrossRef]
• Hill, G. W. 1878c, Am. J. Math, 1, 245 [CrossRef]
• Krasinsky, G. A. 2002, Commun. IAA RAS, 148, 1 (In the text)
• Kudryavtsev, S. M. 2004, J. Geodesy, 77, 829 [CrossRef] (In the text)
• Laskar, J. 1993, Phys. D, 67, 257 [CrossRef] (In the text)
• Laskar, J., Froeschlé, C., & Celletti, A. 1992, Phys. D, 56, 253 [CrossRef] (In the text)
• McCarthy, D. D., & Petit, G. 2003, IERS Conventions, IERS Technical Note 32 (Frankfurt am Main: Verlag des Bundesamts fur Kartografie und Geodasie) (In the text)
• McCuttcheon, R. A. 2003, Prepr. IAM RAS, 32, 7 (In the text)
• Pitjeva, E. V. 2001, Celest. Mech. Dyn. Astron., 80, 249 [NASA ADS] [CrossRef] (In the text)
• Pitjeva, E. V. 2003, Commun. IAA RAS, 156, 1 (In the text)
• Pitjeva, E. V. 2005, Sol. Sys. Res., 39, 176 [NASA ADS] (In the text)
• Schmidt, D. S. 1980, The Moon and the Planets, 23, 135 [NASA ADS] [CrossRef] (In the text)
• Sidlichovský, M., & Nesvorný, D. 1997, Celest. Mech. Dyn. Astr., 65, 137 [NASA ADS] [CrossRef] (In the text)
• Simon, J. L., Bretagnon, P., Chapront, J., et al. 1994, A&A, 282, 663 [NASA ADS] (In the text)
• Standish, E. M. 1998, JPL IOM 312.F-98-048 (In the text)
• Standish, E. M. 2003a, JPL IOM 312N-03-007 (In the text)
• Standish, E. M. 2003b, JPL IOM 312N-03-009
• Standish, E. M. 2006, JPL IOM 343R-06-002 (In the text)
• Standish, E. M., & Williams, J. G. 1981, ftp://nav.jpl.nasa.gov/pub/ephem/export/ (In the text)
• Standish, E. M., Newhall, X. X., Williams, J. G., & Folkner, W. F. 1995, JPL IOM 314.10-127 (In the text)
• Williams, J. G., Newhall, X. X., & Dickey, J. O. 1991, AJ, 241, L9 [NASA ADS] (In the text)