Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://astro.uni-altai.ru/~aw/DeltaT/1999A&A...343..624C.pdf Äàòà èçìåíåíèÿ: Tue Feb 12 11:25:14 2013 Äàòà èíäåêñèðîâàíèÿ: Fri Feb 28 04:44:07 2014 Êîäèðîâêà: Ïîèñêîâûå ñëîâà: apollo 15

Astron. Astrophys. 343, 624­633 (1999)

ASTRONOMY AND ASTROPHYSICS

Determination of the lunar orbital and rotational parameters and of the ecliptic reference system orientation from LLR measurements and IERS data
´ J. Chapront, M. Chapront-Touze, and G. Francou
Observatoire de Paris/DANOF ­ CNRS/URA 1125, 61 avenue de l'Observatoire, F-75014 Paris, France Received 27 April 1998 / Accepted 7 December 1998

Abstract. An analysis of Lunar Laser Ranging (LLR) observations from January 1972 till March 1998 is performed using the lunar theory ELP 2000-96 and the completed Moons' theory of the lunar libration. The LLR station coordinates, polar motion and Universal Time are provided by the International Earth Rotation Service (IERS). In Solution 1 the precession-nutation transformation is given by recent analytical theories, while in Solution 2 it is derived from the IERS daily corrections. Orbital and free libration parameters of the Moon, and coordinates of the reflectors are obtained in both cases. The position of the inertial mean ecliptic of J2000.0 with respect to the equator of the mean Celestial Ephemeris Pole (CEP) of J2000.0 (in Solution 1) and to the International Celestial Reference System (ICRS), the IERS celestial reference system, (in Solution 2) are fit. The position of the mean CEP equator of J2000.0 and of several dynamical reference planes and origins, with respect to ICRS, are derived from these fits (Fig. 1). The leading results are the following: 0. 057 60±0. 000 20 (in the equator) for the separation of the origin of right ascensions in ICRS from the ascending node of the inertial mean ecliptic of J2000.0 on the reference plane of ICRS, -0. 0460 ± 0. 0008 (in the ecliptic) for the separation of the latter point from the inertial dynamical mean equinox of J2000.0, -0. 015 19 ± 0. 000 35 (in the equator) for the separation of the inertial dynamical mean equinox of J2000.0 from the J2000.0 right ascension origin derived from IERS polar motion and Universal Time and from precise theories of precession-nutation, and 23 26 21. 405 22 ± 0. 000 07 for the inertial obliquity of J2000.0. A correction of -0. 3437 ± 0. 0040 /cy to the IAU 1976 value of the precession constant is also obtained (the errors quoted are formal errors). Key words: ephemerides ­ reference systems ­ Moon

and then from R at t2 to the station receiver O at t3 . DT C , the computed value of DT , is given by: DT C = t3 - t1 - DT1 (t3 ) + DT1 (t1 ) with t3 = t2 + 1 | BR(t2 ) - BO (t3 ) | +DT3 + DT4 c 1 t2 = t1 + | BR(t2 ) - BO(t1 ) | +DT3 + DT4 . c

(1)

c is the velocity of light and B stands for the barycenter of the solar system. DT1 is a relativistic time scale correction, DT3 the relativistic propagation correction (one way), and DT4 a tropospheric correction. t1 is derived from the observation time and t2 , t3 from Eqs. (1). In this analysis, the time scale TDB is adopted for time t (e.g. t1 , t2 , t3 ) and TAI (the time scale of DT ) for DT C . DT1 is evaluated from the first three terms of the series TDB-TT in (Fairhead & Bretagnon, 1990). DT3 is given by the following approximation of the classical formula (McCarthy, 1992) in the frame of the general relativity theory: D T3 = | OR | 4GmS 3 c | ST | + | SL | 2GmT 2 | TL | ln + c3 | TO | + | TR | - | OR |

where T is the terrestrial mass center, L the lunar one, and S the solar one, all the vectors being evaluated at time t1 . G is the constant of gravitation, mS and mT are the masses of the Sun and the Earth. DT4 is given by Marini and Murray's formula (McCarthy, 1992). The computation of DT C involves: ­ the components of TL in a celestial barycentric reference system, and the celestial barycentric coordinates of the Earth-Moon barycenter to a lesser accuracy, ­ the components of TO and TO in a terrestrial reference system, ­ the components of LR in a selenocentric reference system, ­ the rotations from the terrestrial and selenocentric reference system axes to the celestial ones,

1. Introduction Lunar laser ranging (LLR) stations provide normal points which may be used as observed values of the light time DT from the LLR station transmitter O at time t1 to a lunar reflector R at t2 ,
Send offprint requests to: J. Chapront


J. Chapront et al.: Determination of lunar parameters and ecliptic reference system orientation Table 1. Observations used Station McDonald 2.70m Old MLRS New MLRS Haleakala CERGA rubis CERGA beginning Jan. 1972 Aug. 1983 Feb. 1988 Apr. 1987 Apr. 1984 Oct. 1987 end Jan. 1985 Jan. 1988 March 1998 Aug. 1990 June 1986 March 1998 number 3249 483 1807 482 1190 5134

625

2. Orbital motion The lunar orbital motion is given by the improved analytical solution ELP 2000-82B plus numerical complements fit to the numerical integration DE245 of the Jet Propulsion Laboratory ´ (JPL), as described in (Chapront & Chapront-Touze, 1997). This solution, which is denoted now as ELP 2000-96, allows us to compute lunar geocentric rectangular coordinates xi referred to the inertial mean ecliptic of J2000.0, under the form: x = R3 (-W1 )X + . x is the 1-column matrix of coordinates xi . Ri is the rotation matrix with respect to the xi axis, with the same conventions as in (Folkner et al., 1994), x1 and x2 axes being arbitrary fixed orthogonal axes in the inertial mean ecliptic of J2000.0. W1 is the lunar mean mean longitude referred to the x1 axis. X and are 1-column matrices whose elements are respectively Poisson series and numerical complements; both do not depend on the precise direction of the x1 axis. Matrix X takes into account the rotation from the mean ecliptic of date (the original reference plane of ELP 2000-82B) to the mean ecliptic of J2000.0. xi are space coordinates of a BRS whose time coordinate is in TDB The expression of DT C involves also the barycentric position of the Earth-Moon barycenter but only through the differences of these vectors at t1 and t2 or t2 and t3 . The amplitudes of these differences do not depend on the precise directions of the xi axes. 3. Transformation from selenocentric coordinates to celestial coordinates The selenocentric axes adopted in this analysis are the principal axes of inertia. The transformation from the components i of LR in the selenocentric axes to its components yi in the ecliptic axes xi defined in Sect. 2 is: y = R3 (-W1 + 180 )M (p1 , p2 , , t) . y and are the 1-column matrices whose elements are respectively yi and i . W1 is the same as in Sect. 2, M (p1 , p2 , , t) is a matrix function of the libration variables p1 , p2 , and of time t. The adopted solution for the libration is Moons' theory (1982, 1984) with analytical and numerical complements as described in (Chapront et al., 1998). The main problem of Moons' theory accounts for a rigid Moon and involves the lunar physical pa2 rameters = (C - A)/B , = (B - A)/C , C /mL RL , and unnormalized Cn,k , Sn,k (0 k n, 3 n 4), under a literal form. The values assigned to these parameters in the present analysis are those of the JPL numerical integration DE245, ex2 cept for C /mL RL . They are quoted in Table 2 except for the values of C4,k , S4,k which are those of Ferrari et al. (1980). 2 We have assumed that the values of , , C /mL RL of Table 2 involve the constant perturbation due to the lunar rotation. The tidal perturbations have been computed, assuming a constant time delay L , with the values of lunar tidal parameters L and k2,L (Love number) of DE245, quoted in Table 2. Similar to the lunar orbital motion, matrix M takes into account the rotation

­ relativistic corrections for the transformations of coordinates in terrestrial and selenocentric reference systems to coordinates in a barycentric reference system (BRS). In this analysis, the coordinates of the LLR stations in the International Terrestrial Reference System (ITRS), the terrestrial reference system of the International Earth Rotation Service (IERS), are derived from ITRF94 (Boucher et al., 1996) and corrected for the Earth's deformations due to Earth tides, ocean tides, and pressure anomaly, following the recommendations of the IERS Standards 1992 (McCarthy, 1992). Files EOP(IERS)97 C04, provided by IERS, yield the orientation of the ITRS axes both with respect to the instantaneous axes defined by the celestial ephemeris pole (CEP) and "true equinox of date" and with respect to the axes of the International Celestial Reference System (ICRS), the IERS celestial reference system. In the first case, precise analytical theories of precession and nutation subsequently allow one to refer to the celestial system of axes defined by the mean CEP equator of J2000.0 and the fixed origin of right ascensions in this plane derived from the IERS data for polar motion and UT1 and from precession-nutation theory (in this paper, MCEP system). We use solutions for the lunar libration and for the lunar orbital motion which preserve the analytic advantages, both referred to the mean ecliptic of J2000.0 in the inertial sense as defined by Standish (1981). This allows one to fit simultaneously parameters of the free libration, orbital parameters of the Moon and of the Earth-Moon barycenter (including the lunar tidal secular acceleration), selenocentric coordinates of the lunar reflectors, and positional parameters of the inertial mean ecliptic of J2000.0 with respect to the ICRS axes and to the MCEP system. In the latter case, a correction to the IAU 1976 value of the precession constant is also fit. This analysis uses the normal points obtained by the LLR stations of McDonald (three locations), Haleakala, and CERGA (two successive instruments at the same location) from January 1972 to March 1998. The number of normal points and the time spans covered are given in Table 1. For all the stations, except Haleakala, the normal points involve the same location for transmitter O and receiver O . We mention that, since the coordinates of the new MLRS and Haleakala receivers are not given in ITRF94, we have computed them by adding respectively to the coordinates of the old MLRS receiver and of the Haleakala transmitter in ITRF94 the differences between the values of the corresponding coordinates in (Newhall et al., 1991).


626

J. Chapront et al.: Determination of lunar parameters and ecliptic reference system orientation

Table 2. Values adopted for the lunar physical parameters (from DE245 2 except for C /mL RL ) = 0.631 619 133 10-3 = 0.227 885 980 10-3 2 C /mL RL = 0.394 872 400 C3,0 = -0.086 802 10-4 C3,1 = 0.307 083 10-4 S3 C3,2 = 0.048 737 10-4 S3 C3,3 = 0.017 161 10-4 S3 k2,L = 0.029 920 L = 0.164 846days

,1 ,2 ,3

= 0.046 115 10-4 = 0.016 975 10-4 = -0.002 844 10-

4

from the mean ecliptic of date (the original reference plane of Moons' theory) to the mean ecliptic of J2000.0. In the present analysis, the relativistic corrections for the transformation of selenocentric coordinates to celestial coordinates have not been taken into account, so i are, in fact, BRS coordinates. 4. Transformation from terrestrial coordinates to celestial coordinates The transformation from the components i of TO or TO in the ITRS axes to its components zi in the ecliptic celestial axes xi defined in Sect. 2 is: z = R1 ()R3 ()P -1 N -1 âR3 (-GS T )R1 (yp )R2 (xp ) + Dx. z and are the 1-column matrices whose elements are respectively zi and i . Dx is the 1-column matrix of the relativistic corrections for the transformation of space coordinates in a terrestrial reference system whose time coordinate is in TCG to coordinates in a BRS whose time coordinate is in TDB. The expression used for Dx is similar to Eq. (20) in (Martin et al., 1985) reduced to the solar contribution and with L = 1.5505 10-8 instead of 1.481 10-8 (since ITRF94 uses TCG instead of TT). xp and yp are the polar motion components giving the tangential coordinates of the CEP with respect to ITRS. GS T is the Greenwich true sidereal time given by: GS T = GM S T (UT1) + cos A + 0. 00264 sin +0. 000063 sin 2 + (p cos A ) t

(2)

where GM S T (UT1) is the mean sidereal time given by Aoki et al. (1982) as a function of Universal Time UT1, A the mean obliquity of date, the nutation in longitude, the tropic mean longitude of the lunar node. (p cos A ) t is a correction induced by the correction p to the IAU 1976 value of the precession constant, when the mean motions of the theory are referred to a J2000.0 fixed origin of longitudes, in order not to introduce a drift in UT (Williams & Melbourne, 1982; Zhu & Mueller, 1983; Williams, 1994); t is the time reckoned from J2000.0. Matrix R3 (-GS T )R1 (yp )R2 (xp ) rotates the ITRS axes to the celestial instantaneous axes, two of them pointing respectively towards the CEP and the "true equinox of date". In this

analysis xp , yp , and the differences UT1 - UTC have been computed by interpolating the daily values of the IERS files EOP(IERS)97 C04. Matrix P -1 N -1 rotates the celestial instantaneous axes to a J2000.0 fixed celestial "equatorial" system of axes, i.e. a system of axes in which the reference plane (the plane of the first two axes) is close to the mean CEP equator of J2000.0. We have considered two cases. In the first case, N is the nutation matrix provided by an analytical theory and P is the equatorial precession matrix, between J2000.0 and the date of observation, yielded by an analytical theory in which the precession constant is a fitted parameter. In Eq. (2), p is the difference between this fitted parameter and the IAU 1976 value, and is provided by the same nutation theory as N . So the fixed celestial "equatorial" system of axes is the MCEP system. For nutation, the present analysis uses the ZMOA 1990 solution (Herring, 1991). For precession, we have adopted Williams' expressions (1994) completed by the derivatives with respect to the precession constant and obliquity from (Simon et al., 1994). The set of fitted parameters obtained in this way is denoted as Solution 1 (Sol. 1). In the second case, P is the equatorial precession matrix between J2000.0 and the date of observation from (Lieske et al., 1977), involving the IAU 1976 value of the precession constant, and N is the nutation matrix in which the nutation in longitude and the nutation in obliquity are computed by adding to the IAU 1980 expressions (Seidelmann, 1982) the corrections d and d provided in files EOP(IERS)97 C04. The same is also introduced in Eq. (2). d and d involve nutation corrections, precession corrections (e.g. a correction to the IAU 1976 value of the precession constant), and a small rotation, so that P -1 N -1 rotates the celestial instantaneous axes to the ICRS axes, which are slightly different from those of the MCEP system. In this second case, the precession constant is not fit and p is set to zero in Eq. (2), the correction (p cos A ) t being included in d . The so obtained set of fitted parameters is denoted as Solution 2 (Sol. 2). In Sol. 1, is the inclination (MCEP) of the J2000.0 inertial mean ecliptic to the mean CEP equator of J2000.0 and I is the arc o(MCEP)2000 (MCEP). o(MCEP) is the origin of right ascensions in the MCEP system defined in Sect. 1; I 2000 (MCEP) is the ascending node of the J2000.0 inertial mean ecliptic on the mean CEP equator of J2000.0 (i.e. the inertial dynamical mean equinox of J2000.0). In Sol. 2, is the inclination (ICRS) of the J2000.0 inertial mean ecliptic to the I ICRS reference plane and is the arc o(ICRS)2000 (ICRS) between o(ICRS), the origin of right ascensions in ICRS, and I 2000 (ICRS), the ascending node of the J2000.0 inertial mean ecliptic on the ICRS reference plane. So the arbitrary ecliptic x1 I axis defined in Sect. 2 points towards 2000 (MCEP) in Sol. 1 I and towards 2000 (ICRS) in Sol. 2 (see Fig. 1 in Sect. 6.4). 5. The fits In addition to parameters and described in Sect. 4, each fit provides:


J. Chapront et al.: Determination of lunar parameters and ecliptic reference system orientation Table 3. Number of observations n retained in the groups and RMS of their post fit residuals in cm for one way range. Groups McDonald 1972­1986 McDonald 1987­1998 CERGA, rubis 1984­1986 CERGA 1987­1998 Haleakala 1987­1990 n 3513 1823 1166 4979 462 RMS 34.7 5.0 18.2 4.8 11.1
(0) (0) (0)

627

Table 4. Corrections to values of orbital parameters fit to DE200. Units: /cy for and n , and arcsecond for the other variables Variable
(0) W1 (0) W2 (0) W3

Sol. 1 -0. -0. -0. -0. 0. 0. -0. -0. 0. 0. 120 066 106 391 000 000 075 060 033 000 35 13 76 65 61 18 89 20 19 06 ± ± ± ± ± ± ± ± ± ± 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 000 000 000 001 000 000 000 000 001 000 34 33 96 15 03 00 33 51 16 01 0. 0. 0. 0. 0. 0. -0. -0. 0. 0. - - - - 074 020 063 390 000 000 029 013 032 000 31 04 10 40 69 18 81 12 14 06

Sol. 2 ± ± ± ± ± ± ± ± ± ± 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 000 000 000 001 000 000 000 000 000 000 21 20 93 00 03 00 20 41 94 01

­ the geocentric lunar orbital parameters W1 , W2 , W3 (values of the mean mean longitude, mean longitude of perigee, and mean longitude of node in J2000.0), , , E (sidereal mean motion in J2000.0, constant relative to the sine of half inclination, eccentricity constant); ­ the heliocentric orbital parameters of the Earth-Moon barycenter T (0) , (0) (values of the mean mean longitude and mean longitude of the perihelion in J2000.0), n , e (sidereal mean motion in J2000.0 and eccentricity constant); (2) (1) (1) ­ the bias parameters W1 , W2 , W3 (observed corrections to the computed coefficient of the quadratic term of the lunar mean longitude and to the computed mean motions of perigee and node); ­ the free libration parameters 2P , 2Q, 2R, p(0) , q (0) , r(0) (parameters tied to the coefficients of the main free libration terms and values of the free libration arguments in J2000.0 in Moons' theory). , the tidal part W1 yields an observed value of W1 of the coefficient of the quadratic term of the mean longitude (half tidal secular acceleration), the other contributions to this quadratic term being given with enough precision by the theory. All the fitted angles, except the free libration ones, are reI I ferred to 2000 (MCEP) in Sol. 1, and to 2000 (ICRS) in Sol. 2. (0) (2) The introduction of the fitted values of W1 , , and W1 in the theory allows us to compute the values of W1 at the weighted mean epoch of observations, referred respectively to I I 2000 (MCEP) and 2000 (ICRS). The difference of these valI I ues yields 2000 (MCEP)2000 (ICRS), measured in the ecliptic, which, combined with the determinations of and , allows us to tie the MCEP axes to the ICRS axes. The fits are iterative weighted fits performed by the leastsquares method. The same weight is assumed for all the observations of each group given in Table 3. It is computed according to the RMS of the post-fit residuals of the group at the previous iteration. Every observation whose absolute value of the residual is greater than three times the RMS of the residuals of its group is disregarded in the next iteration. The groups depend on the time span and on the LLR station. Table 3 gives the number of observations retained in each group and the RMS of their post fit residuals in centimeter for one way range (i.e. | DT - DT C | /2c) in Sol. 1. The values for Sol. 2 are very similar. The mean epoch of observations is 1 February 1988, the weighted mean epoch is 1 April 1992.
(2) (2,T)

E T0
(0)

n e

Table 5. Corrections to the values of the orbital parameters of the theory. Units: /cy for and n , and arcsecond for the other variables Variable Sol. 0.164 -0.080 E 0.018 n -0.031 e -0.128 1 Sol. 39 0.165 05 -0.079 07 0.018 01 -0.032 73 -0.128 2 64 97 07 06 73

Table 6. Fitted value of the tidal part of the quadratic term of the mean longitude (in /cy2 ) and observed corrections to the mean motions of perigee and node (in /cy) Variable
(2,T) W1 (1) W2 (1) W3

Sol. 1

Sol. 2

-12.900 7 ± 0.002 3 -12.890 8 ± 0.002 2 0.026 95 ± 0.001 37 0.029 24 ± 0.001 32 -0.194 53 ± 0.013 65 -0.266 42 ± 0.013 47

6. The results 6.1. Orbital motion Table 4 gives, for Sol. 1 and Sol. 2, the fitted values of the parameters of the orbital motions of the Moon and of the EarthMoon barycenter in the form of corrections to previous values fit to the JPL numerical integration DE200 and adopted in the lunar ´ ephemeris ELP 2000 (Chapront-Touze & Chapront, 1983). The I I origin of angles is 2000 (MCEP) in Sol. 1, and 2000 (ICRS) in Sol. 2. Table 5 gives, for , , E , n and e , the differences between the new fitted values and the values of the theory, i.e. the corrections to be added to the constants introduced in the main problem series and mean motions of ELP 2000-82B by means of the derivatives of the series coefficients and mean motions. Table 6 gives the observed corrections to the mean motions of perigee and node computed with the values of , , E , n , and e derived from Tables 4 and 5, and the values of the tidal quadratic term of the mean longitude obtained in the two fits. The value of this quadratic term computed with our analytic model of tidal perturbations and the values of the physical parameters involved in the JPL numerical integration 2 DE403 is -12. 7898 /cy . From the results of Tables 4, 5, and


628

J. Chapront et al.: Determination of lunar parameters and ecliptic reference system orientation Table 7. Differences between the values of orbital elements fit in this paper and values fit to DE403, after correction of the origin. Units: /cy2 (2) for W1 , /cy for and n , and arcsecond for the other variables Variable
(0) W1 (0) W2 (0) W3

´ 6, and from (Chapront-Touze & Chapront, 1988), we derive the following expressions for the angular mean elements of the Moon and of the Earth-Moon barycenter, referred to the inertial mean ecliptic of J2000.0. t is the time (TDB scale) reckoned from J2000.0 (Julian date 2 451 545.0) in Julian centuries. W1 = W2 = W3 = T= =
I Sol. 1 (origin = 2000 (MCEP)) 218 18 59. 835 36 + 1 732 559 343. 344 39 t -6. 841 7 t2 + 0. 006 604 t3 - 0. 000 031 69 t4 , 83 21 11. 608 62 + 14 643 420. 337 8 t - 38. 263 9 t2 -0. 045 047 t3 + 0. 000 213 01 t4 , 125 02 40. 291 40 - 6 967 919. 720 7 t + 6. 359 3 t2 +0. 007 625 t3 - 0. 000 035 86 t4 , 100 27 59. 144 70 + 129 597 742. 309 0 t - 0. 020 2 t2 +0. 000 009 t3 + 0. 000 000 15 t4 , 102 56 14. 367 33 + 1 161. 228 3 t + 0. 532 7 t2 -0. 000 138 t3 . I Sol. 2 (origin = 2000 (ICRS)) 218 18 59. 881 40 + 1 732 559 343. 345 64 t -6. 831 8 t2 + 0. 006 604 t3 - 0. 000 031 69 t4 , 83 21 11. 654 71 + 14 643 420. 3367 t - 38. 263 9 t2 -0. 045 047 t3 + 0. 000 213 01 t4 , 125 02 40. 335 06 - 6 967 919. 7920 t + 6. 359 3 t2 +0. 007 625 t3 - 0. 000 035 86 t4 , 100 27 59. 190 78 + 129 597 742. 3079 t - 0. 020 2 t2 +0. 000 009 t3 + 0. 000 000 15 t4 , 102 56 14. 414 41 + 1 161. 2283 t + 0. 532 7 t2 -0. 000 138 t3 .

Sol. 1 -0. 0. -0. -0. -0. 0. 0. -0. -0. 0. -0. 0. 0. 000 000 006 011 000 000 003 175 064 000 004 006 000 38 31 87 85 70 04 9 2 5 27 59 30 06

Sol. 2 -0. 0. -0. -0. -0. 0. 0. -0. -0. 0. -0. 0. 0. 000 000 009 010 000 000 002 246 054 000 003 005 000 34 40 21 60 62 04 8 5 6 35 51 25 06

E (1) W2 (1) W3 (2) W1 T (0)
(0)

n e

W1 = W2 = W3 = T= =

Table 7 gives the differences between the values of this paper and values fit to the JPL numerical integration DE403 ´ from (Chapront & Chapront-Touze, 1997). For angular vari(0) (0) (0) (0) (0) , the differences have been ables, W1 , W2 , W3 , T , I corrected for the separation between the origins 2000 (MCEP) I I or 2000 (ICRS) and 2000 (DE403), derived from Table 11. For (1) (1) (2) W2 , W3 , and W1 the differences are those of the total values (computed values + observed corrections or bias). The values of , and n obtained in this analysis are closer to the values fit to DE403 than to the older values fit to DE200 because both DE403 and the present analysis are based on sets of LLR observations covering much larger time spans than DE200. The errors quoted in Tables 4 and 6 are 1- errors provided by the least-squares fits. The errors of the angular variables, (0) except W3 , are significantly smaller in the second solution because the origin is more accurately determined if the precession constant is not fit; the errors of a test fit similar to Sol. 2 but with a residual precession constant included in the solution (see Sect. 6.4) are similar to those of Sol. 1. Also, the errors were larger in earlier fits performed with the method described in this paper but with a less precise libration theory (i.e. without the numerical complements).

All the variables are not determined with the same accuracy. In the series of the lunar center motion (longitude, latitude and distance) W2 , W3 and T appear through their differences with respect to W1 in Delaunay arguments l, F , and D, and appears through its difference with respect to T in l . The range station-reflector depends mainly on the distance, but also on the sines and cosines of longitude and latitude multiplied by the equatorial parallax; then, considering the mean Earth-Moon distance as a scale factor, the larger trigonometric terms involving l, W1 , D, F , and l have amplitudes of 0.055, 0.017, 0.0096, 0.0015, and 0.0005 respectively; E , , and e are multiplied by trigonometric terms in l, F and 2D - l whose amplitudes are 1, 0.034, and 0.032 respectively. The errors follow the relative magnitudes of the above amplitudes, but one must also take into account the length of the time span covered by the observations and the weights. (1) We note that the observed correction W2 to the computed mean motion of perigee is smaller by about 0 .025 /cy ´ than the bias obtained in Chapront & Chapront-Touze (1997) by comparison to the JPL numerical integration DE403, though in Table 7 the total values of the mean motion differ by less than 0 .0040 /cy. The reason seems to be that the total mean motion (1) (1) W2 is well determined and that W2 corrects errors in the (1) computed value of W2 resulting from errors in the metric variables, especially in : the value 0 .025 /cy corresponds to an error of 0 .000 66 in . So the realistic error in could be much larger than the formal one, even though the factor of 20 resulting from the rough calculation above seems too large. If we suppose a factor of 10 between the formal error and the realistic one in , (0) (1) this factor also applies to the errors in W3 and W3 , and (1) W3 being tied, and the large observed correction W3 becomes unsignificant. The fact that the time span of observations covers only 1.5 period of the node and the large weights given to the most recent observations certainly contributes to the large (1) realistic error in W3 and consequently to the realistic errors (0) in W3 and .


J. Chapront et al.: Determination of lunar parameters and ecliptic reference system orientation Table 8. Fitted values of the libration parameters Variable p q (0) r(0) 2P 2Q 2R
(0)

629

Sol. 1 224. 342 1 ± 0. 01 161 290 8 ± 0 01 . . 94 039 2 ± 1 10 . . 0.292 61 ± 0.000 5.218 74 ± 0.002 0.023 84 ± 0.000


Sol. 2 2 4 4 0 1 4 1 7 9 6 1 0 224. 161 . 97 . 0.29 5.22 0.02


3 2 0 2 3 5

41 57 00 60 39 52

5 ± 0 01 . 7 ± 0 01 . 8 ± 1 09 . ± 0.000 ± 0.002 ± 0.000

2 4 2 0 1 3

1 8 1 6 2 9

of this paper in Moons' literal series, we obtain the following expressions for the parts of the libration series which depend on the free libration (terms greater than 0. 005). The values of the lunar physical parameters involved are given in Table 2. Sol. 1 = -3. 312 sin(q - 0 01) + 0. 034 sin(q - l) . +0. 026 sin(p - F ) + 0. 023 sin(p + F ) +0. 024 sin(F + r) - 0. 022 sin(q + l), = 8. 197 cos(q + 0 01) - 0. 035 cos(q - l) . -0. 026 cos(p - F ) + 0. 023 cos(p + F ) +0. 024 cos(F + r) - 0. 022 cos(q + l). Sol. 2 = -3. 315 sin(q - 0 01) + . +0. 026 sin(p - F ) + 0. +0. 023 sin(p + F ) - 0. = 8. 204 cos(q + 0 01) - 0 . -0. 026 cos(p - F ) + 0. +0. 025 cos(F + r) - 0.

p

F 1

p Except for the difference of origin, the differences between Sol. 1 and Sol. 2 are mostly due to the inaccuracies of the nutation model ZMOA 90 involved in Sol. 1. The comparison of Columns 2 and 3 of Table 7 with the formal errors of Tables 4 and 6 shows that the differences between Sol. 1 and Sol. 2 do (0) (1) not exceed the formal errors except for , W3 , and W3 , for which the realistic errors are much larger than the formal (2) ones, and except for W1 and . For these last two variables, the realistic errors are also probably larger than the formal ones, (2) and from this comparison, we propose a factor of 5 for W1 and 2 for . For the others variables, the formal errors seem to be realistic. 6.2. Free libration parameters Table 8 gives the fitted values of the libration parameters. The quoted errors are 1- errors provided by the least-squares method; they are almost the same in the two fits. For the angu obtained in the two fits lar parameters and for 2P , the values Q differ by less than the 1- errors. For 2 , the coefficient of a term whose period is 75 years, and for 2R, the coefficient of a term whose period is close to the draconitic one, the values differ respectively by twice three times the errors. and The values of 2P , p0 , 2R, r0 are dependent on the solution used for the libration. The introduction of numerical com plements has modified the fitted values of p0 and 2P , and of r0 by about 150 times and 20 times the 1- errors respectively, because of the introduction of terms with frequencies close to p and F + r frequencies in the numerical complements. In the opposite, the values of q0 and 2Q have been changed by only four times the errors quoted in Table 8. In fact, though the period of q is long, 2Q is well deter mined. The value of 2Q derived from a comparison of the libration solution used in this paper with the libration part of the JPL numerical integration DE403 over a time span of 3 centuries (Chapront et al., 1998) is 5.2095. It differs by 7 times the 1- error from the value of Sol. 2 which is probably better because the time span of the LLR observations involved is larger in the present analysis than in DE403. Furthermore, we note that the value derived from the coefficient of the cos q term F of p2 in (Calame, 1977) is 5.0, based on the first six years of LLR observations only. This good agreement may be due to the fact that the frequency of q is known by the theory and that the other terms with long periods in the libration variables p1 and p2 are very small. By introducing the values of 2P , 2Q, 2R

F 2

p

F 1

p

F 2

0. 034 sin(q - l) 026 sin(F + r) 022 sin(q + l), . 035 cos(q - l) 023 cos(p + F ) 022 cos(q + l).

F

Sol. 1 and Sol. 2 = 1. 819 sin p + 0. 088 sin(q + 2F - 2l + 36 ) +0. 077 sin(q + F ) + 0. 069 sin(q + F - l) -0. 033 sin(q - F ) + 0. 015 sin(q - F + l).

6.3. Reflector coordinates Table 9 gives the fitted values of the reflector coordinates i referred to the lunar principal axes of inertia (PA coordinates). These values must be considered as space coordinates in a BRS whose time coordinate is in TDB, i.e. no scale factor and Lorentz contraction to convert to a selenocentric reference system has been applied. The errors quoted are 1- errors of the least-squares method. The realistic errors estimated from the differences between Sol. 1 and Sol. 2 are probably about 15 centimeters. Reflector coordinates are commonly given in the mean Earth/rotation axes system (MA coordinates). Following (Ferrari et al., 1980), the MA axes are shifted from the PA axes by means of the constant terms of the libration variables p1 , p2 , . This means that p1 and p2 are the components of the unit vector of the polar MA axis on the equatorial PA axes, and is the angle measured from the first equatorial MA axis to the ascending node of the PA equator on the MA equator and then from the node to the first PA axis. The libration theory and the values of the lunar physical constants adopted in the present analysis give p1 = -78. 931 611 286, p2 = 0. 290 177 293, = 66. 189 835 222, from which we derive: = R1 (-0. 264 848 351)R2 (-78. 931 702 266) R3 (-66. 189 784 547)


630

J. Chapront et al.: Determination of lunar parameters and ecliptic reference system orientation

Table 9. Fitted values of the reflector coordinates i referred to lunar principal axes of inertia. Unit: meter Reflector Sol. 1 Apollo 11 1 1 591 963.113 ± 0.001 2 690 708.509 ± 0.037 3 21 006.379 ± 0.025 Sol. 2 1 591 963.112 ± 0.001 690 708.386 ± 0.037 21 006.461 ± 0.025

squares method. The results of Table 11 and other results of this section are represented in Fig. 1. I 2000 (MCEP) is the inertial dynamical mean equinox of J2000.0; the rotational dynamical mean equinox of J2000.0 R 2000 (MCEP) is the ascending node of the J2000.0 mean ecliptic, in the rotational sense as defined by Standish (1981), on the mean CEP equator of J2000.0. Following Standish, the separation between the two equinoxes is given by:
I R 2000 (MCEP)2000 (MCEP) = 0. 09366

Lunakhod 2 1 1 339 357.864 ± 0.006 2 801 879.560 ± 0.044 3 756 361.110 ± 0.034

Apollo 15 1 1 554 676.642 ± 0.005 2 98 104.053 ± 0.034 3 765 008.184 ± 0.016

Apollo 14 1 1 652 692.801 ± 0.002 1 652 692.801 ± 0.002 -520 989.116 ± 0.037 -520 989.243 ± 0.037 2 3 -109 727.749 ± 0.024 -109 727.677 ± 0.024 1 554 676.629 ± 0.005 98 103.926 ± 0.034 765 008.265 ± 0.016 1 339 357.854 ± 0.006 801 879.410 ± 0.044 756 361.187 ± 0.034

measured in the mean CEP equator of J2000.0. So we deduce from Table 11: o(MCEP)
R 2000

(MCEP) = 0. 078 47 ± 0. 000 35

Denoting by 2I 00 (MCEP) and 2R00 (MCEP), the projections 0 0 R I of 2000 (MCEP) and 2000 (MCEP) on the reference plane of ICRS, the separation on the equator between the origin of right ascensions in ICRS and the rotational dynamical mean equinox of J2000.0 is given by: o(ICRS)2R00 (MCEP) 0 I = o(ICRS)2000 (ICRS) I I +2000 (ICRS)2000 (MCEP) cos I R +2000 (MCEP)2000 (MCEP) This quantity is also the component of the unit vector pointing towards the rotational dynamical mean equinox of J2000.0 on the x2 axis of ICRS. The values of Table 11 give: o(ICRS)2R00 (MCEP) = 0. 0783 ± 0. 0009, 0 o(ICRS)2I 00 (MCEP) = -0. 0154 ± 0. 0009. 0 Our value of the separation on the equator between the origin of right ascensions in ICRS and the rotational dynamical mean equinox of J2000.0 is very close to the result 0. 078 ± 0. 010 of Folkner et al.(1994). Note that the origin of right ascensions in ICRS is closer to the projection of the inertial mean equinox of J2000.0 than to the projection of the rotational one by a factor of five. The separation on the equator between o(MCEP) (the origin of right ascensions in the mean CEP equator of J2000.0 derived, by means of precession-nutation transformation, from the "true equinox of date" as it results from the polar motion and UT1 provided by IERS) and o(ICRS), derived from Table 11, is almost zero. The test fit 1 is a similar fit as that of Sol. 1 but using the precession expressions of Simon et al. (1994) instead of those of Williams (1994). The leading difference between the two sets of expressions concerns the obliquity, with: A (Williams, t) - A (Williams, J 2000.0) = A (Simon, t) - A (Simon, J 2000.0) - 0. 024 40t

Table 10. Reflector coordinates i referred to mean Earth/rotation axes for Sol. 1. Unit: meter Reflector Apollo 11 1 591 Apollo 14 1 652 Apollo 15 1 554 Lunakhod 2 1 339 1 749. 817. 937. 389. 307 691 789 - 520 714 98 814 802 2 219. 458. 601. 308. 304 20 603 - 110 958 764 344 755 3 398. 360. 413. 849. 151 912 336 649

and are respectively the 1-column matrix of the PA and MA coordinates. Table 10 gives the MA coordinates of the reflectors for Sol. 1. The rough comparison of the values of Table 10 with the values given by Williams et al. (1996) yields a distance of 2.8 meters between the two positions on the lunar surface of Apollo 15, which is the most observed reflector, and distances of 2.6, 2.0, and 3.0 respectively between the positions of the other reflectors. Nevertheless we mention that probably the MA axes are not rigorously the same in the two papers. The fitted positions of the reflectors depend on the libration solution, and the introduction of the numerical complements in the solution of the lunar libration has changed the locations of the first two reflectors by 0.3 and 0.6 meter respectively, and the locations of the last two ones by 1.8 meter. 6.4. Orientation of celestial axes Columns 2 and 3 of Table 11 give the fitted values of I (R) and of o(R)2000 (R) measured in the reference plane of R, R standing for the MCEP system (Sol. 1) or ICRS (Sol. 2) as described in Sect. 4. Column 4 gives the separaI I tion 2000 (ICRS)2000 (MCEP) measured in the inertial mean ecliptic of J2000.0 and derived by the method of Sect. 5. The quoted errors are (or are derived from) 1- errors of the least-

(3)

(t reckoned from J2000.0 in Julian centuries), and is due to the motion of the equator. If we assume that LLR observations provide an accurate position of the mean CEP equator at the weighted mean epoch of observations, the obliquity at


J. Chapront et al.: Determination of lunar parameters and ecliptic reference system orientation

631

Fig. 1. Positions of several reference planes and origins of right ascensions or longitudes in the tangential plane to the celestial sphere Table 11. Position angles of the inertial mean ecliptic of J2000.0 with respect to "equatorial" celestial systems R (arcseconds) and fitted correction to the IAU 1976 value of the precession constant ( /cy) R (R) - 23 26 21 405 407 410 410 409 408 22 06 81 66 28 83 ± ± ± ± ± ± 0. 0. 0. 0. 0. 0. 000 000 000 000 000 000 o(R) - - - - - - 0. 0. 0. 0. 0. 0. 015 014 057 056 052 092
I 2000

(R) 000 000 000 000 000 000 35 35 20 35 01 16



I 2000

(ICRS)

I 2000

(R)

p -0.3437 ± 0.0040 -0.3382 ± 0.0040

MCEP 0. /Test fit 1 0. ICRS 0. /Test fit 2 0. DE403 0. DE200 0.

07 08 07 08 00 06

19 67 60 19 94 45

± ± ± ± ± ±

0. 0. 0. 0. 0. 0.

0.0460 ± 0.0466 ± 0 0.0015 ± 0.0069 ± -0.0339 ±

0.0008 0.0008 0.0007 0.0004 0.0011

this epoch (1 April 1992) is well determined and independent of the precession expressions. So Eq. (3) gives at J2000.0: A (Williams, J 2000.0) - A (Simon, J 2000.0) = -1.90 mas This quantity differs from the "observed value" derived from Table 11: (MCEP) - (MCEP/Test fit 1) = -1.84 ± 0.15 mas by less than 1- . This result yields a verification of the internal consistency of our method. Note that the test fit 1 yields for the separation on the equator between the origin of right ascensions in ICRS and the inertial dynamical mean equinox of J2000.0 o(ICRS)
I 2000

we have:
II D 2000 (ICRS/Test fit 2) = t(p + dp) II D 2000 (ICRS) = tp

and by subtraction:
I I 2000 (ICRS)2000 (ICRS/Test fit 2) = t dp

(4)

(MCEP/Test fit 1) = -0. 0148 ± 0. 0009

which differs by less than 1- from the value derived from Sol. 1. Column 5 of Table 11 gives the corrections to the IAU 1976 value of the precession constant obtained in Sol. 1 and by test fit 1; the two values differ by less than 2- . No precession constant value has been fit in Sol. 2, but in the test fit 2 a correction dp to the precession constant value p included in the IERS d and d has been determined in addition to the terms fit in Solution 2. The value obtained is -0. 0192 ± I I 0. 0040 /cy. The value of 2000 (ICRS)2000 (ICRS/Test fit 2) given in Table 11 is correlated with dp for the following reason. LLR observations provide an accurate position of the "mean plane" derived from the ICRS reference plane by means of the precession transformation between J2000.0 and the mean epoch of observations. The ascending node of the J2000.0 inertial I mean ecliptic on this mean plane D is accurately determined and is independent of the value of the precession constant. Then

On 1 April 1992, the right hand member of Eq. (4) amounts to 0. 0015. This value is equal to the value given in Table 11. This result yields another verification of the internal consistency of our method. It is also consistent with the precision of d (0.6 mas for d sin = 1.5 mas for d ) stated by IERS for 1992. ´ ´ (Chapront & Chapront-Touze, 1997) and (Chapront-Touze I & Chapront, 1983) give fitted values of (R),of o(R)2000 (R), I and of the lunar mean longitude referred to 2000 (R), R standing for the reference frames defined by the JPL numerical integrations DE403 and DE200. (R) is the inclination of the inertial mean ecliptic of J2000.0 on the reference plane of R, I o(R) is the origin of right ascensions in R and 2000 (R) is the node of the inertial mean ecliptic of J2000.0 on the refI erence plane of R. The values of (R) and o(R)2000 (R) are quoted in Table 11. From the lunar mean longitudes referred I I to the origins 2000 (DE403) and 2000 (DE200), evaluated at the mean epoch of the observations involved in the numerical integrations (set here to 1 January 1985 and 1 January 1975 respectively), and from Sol. 2 of this paper, we obtain, by the I I method of Sect. 5, the values of 2000 (ICRS)2000 (DE403) I I and 2000 (ICRS)2000 (DE200), measured in the inertial mean ecliptic of J2000.0, as quoted in Table 11. If we denote by o (DE200) the projection of o(DE200) on the mean CEP equator of J2000.0, we derive from Table 11 the


632

J. Chapront et al.: Determination of lunar parameters and ecliptic reference system orientation Table 12. Rotations between ICRS axes and other "equatorial" celestial systems of axes R, and comparison with the results of other authors: (F) Folkner et al. (1994), (I) IERS conventions 1996 R MCEP (I) MCEP/Test fit 1 (F) DE403 DE200 (F) 1 (mas) 5. 5. 3. 1. 1. 2. 2 6 1 8 9 5 0 ± ± ± ± ± ± ± 0. 0. 0. 0. 0. 0. 2 2 2 2 2 1 2 - - - - - 2 (mas) 18. 17. 18. 19. 2. 13. 12 3 3 5 2 7 5 ± ± ± ± ± ± ± 0. 0. 0. 0. 0. 0. 3 3 (mas) 4 -0.2 2 4 -0.2 3 2 1. 7 5 3. 7 6 ± 1.2 ± 1.2 ± 0.6 ± 1.4 ±3

following values of the separations on the equator between the origin of right ascensions of DE200 and the following origins: o(MCEP), the inertial dynamical mean equinox of J2000.0, and the rotational dynamical mean equinox of J2000.0: o (DE200)o(MCEP) = -0. 0039 ± 0. 0023 R o (DE200)2000 (MCEP) = 0. 0746 ± 0. 0019 I o (DE200)2000 (MCEP) = -0. 0191 ± 0. 0019 So, though o(DE200) has been constructed as the ascending node of the rotational mean ecliptic of J2000.0 on the reference plane of DE200 in order to represent the rotational dynamical mean equinox of J2000.0 (Standish, 1982), its projection on the mean CEP equator of J2000.0 is closer to the present position of the inertial dynamical mean equinox of J2000.0 than to that of the rotational one. This is due to the improvment of the precession constant value which makes the mean CEP equator of J2000.0 different from the reference plane of DE200. The value of the precession constant used for the reduction of the observations involved in DE200 was the IAU 1976 value, and, if we assume that the reference plane of DE200 results from the LLR observations only, the method leading to Eq. (4) yields here:
I 2000 I (MCEP)2000 (DE200) = -tp

at the time t of the mean epoch of LLR observations involved in DE200 (set here to 1 January 1975). The value of I I 2000 (MCEP)2000 (DE200) obtained is -0. 0859 ± 0. 0010 which differs from the value -0. 0799 ± 0. 0019 derived from Table 11 by about three times the greater error. If we now assume that the precession expressions of Lieske et al. (1977) have been used in the reduction of observations involved in DE200, an expression similar to Eq. (3), by changing Simon to Lieske and -0. 02440 to -0. 01896, gives for (MCEP) - (DE200) the value -4.74 mas which differs from the value -3.61 ± 0.13 mas, derived from Table 11, by about 9 . This comparison and the previous one yield less "good" results than the two internal comparisons above, and the large difference between the mean epochs of observations in DE200 and in this paper is probably a contributing factor. Nevertheless the formal errors concerning the line DE200 in Table 11 are probably smaller than the realistic ones by factors of 5 (Column 2) and 2 (Column 4). Table 12, derived from Table 11, gives the angles of the rotations transforming the ICRS axes to an "equatorial" celestial system of axes R with: u(R) = R1 (1 )R2 (2 )R3 (3 )u(ICRS) (u(R) and u(ICRS) are the vector of coordinates in R and ICRS, respectively), R standing for the MCEP system or for the reference axes of DE403 or DE200. 1, 3 , -2 are approximately the components, on the ICRS axes, of the unit vector pointing towards the origin of right ascensions in R and 2 , -1 , 1 are those of the unit vector of the polar axis of R. The quoted errors are derived from the 1- errors quoted in Table 11. For comparison, Table 12 mentions also the results obtained by Folkner et al. (1994) for the angles of the transformation of

the ICRS axes to the DE200 axes and for the angles 1 and 2 of the transformation of the ICRS axes to the MCEP system. For the latter, the values are closer to our values of the test fit 1 because they do not take into account Williams' correction quoted in Eq. (3). Table 12 gives also the values of 1 and 2 for the transformation of the ICRS axes to the MCEP system, mentioned in the IERS conventions 1996 (McCarthy, 1996). These values are comparable to our values from Sol. 1 but are certainely more accurate. It leads us to apply a factor of 2 to our 1- errors in the line MCEP of Table 12 in order to obtain realistic errors which insure the consistency of the two sets of values. This factor may also be applied to all the other lines of Table 12 which concern our results and to Columns 2 and 3 of Table 11, except for 1 and of DE200 which need a factor of 5. The 1- errors quoted in Column 4 of Table 11 are probably realistic. We note that our results concerning DE200 are consistent with those of Folkner et al. Our value of 1 from the test fit 1 differs from their value by about 5 times our estimated realistic error. From the results of Table 12, we derive the angles of the transformation rotating the reference axes of DE403 to the MCEP system, respectively 4.1 ± 0.3 mas, -15.6 ± 0.6 mas, -1.9 ± 1.8 mas (1- errors). Previous results have been given ´ in (Chapront & Chapront-Touze, 1997), but the present ones are much more precise for several reasons: the computation of the different corrections introduced in DT C , especially relativistic corrections and tropospheric effects, has been considerably improved, weights have been introduced in the fits, the libration solution has been improved, and the time span covered by LLR observations is larger. Our present values of the first two angles are in good agreement with the corresponding values (3.99 mas and -15.36 mas) adopted in the construction of DE403 (Standish et al., 1995). 7. Conclusion Beside values of the orbital and rotational parameters of the Moon, this paper gives the position of the inertial dynamical mean ecliptic of J2000.0 with respect to several celestial "equatorial" systems of axes, and derives from these results the transformations between the "equatorial" systems themselves. Several tests checking the validity of those results are performed


J. Chapront et al.: Determination of lunar parameters and ecliptic reference system orientation

633

and comparisons to the results of other authors are, for the most part, satisfying. This paper gives, in particular, the position of the inertial dynamical mean equinox of J2000.0, i.e. the ascending node of the inertial mean ecliptic of J2000.0 (the reference plane of the modern analytical theories for the Moon and the planets) on the mean CEP equator of J2000.0. It is shown that this equinox is closer than its rotational version, by a factor of 5, to the J2000.0 right ascension origin defined by polar motion and UT1 provided by IERS and by modern analytical theories of precession and nutation. The inertial dynamical mean equinox of J2000.0 is also closer to the projections, on the mean CEP equator of J2000.0, of the origins of right ascensions in the IERS celestial reference system and in the JPL numerical integration DE200. So we suggest that it replaces the rotational one in reference texts.
Acknowledgements. The authors are grateful to the LLR staff of CERGA for providing the observations used. They thank F. Mignard, director of CERGA, for fruitful discussions, M. Feissel, the former director of Central Bureau of IERS, for her encouragement to the part of this work which concerns reference systems, and D.D. McCarthy for helpful comments on the manuscript. They wish also to thank N. Capitaine, director of DANOF, who warmly welcomed them in her department.

References
Aoki S., Guinot B., Kaplan G.H., et al., 1982, A&A 105, 359 Boucher C., Altamimi Z., Feissel M., Sillard P., 1996, IERS Technical Note n 20, Observatoire de Paris Calame O., 1977, In: Mulholland J.D. (ed.) Scientific applications of Lunar Laser Ranging. Reidel, Dordrecht, p. 53

´ Chapront J., Chapront-Touze M., 1997, Celest. Mech. 66, 31 ´ Chapront J., Chapront-Touze M., Francou, G., 1998, submitted to Celest. Mech. ´ Chapront-Touze M., Chapront J., 1983, A&A 124, 50 ´ Chapront-Touze M., Chapront J., 1988, A&A 190, 342 Fairhead L., Bretagnon P., 1990, A&A 229, 240 Ferrari A.J., Sinclair W.S., Sjogren W.L., Williams J.G., Yoder C.F., 1980, JGR 85, 3939 Folkner W.M., Charlot P., Finger M.H., et al., 1994, A&A 287, 279 Herring T.A., 1991, In: Hughes J.A., Smith C.A., Kaplan G.H. (eds.) Proceedings of the 127th Colloquium of the IAU. USNO, Washington D.C., p. 157 Lieske J.H., Lederle T., Fricke W., Morando B., 1977, A&A 58, 1 Martin C.F., Torrence M.H., Misner L.W., 1985, JGR 90, 9403 McCarthy D.D. (ed.), 1992, IERS Technical Note n 13, Observatoire de Paris McCarthy D.D. (ed.), 1996, IERS Technical Note n 21, Observatoire de Paris Moons M., 1982, The Moon and the Planets 27, 257 Moons M., 1984, Celest. Mech. 34, 263 Newhall XX, Williams J.G., Dickey J.O., 1991, In: IERS Technical Note n 8, Observatoire de Paris Seidelmann P.K., 1982, Celest. Mech. 27, 79 Simon J.L., Bretagnon P., Chapront J., et al., 1994, A&A 282, 663 Standish E.M., 1981, A&A 101, L17 Standish E.M., 1982, A&A 114, 297 Standish E.M., Newhall XX, Williams J.G., Folkner W.M., 1995, JPL IOM 314.10-127 Williams J.G., 1994, AJ 108, 711 Williams J.G., Melbourne W.G., 1982, In: Calame O. (ed.) Highprecision Earth rotation and Earth-Moon dynamics. Reidel, Dordrecht, p. 283 Williams J.G., Newhall XX, Dickey J.O., 1996, Planet. Space Sci. 44, 1077 Zhu S.Y., Mueller I.I., 1983, Bull. Geod. 57, 29