EDP Sciences
Free Access

This article has an erratum: [erratum]

Issue
A&A
Volume 534, October 2011
Article Number A22
Number of page(s) 19
Section Celestial mechanics and astrometry
DOI https://doi.org/10.1051/0004-6361/201117274
Published online 23 September 2011

© ESO, 2011

1. Introduction

Precession models are designed for two different phenomena: the precession of the ecliptic due to planetary perturbations and the precession of the equator due to the luni-solar and planetary torques on the oblate Earth. In both cases, precession represents the secular part of the motion. The term “secular” will be used throughout the paper to designate quasi periodic motions with very long periods. The motion of the Celestial Intermediate Pole (CIP), or equivalently of the equator of the CIP, with respect to the Geocentric Celestial Reference System (GCRS), is composed of precession and nutation, which are differentiated by a convention. Here we define the precession of the equator as that part of the motion of the equator that covers periods longer than 100 centuries, while terms of shorter periods are presumed to be included in the nutation.

In this connection, it is necessary to mention that the IAU 2000 model of nutation includes several terms with longer periods: 105 cy, 209 cy for the luni-solar terms and 933 cy, 150 cy, 129 cy, 113 cy for the planetary terms. The amplitudes of these terms are, however, very small (less than 4 mas for one term and less than 0.1 mas for the others), much smaller than the amplitudes of precession with similar periods.

Similarly, the osculating elements of the Earth-Moon barycenter (EMB) orbit are quasi-periodic functions of the time that can be expressed in the form of Poisson series whose arguments are linear combinations of the mean planetary longitudes. Here we define the precession of the ecliptic as being the secularly-moving ecliptic pole (i.e. mean EMB) orbital angular momentum) vector in a fixed ecliptic frame. The IAU 2006 precession of the ecliptic was computed as the part of the motion of the ecliptic covering periods longer than 300 centuries, while shorter ones are presumed to be included in the periodic component of the ecliptic motion. Almost all models of precession in use, including the most recent one, IAU 2006 (Capitaine et al. 2003; Hilton et al. 2006), are expressed in terms of polynomial developments of all the various precession parameters, which are intended for high-accuracy applications over a time span of a few centuries.

This paper is a continuation of our preceding studies (Vondrák et al. 2009, 2011), in which we used comparisons with numerical integrations to show that the IAU 2006 model is usable only for a limited time interval (depending on the parameterization used, as short as a few centuries around the epoch J2000.0), and its errors rapidly increase with longer time spans. An independent comparison by Fienga et al. (2008) of IAU 2006 with the numerical INPOP06 solution for the motion of the Earth’s axis over a few millennia shows discrepancies increasing from less than 0.2′′ to about 3′′ between  − 1 to  − 5 millennia from J2000.0. In reality, precession represents a rather complicated, very long-periodic process, whose periods are hundreds of centuries long. This behavior can clearly be seen in numerically integrated equations of motion of the Earth in the solar system and of its rotation.

The goal of the present study is to find relatively simple expressions for all precession parameters (listed, e.g., by Hilton et al. 2006), the primary ones being the orientation parameters of the secularly-moving ecliptic and equator poles with respect to a fixed celestial frame. We require that the accuracy of these expressions is comparable to the IAU 2006 model near the epoch J2000.0, while lower accuracy is allowed outside the interval  ± 1000 years, gradually increasing up to several arcminutes at the extreme epochs  ± 200 millennia. Thus the new expressions should be more universal, serving both for modern applications and a number of problems which require a precession model valid over long time intervals, in particular archaeoastronomy, e.g. predicting star alignments with ancient monuments, or calendrical studies involving the seasons.

2. Numerical representation of precession

2.1. Precession of the ecliptic

The numerical representation of the precession of the ecliptic was obtained as follows:

  • a)

    Integration of the solar system motion was performed using the Mercury 6 package (Chambers 1999), in the interval  ± 200 millennia from J2000.0 with a 1-day step. Mercury 6 is a general-purpose software package designed to calculate the orbital evolution of objects moving in the gravitational field of a large central body, such as the motion of the planets or the Earth-Moon barycenter in the solar system. The numerical integrator is based on a second-order mixed-variable symplectic (MVS) algorithm incorporating simple symplectic correctors. To compute the orbital evolution of the Earth-Moon barycenter we used the values of planetary masses and initial positions/velocities of the “big” bodies of the solar system as recommended by Chambers in Mercury 6. Only the Sun, Mercury, Venus, Earth-Moon barycenter, Mars, Jupiter, Saturn, Uranus, Neptune and Pluto were taken into account. Since the Mercury 6 package does not include general relativity, it was not applied. The comparison (see below) with the solution of Laskar et al. (1993), which includes relativity, demonstrated that it has no visible consequence on the computation of PA and QA. The integration was performed for the precession parameters PA = sinπAsinΠA, QA = sinπAcosΠA, which represent the orientation parameters of the secularly-moving ecliptic pole. The reference system in which these parameters are described is the mean equinox and ecliptic at J2000.0. The elements of the EMB orbit were then smoothed and interpolated with a 100-year step in order to retain only variations with periods longer than 3 millennia. The values of PA, QA were then checked against the results of Laskar (1993), who lists the values p = sini/2sinΩ, q = sini/2cosΩ in 1000-year steps from  − 20 to  + 10 Myr. We used the obvious relations , and interpolated the values, using a cubic spline function, at 100-year intervals. The comparison displays very good consistency between the two series – the differences are only a few milliarcseconds near the epoch J2000.0, reaching 20 arcseconds for the most distant epochs.

  • b)

    Integrated values for PA and QA were replaced, inside the interval  ± 1000 years around J2000.0, with the values computed from the IAU 2006 model for the precession of the ecliptic which, in turn, is based on the semi-analytical theory VSOP87 (Bretagnon & Francou 1988) and JPL DE406 (Standish 1998) ephemerides.

2.2. Precession of the equator with respect to the ecliptic of date

The precession of the equator was represented by the general precession in longitude, pA, and mean obliquity of date, εA, which are the orientation angles of the mean equator of date with respect to the ecliptic of date.

The numerical representation of the precession of the equator was obtained as follows:

  • a)

    The numerically integrated solution La93 by Laskaret al. (1993) ofpA and εA, available in the interval  ± 1 million years with a 1000-year step, was interpolated to give values in 100-year steps. The La93 solution for these precession quantities was obtained from a numerical integration of the precession equations (derived from the equations of the rigid Earth theory) based on the La93 precession of the ecliptic and on numerical values adopted at the reference epoch for the mean obliquity of the ecliptic, ϵ0, the speed of precession, and the geodesic precession. The expression for the general precession, pA, combines the precession in longitude of the equator and the precession of the ecliptic, the former being a function of the Earth’s dynamical flattening, of the masses of the Sun, the Earth and the Moon and other constants related to the orbital elements of the Moon and the EMB. The La93 solution for the quantity (εA − ϵ0) results from the combined effect of the general precession and precession of the ecliptic. Note that although an improved version of the La93 precession solution has been obtained by Laskar et al. (2004), intended for representing variations over several Myr, the differences between the La2004 and La93 solutions are insignificant in the interval  ± 200 millennia from J2000.0. As the form in which the La93 solution is available is the most suitable for the purpose of our work, this solution was preferred. Several corrections were applied to the La93 nominal solution in order to make it consistent with the constants and models upon which the IAU 2006 precession is based:

    • a linear correction of  −0.295767′′/cy to pA to account for the slightly different La93 value for the dynamical ellipticity of the Earth (HD = 0.00328005) compared with the IAU 2006 value (HD = 0.003273795), as well as for other small differences between the La93 and IAU 2006 dynamical models used in the precession equations;

    • a quadratic correction of  − 0.0071′′/cy2 to pA to account for effects considered in the IAU 2006 model and not in the La93 nominal solution, namely (i) the effect of the time rate of change, dJ2/dt =  −3.001 × 10-9 cy-1, in the dynamical form factor, (ii) the J2 and planetary tilt effects and (iii) the tidal effects, from Williams (1994);

    • a constant correction of  − 0.042′′ to εA to account for the slightly different value for the obliquity between the IAU 2006 and La93 values;

    • a linear correction of  − 0.02575′′/cy to εA, to account for the IAU 2006 value for the obliquity rate (i.e. the rate of precession in obliquity of the equator with respect to the ecliptic of epoch), which was zero in the La93 solution; the obliquity rate is due to the planetary tilt effect, the direct planetary effect and the tidal effect in obliquity (Williams 1994).

    The corrections just described can be summarized by Eq. (1); the numerical values of these total corrections for the dates of the numerical La93 solution were used to correct the latter. (1)The parameter T used in the above expression, as well as in those below, is the elapsed time in Julian centuries since J2000.0 TT, defined by: T = (TT − 2000  January  1d  12h  TT)/36525, with TT in days.

  • b)

    Inside the interval  ± 1000 years around J2000.0, the integrated values for pA and ϵA were replaced with the values computed from the IAU 2006 expressions of these quantities.

3. Calculation of the precession parameters

From the values of the precession parameters PA, QA, pA and ϵA, different precession parameters were calculated in the interval  ± 200 millennia from J2000.0.

The relations of the four above mentioned parameters to the other ones describing precession are shown in Fig. 1. To calculate them, we obtain first the auxiliary angles α,β,μ from the triangle ΥΥoN:

thumbnail Fig. 1

Precession parameters.

Open with DEXTER

(2)then the angles η, δ by solving the triangle ΥΥoPt:

(3)

and, from triangle ΥoPtPo, we get the precession parameters θA,ζA:

(4)

From the triangle PoPtCo the precession parameters ωA,ψA then follow:

(5)

and from the triangles PtCCo, PoPtCo the parameters χA,zA:

(6)

Instead of parameters θA,zA,ζA we use their combinations XA = sinθAcosζA, YA =  −sinθAsinζA, VA = sinθAsinzA, WA = sinθAcoszA that, unlike the former, are continuous functions of time. θA, ζA and zA exhibit rather large discontinuities (of about 94° for θA, 180° for ζA and zA) at intervals that are not quite regular: there is a change of sign approximately each 26 000 years. This would make the long-term analytical approximation of these parameters extremely difficult.

Note that XA and YA are the precession part of the coordinates X,Y of the CIP unit vector in the GCRS, or equivalently the coordinates of the secularly-moving CIP with respect to the mean equator and equinox at J2000.0, while WA and VA are the coordinates of the mean pole at J2000.0 with respect to the mean equator and equinox of date.

There are three more parameters ϕ, γ and ψ that can be calculated from the triangles PoCoC and PoCPt: (7)

We used Eqs. (2) through (7) to calculate time series of all the above precession parameters in the interval  ± 2000 cy with 1-cy steps. It is important to note that all these parameters are referred to the reference frame defined by the mean equator and equinox of standard epoch J2000.0, not to the slightly-offset GCRS (both pole and origin of the GCRS are not exactly identical with the pole and equinox of J2000.0; the offset, called “frame bias” is an unintended difference in orientation of about 23 mas between both frames – see subroutine ltp_PBMAT in the Appendix). Note also that the precession quantities are completely determined by the motions of the ecliptic pole and the equator pole, represented by the precession parameters PA,QA for the ecliptic and XA,YA for the equator. These four quantities will therefore be considered in the following as being the “primary precession parameters”.

4. Long term analytical approximations

4.1. General method

To find the long-term analytical approximation of each of the precession parameters, we performed the following steps:

  • spectral analysis of integrated values to find hiddenperiodicities, using the Vaníček (1969) method, based onleast-squares approximation, as modified by Vondrák (1977);

  • periods found were compared with those found by Laskar et al. (1993, 2004) and where there was a match Laskar’s values were adopted;

  • sine/cosine amplitudes of the terms found in the preceding step, plus  cubic polynomials, were fitted to the numerical integration, using a least-squares method. The weights used in the fit were chosen to avoid significant disagreement with IAU 2006 in its region of applicability while remaining consistent with the long-term character of the numerical integration: very high close to J2000.0 (equal to 104 inside the interval  ± 100 yr), decreasing quadratically with time (equal to 1/T2 outside the interval  ± 100 yr);

  • small additional corrections were then applied to the constant and linear terms to ensure that the function value and its first derivative at the epoch J2000.0 were identical with those of the IAU 2006 model.

The method used leads to results that are not unique; there are strong correlations between the pairs of estimated sine/cosine terms (often due to close periods that cannot be resolved in a given time interval), so that a tiny change of a period brings about rather large changes of all estimated amplitudes. However, the function values computed from different models within the interval  ± 200 millennia from J2000.0 are almost insensitive to these changes. One must keep in mind that the model is valid only within this interval, whereas outside the errors diverge rapidly. Here and in the following, these long-term expressions will be called “long-term model” or “new model”, with the understanding that the model is empirical, not physical.

Below we present long-term expressions for all the precession parameters. In these models, T is the time in Julian centuries running from J2000.0, as defined in Sect. 2.2. Periodic terms, whose cosine and sine amplitudes in arcseconds C,S and periods in centuries P are given in corresponding tables, have the general form  ∑ (Cicos2πT/Pi + Sisin2πT/Pi). In all these tables, according to Laskar et al. (1993), gi and si (i = 1,   6) denote the secular frequencies of the perihelions and nodes of the first six planets of the solar system, respectively, while p designates the main precession frequency, νk secular frequencies of the EMB orbital motion, and σk other frequencies of the La93 solution. Their values are provided in Tables 3 and 5 of Laskar et al. (2004) and in Table 2 of Laskar et al. (1993).

Comparisons of the new model with integrated values and the IAU 2006 model for all derived parameters are graphically depicted in Figs. 2 through 9. In all these figures, the new model is represented as a solid line, while the values corresponding to the numerical integration values described in Sect. 2 are represented as a dotted line. We note that, in all cases, the two curves are so close that they are graphically indistinguishable in the whole interval; they appear as a single line. On the other hand, the IAU 2006 model (dashed line) fits well both to the integrated values and to the new model near the epoch J2000.0, but it rapidly diverges from both for more distant epochs.

thumbnail Fig. 2

Long-term model of precession parameters PA,QA – new model and integrated values (solid), IAU2006 (dashed).

Open with DEXTER

4.2. Primary precession parameters

Table 1

Periodic terms in PA,QA.

The long-term expressions for the precession of the ecliptic, PA, QA, are given as (8)where the cosine/sine amplitudes of the periodic parts  ∑ P,  ∑ Q are given in Table 1. Names of some of the terms in Col. 1 come from Laskar (1993, 2004). The comparisons of the long-term models of the precession of the ecliptic, PA (top), QA (bottom) with integrated values and the IAU 2006 model are shown in Fig. 2.

thumbnail Fig. 3

Long-term model of precession parameters XA, YA – new model and integrated values (solid), IAU2006 (dashed).

Open with DEXTER

Table 2

Periodic terms in XA, YA.

The long-term expressions for the precession angles XA = sinθAcosζA, YA =  −sinθAsinζA, are given as (9)where the cosine/sine amplitudes of the periodic parts  ∑ X,  ∑ Y are given in Table 2. The comparisons of the long-term models of precession angles XA(top) and YA (bottom) are shown in Fig. 3.

As PA,QA and XA,YA represent direction cosines, it is easy to understand that Eqs. (8) and (9), having small coefficients in T and T2, are more appropriate expressions for representing those precession quantities in the long term than the usual IAU polynomial of times with large coefficients in T and T2.

thumbnail Fig. 4

Long-term model of precession parameters pA,εA – new model and integrated values (solid), IAU2006 (dashed).

Open with DEXTER

4.3. Other precession parameters

Table 3

Periodic terms in pA,εA.

Similarly, the long-term expressions for the general precession and obliquity, pA, εA, are given as (10)where the cosine/sine amplitudes of the periodic parts  ∑ p,  ∑ ε are given in Table 3. The comparisons of the long-term models of general precession pA (reduced by a conventional rate 5045′′/cy for clarity) and obliquity εA (bottom) are shown in Fig. 4.

thumbnail Fig. 5

Long-term model of precession parameters ψA,ωA – new model and integrated values (solid), IAU2006 (dashed).

Open with DEXTER

Long-term approximations of the precession angles ψA, ωA, are (11)where the cosine/sine amplitudes of the periodic parts  ∑ ψ,  ∑ ω are given in Table 4. The comparison of the long-term models with integrated and IAU 2006 values of precession angles ψA (top, again reduced by a conventional rate 5045′′/cy for clarity) and obliquity ωA (bottom) are shown in Fig. 5.

Table 4

Periodic terms in ψA, ωA.

Long-term expressions for the precession angles VA = sinθAsinzA, WA = sinθAcoszA, are given as (12)where the cosine/sine amplitudes of the periodic parts  ∑ V,  ∑ W are given in Table 5. The comparisons of the long-term models of precession angles VA (top) and WA (bottom) are shown in Fig. 6.

Table 5

Periodic terms in VA, WA.

thumbnail Fig. 6

Long-term model of precession parameters VA,WA – new model and integrated values (solid), IAU2006 (dashed).

Open with DEXTER

The parameter χA is approximated as (13)where the cosine/sine amplitudes of the periodic part  ∑ χ are given in Table 6. The comparison of the long-term model of precession angle χA is shown in Fig. 7.

Table 6

Periodic terms in χA.

thumbnail Fig. 7

Long-term model of precession parameter χA – new model and integrated values (solid), IAU2006 (dashed).

Open with DEXTER

The long-term expressions for precession parameters ϕ, γ are (14)where the cosine/sine amplitudes of the periodic parts  ∑ ϕ,  ∑ γ are given in Table 7. The comparisons of the long-term models of precession parameters ϕ and γ are illustrated in Fig. 8.

Table 7

Periodic terms in ϕ,γ.

thumbnail Fig. 8

Long-term model of precession parameters ϕ,γ – new model and integrated values (solid), IAU2006 (dashed).

Open with DEXTER

The last parameter ψ (not to be confused with a different parameter ψA, see Fig. 1) is approximated as (15)where the cosine/sine amplitudes of the periodic part  ∑ ψ are given in Table 8. The comparison of the long-term model of precession angle ψ, with a nominal reduction by 5045′′/cy, is shown in Fig. 9.

thumbnail Fig. 9

Long-term model of precession parameter ψ – new model and integrated values (solid), IAU2006 (dashed).

Open with DEXTER

Table 8

Periodic terms in ψ.

5. Precession matrix, and parameterization options

There are several ways to compute the precession matrix P, which is used to transform the rectangular coordinates of a celestial target from the coordinate system based upon the mean equator and equinox of epoch J2000.0 to that of an arbitrary date. We can use combinations of various parameters to do this; the different possibilities are discussed by Hilton et al. (2006). These possibilities include a form of the matrix that refers right ascensions to a point of the mean equator of date that is independent of the precession of the equinox (see Sect. 5.5).

5.1. Matrix based on the equatorial precession parameters

The classical parameterization of the precession matrix (Lieske et al. 1977) is obtained by a sequence of three rotations: by  − ζA (around the z-axis), θA (around the new position of the y-axis), and  − zA (around the new position of the z-axis): (16)where Ri stands for the rotation matrix around the ith axis. Thus the individual elements of P are given as (17)As mentioned earlier (Sect. 2) the properties of the three angles make the development of long-term series problematical. These difficulties are avoided by instead re-expressing the matrix elements in terms of the angles VA, WA, XA and YA: (18)where the developments of XA, YA and VA, WA are given by (9) and (12), respectively, and ZA = (1 − r2)1/2.

An obvious disadvantage of this parameterization is the singularity at the epoch J2000.0, when XA = YA = VA = WA = r2 = 0; in this special case, P is the identity matrix.

5.2. Matrix using the precession of the equator relative to the fixed ecliptic

The parameterization used e.g. by Capitaine et al. (2003) provides a clear separation between the precession of the equator and precession of the ecliptic, by using a sequence of four rotations: around the x-axis by ε0, around the z-axis by  − ψA, around the x-axis by  − ωA, and around the z-axis by χA (the first three describing the precession of the equator with respect to the fixed ecliptic, the last one precession of the ecliptic), i.e. (19)The angle ε0 = 84   381.406″; the remaining three parameters are given by developments (11) and (13). Elements of the precession matrix can then be computed from (20)

5.3. Matrix using the precession of the equator relative to the moving ecliptic

An alternative parameterization was proposed by Williams (1994) and Fukushima (2003). It is comprised of a sequence of four rotations: by γ around the z-axis, by ϕ around the x-axis, by  − ψ around the z-axis, and by  − εA around the x-axis (the first one describing the precession of the ecliptic, the last three precession of the equator relative to the moving ecliptic), i.e. (21)where the expressions to compute parameter εA are given by Eq. (10), ψ by Eq. (15), and ϕ, γ by Eqs. (14). The elements of matrix P can then be calculated as (22)

5.4. Matrix based on the ecliptic and mean equator poles

Fabri (1980) bypassed the conventional precession angles, and derived polynomials to deliver directly the two unit vectors representing the ecliptic and mean equator poles. Murray (1983, Sect. 5.4.2), used these vectors to generate the precession matrix: (23)where the vectors k and are the ecliptic pole and mean equator pole, respectively. In the Appendix we provide two computer subroutines, ltp_PECL and ltp_PEQU, to generate the two vectors (with respect to the J2000.0 equator and equinox), and in a third, ltp_PMAT, we use them to generate the precession matrix. A fourth subroutine, ltp_PBMAT, applies a small rotation (the “frame bias”) to the matrix, for use when the starting vector is with respect to the GCRS.

5.5. Matrix for the precession of the equator

In 2000/2003, the IAU introduced a new zero point for right ascension. The celestial intermediate origin (CIO) is a kinematically defined point that has no connection with equinox or ecliptic. Its use has the advantage that precession-nutation and Earth rotation are kept separate, and in particular means that sidereal time is replaced by Earth rotation angle (ERA), which is simply a linear transformation of UT1. It also means that precession-nutation is defined only by the motion of the celestial pole (cf. Sect. 6), which can be taken into account by the GCRS coordinates of the CIP unit vector, and consequently that we can construct a precession matrix using only the XA,YA series, Eq. (9). Using Capitaine & Wallace (2006), Eqs. (3)–(5), and neglecting nutation and bias, the rotation matrix from GCRS to current coordinates can be written as: (24)with (25)The quantity sA is derived from the kinematical definition of the CIO and involves an integral: Capitaine & Wallace (2006), Sect. 4, describe methods to do this. The development of sA derived from the expressions for XA and YA (Eq. (9)) must be taken into account for measuring the ERA in the long term. It is approximated by Eq. (26), in which  ∑ s is the periodic part, given in Table 9, (26)The comparison of the long-term model of locator sA with its IAU 2006 value, after a nominal reduction by  − 413′′/cy, is shown in Fig. 10. The complete value of s can be obtained by adding the effects of nutation, which consist of a small negative linear term plus a short-periodic part, the amplitude of which does not exceed a few arcseconds.

Table 9

Periodic terms in sA.

thumbnail Fig. 10

Locator sA due to precession – new model and integrated values (solid), IAU2006 (dashed).

Open with DEXTER

In contrast, only the second part of the matrix product in Eq. (24) is necessary to provide the celestial pole. Writing Eq. (24) as (27)defines Σ as the point on the mean equator such that ΥoN = ΣN, N being the node of the mean equator of date on the mean equator at J2000.0. PΣ can be expressed in terms of unit row vectors: (28)which identifies the top row of the PΣ matrix as the Σ vector, the bottom row as the mean equator pole vector and the middle row as the corresponding y-axis. As the point Σ is independent of the motion of the ecliptic, the matrix PΣ represents only the precession of the equator.

6. Accuracy comparisons

The ecliptic part of a precession model is essentially a matter of convention, and hence the position of the equinox is as well. In fact the position of the zero point of right ascension is irrelevant for many historical studies, because star visibilities and alignments (rise azimuth for example) depend only on the celestial pole. Consequently, most of the comparisons below (Sects. 6.1 to 6.4) are based solely on the angular separation between the given equator pole and that predicted by the XA,YA model (Eq. (9)), as implemented in the ltp_PEQU subroutine (Sect. A.2). However, for the comparison between the new precession model and reality (Sect. 6.5), we consider both the ecliptic and the equator.

6.1. New model versus numerical integration

Figure 11 summarizes how well the new model fits the numerical integrations (and IAU 2006 in the vicinity of J2000), based on the RMS accuracies for the individual parameters. The average RMS unit-weight error for all parameters was used, in combination with individual weights at different epochs, to compute the plot. The upper plot shows the full time span, the lower one holds for the central  ± 10 millennia.

thumbnail Fig. 11

Overall agreement between the new model and the numerical integration.

Open with DEXTER

thumbnail Fig. 12

New-model poles compared, using Eq. (9) as the reference model.

Open with DEXTER

6.2. Internal consistency of new model

The series given earlier provide not just one precession model, but a choice of several. Each relies on series that have been independently fitted to the numerical integrations, and so each can be expected to behave differently in some respects. It is reasonable to expect the different parameterizations to deliver the accuracies suggested by Fig. 11, but at the same time the statistical basis of each fit is the angle being fitted, not the consequences on the sky, and so differences in character can be expected.

Figure 12 compares the results from three parameterizations, with the ltp_PEQU subroutine (which uses the XA,YA parameterization) used as the reference model. The solid line is for the ε0,ψA,ωA,χA method. The dotted line is for the ζA,θA,zA method, these angles being deduced from VA,WA,XA,YA. The dashed line is for the γ,ϕ,ψ,ε0 method.

All of the parameterizations agree well within a few millennia of J2000.0, and are broadly consistent with the error estimates given above over most of the time span. At the most extreme epochs the γ,ϕ,ψ,ε0 and XA,YA methods agree less well. The sporadic improvements in the ζA,θA,zA case are probably because it and the reference method both use the XA,YA series as a basis.

6.3. Comparisons with existing formulations

Figure 13 compares the ltp_PMAT pole and several standard precession models, past and present, over a time interval that is slightly more than one precessional cycle. The dotted line shows the differences with respect to the IAU 1976 model, which uses the ζA,θA,zA parameterization. The finer of the two dashed lines is for the IAU 2000 model, using the ε0,ψA,ωA,χA parameterization, where ψA and ωA are from Lieske et al. 1977 with the IAU 2000 rate adjustments added. The solid line is the IAU 2006 model, implemented using the ε0,ψA,ωA,χA parameterization (and in fact giving almost exactly the same performance as IAU 1976 in the millennia of most historical interest). The coarser of the two dashed lines is also the IAU 2006 precession, but using direct series for CIP XA,YA, which are developments of the coordinates of the pole unit vector in the mean equator and equinox frame of epoch as polynomial of time around J2000.0; this method (denoted “IERS” in Fig. 13) combines precision and convenience in the modern era, but is not intended for remote epochs.

thumbnail Fig. 13

Comparison of Eq. (9) pole with other models, medium-term.

Open with DEXTER

thumbnail Fig. 14

Comparison of Eq. (9) pole with IAU 2006, short-term.

Open with DEXTER

6.4. New model versus IAU 2006, short term

Figure 14 compares the ltp_PMAT pole with the IAU 2006 model (via ε0,ψA,ωA,χA), for 2000 years centered on J2000. The central portion of the plot shows that replacing the precession portion of the IAU 2006/2000A CIP model would lead to changes in the 20th and 21st centuries that are rather less than 100   μas (comparable with current VLBI uncertainties) and may be acceptable in some multi-purpose computer applications.

6.5. New model versus reality

The evaluation of the accuracy of a model requires comparison with real observations. However, this is possible only for recent epochs, when precise observations are available. In the case of the precession models developed in this paper, the accuracy in the vicinity of J2000.0 has been ensured by the use, inside the interval  ± 1000 years around that epoch, of the IAU 2006 precession model, the basic constants of which were fitted to the DE406 ephemerides for the ecliptic and to VLBI observations for the equator. The long-term accuracy can be evaluated by comparison with other numerical integrations. This is possible for the motion of the ecliptic, for which independent integrated solutions exist (La93 and La2004 from Laskar et al. 1993; and 2004, respectively). The good consistency with the La93 solution for the ecliptic has been reported in Sect. 2.1. Note that, as the La2004 and La93 ecliptic solutions are nearly identical over  ± 200 millennia from J2000.0, this property is also valid with the La2004 solution.

For the equator, a comparison of the long term expression (Eq. (10)) for εA with the corresponding independent La2004 integrated solution shows discrepancies that are below 0.1′′ around the central J2000.0 epoch and reach about 200′′ for the most distant epochs. A similar comparison cannot be performed for the precession quantity pA, the La2004 solution of which is not available; however, an evaluation of the uncertainty can be reflected through the differences between two forms of the La93 solutions based on different tidal dissipation values, which appear to be within an order of magnitude of the discrepancies reported above for εA.

Another way to evaluate the accuracy is to evaluate the possible errors in the integration, i.e. both the accumulated errors in the long term integration and the expected uncertainties in the initial conditions and the numerical values of key parameters of the models upon which the solutions are based. The most important uncertainties for representing the EMB orbital angular momentum in the long term concern the physical parameters in the solar system and the knowledge of the model which is inevitably imperfect.

The most important uncertainties for the precession of the equator are in the change in the Earth’s dynamical ellipticity (or equivalently the dynamical form factor J2) and in the evolution of the dissipative effects, especially the tidal dissipation in the Earth-Moon system and the corresponding evolution of the orbit of the Moon. Note in particular that, even if the IAU 2006 values for dJ2/dt and the tidal dissipation are in agreement with the values derived from long term studies of the Earth rotation variations by Morrison & Stephenson (1997), based upon eclipse data over two millennia, the predictive knowledge of those values cannot be extended to a time span 200 times longer without increasing considerably the uncertainty.

7. Conclusions

The present study provides precession expressions that are appropriate for multi-purpose applications covering the  ± 200 millennia time span around J2000.0. These expressions are based on the IAU 2006 solutions close to J2000.0 and, for more distant epochs, a numerically integrated solution, made using the Mercury 6 package (Chambers 1999) to represent the ecliptic and La93 (Laskar et al. 1993) solution with upgraded constants, to represent general precession and obliquity. From them, different precession parameters have been calculated in the interval  ± 200 millennia from J2000.0, and analytical approximations have been found to obtain a good fit for the whole interval. The precession of the equator and the precession of the ecliptic have been clearly distinguished, by providing separately the expressions for the “primary precession parameters” (PA,QA for the ecliptic and XA,YA for the equator) that are sufficient to obtain the directions of the ecliptic pole and equator pole, respectively, in a fixed reference frame. Expressions for a number of other precession quantities have also been provided, to support various other kinds of calculations.

We have shown that these precession expressions, each of which comprises a cubic polynomial plus from 8 to 14 periodic terms, allow precession to be computed with accuracy comparable to IAU 2006 around the central epoch J2000.0, to a few arcseconds throughout the historical period, and to better than a degree at the ends of the  ± 200 millennia time span.

Software is provided that computes the ecliptic and mean equator poles, and the classical precession matrix.


Acknowledgments

This work was supported for one of us (J.V.) by the project LC506 “Recent dynamics of the Earth” of the Ministry of Education, Youth and Sports of the Czech Republic. We thank M. Šidlichovský for helping us with using Mercury 6 package, and also to the referee, D. D. McCarthy, for his valuable comments on the text.

References

Online material

Appendix:  Implementation in software

Here we present four computer subroutines that implement the long-term precession model. The chosen parameterization is PA,QA for the ecliptic pole and XA,YA for the equator pole. These results enable the precession matrix P to be computed. Two versions of the precession matrix are presented, the first for use with J2000.0 mean place and the second for use with Geocentric Celestial Reference System (GCRS) coordinates. Fortran code is used; implementations in other languages would be along very similar lines. The time argument for all four subroutines is Julian epoch (TT).

A.1. Precession of the ecliptic

The Fortran subroutine ltp_PECL generates the unit vector for the pole of the ecliptic, using the series for PA,QA (Eq. (8) and Table 1):

 SUBROUTINE ltp_PECL ( EPJ, VEC )
*+
* - - - - -
* P E C L
* - - - - -
*
* Long-term precession of the ecliptic.
*
* Given:
* EPJ d Julian epoch (TT)
*
* Returned:
* VEC d ecliptic pole unit vector
*
* The vector is with respect to the J2000.0 mean equator and equinox.
*
* Reference: Vondrak et al., A&A (2011), Eq.8, Tab.1
*
* Date: 2011 May 14
*
* Authors: J.Vondrak, N.Capitaine, P.Wallace
*
*-----------------------------------------------------------------------

 IMPLICIT NONE
 DOUBLE PRECISION EPJ, VEC(3)

* Arcseconds to radians
 DOUBLE PRECISION AS2R
 PARAMETER ( AS2R = 4.848136811095359935899141D-6 )

* 2Pi
 DOUBLE PRECISION D2PI
 PARAMETER ( D2PI = 6.283185307179586476925287D0 )

* Obliquity at J2000.0 (radians).
 DOUBLE PRECISION EPS0
 PARAMETER ( EPS0 = 84381.406D0 * AS2R )

* Number of polynomial terms
 INTEGER NPOL
 PARAMETER ( NPOL = 4 )

* Number of periodic terms
 INTEGER NPER
 PARAMETER ( NPER = 8 )

* Miscellaneous
 INTEGER I, J
 DOUBLE PRECISION T, P, Q, W, A, S, C, Z

* Polynomial and periodic coefficients
 DOUBLE PRECISION PQPOL(NPOL,2), PQPER(5,NPER)

* Polynomials
 DATA ((PQPOL(I,J),I=1,NPOL),J=1,2) /
 :
 : +5851.607687 D0,
 : -0.1189000 D0,
 : -0.00028913 D0,
 : +0.000000101 D0,
 :
 : -1600.886300 D0,
 : +1.1689818 D0,
 : -0.00000020 D0,
 : -0.000000437 D0 /

* Periodics
 DATA ((PQPER(I,J),I=1,5),J=1,NPER) /
 :
 : 708.15D0, -5486.751211D0, -684.661560D0,
 : 667.666730D0, -5523.863691D0,
 : 2309.00D0, -17.127623D0, 2446.283880D0,
 : -2354.886252D0, -549.747450D0,
 : 1620.00D0, -617.517403D0, 399.671049D0,
 : -428.152441D0, -310.998056D0,
 : 492.20D0, 413.442940D0, -356.652376D0,
 : 376.202861D0, 421.535876D0,
 : 1183.00D0, 78.614193D0, -186.387003D0,
 : 184.778874D0, -36.776172D0,
 : 622.00D0, -180.732815D0, -316.800070D0,
 : 335.321713D0, -145.278396D0,
 : 882.00D0, -87.676083D0, 198.296071D0,
 : -185.138669D0, -34.744450D0,
 : 547.00D0, 46.140315D0, 101.135679D0,
 : -120.972830D0, 22.885731D0 /


* Centuries since J2000.
 T = (EPJ-2000D0)/100D0

* Initialize P_A and Q_A accumulators.
 P = 0D0
 Q = 0D0

* Periodic terms.
 DO I=1,NPER
 W = D2PI*T
 A = W/PQPER(1,I)
 S = SIN(A)
 C = COS(A)
 P = P + C*PQPER(2,I) + S*PQPER(4,I)
 Q = Q + C*PQPER(3,I) + S*PQPER(5,I)
 END DO

* Polynomial terms.
 W = 1D0
 DO I=1,NPOL
 P = P + PQPOL(I,1)*W
 Q = Q + PQPOL(I,2)*W
 W = W*T
 END DO

* P_A and Q_A (radians).
 P = P*AS2R
 Q = Q*AS2R

* Form the ecliptic pole vector.
 Z = SQRT(MAX(1D0-P*P-Q*Q,0D0))
 S = SIN(EPS0)
 C = COS(EPS0)
 VEC(1) = P
 VEC(2) = - Q*C - Z*S
 VEC(3) = - Q*S + Z*C

 END

A.2. Precession of the equator

The Fortran subroutine ltp_PEQU generates the unit vector for the pole of the equator, using the series for XA,YA (Eq. (9) and Table 2):

 SUBROUTINE ltp_PEQU ( EPJ, VEQ )
*+
* - - - - -
* P E Q U
* - - - - -
*
* Long-term precession of the equator.
*
* Given:
* EPJ d Julian epoch (TT)
*
* Returned:
* VEQ d equator pole unit vector
*
* The vector is with respect to the J2000.0 mean equator and equinox.
*
* Reference: Vondrak et al., A&A (2011), Eq.9, Tab.2
*
* Date: 2011 May 14
*
* Authors: J.Vondrak, N.Capitaine, P.Wallace
*
*-----------------------------------------------------------------------

 IMPLICIT NONE
 DOUBLE PRECISION EPJ, VEQ(3)

* Arcseconds to radians
 DOUBLE PRECISION AS2R
 PARAMETER ( AS2R = 4.848136811095359935899141D-6 )

* 2Pi
 DOUBLE PRECISION D2PI
 PARAMETER ( D2PI = 6.283185307179586476925287D0 )

* Number of polynomial terms
 INTEGER NPOL
 PARAMETER ( NPOL = 4 )

* Number of periodic terms
 INTEGER NPER
 PARAMETER ( NPER = 14 )

* Miscellaneous
 INTEGER I, J
 DOUBLE PRECISION T, X, Y, W, A, S, C

* Polynomial and periodic coefficients
 DOUBLE PRECISION XYPOL(NPOL,2), XYPER(5,NPER)

* Polynomials
 DATA ((XYPOL(I,J),I=1,NPOL),J=1,2) /
 :
 : +5453.282155 D0,
 : +0.4252841 D0,
 : -0.00037173 D0,
 : -0.000000152 D0,
 :
 : -73750.930350 D0,
 : -0.7675452 D0,
 : -0.00018725 D0,
 : +0.000000231 D0 /

* Periodics
 DATA ((XYPER(I,J),I=1,5),J=1,NPER) /
 :
 : 256.75D0, -819.940624D0, 75004.344875D0,
 : 81491.287984D0, 1558.515853D0,
 : 708.15D0, -8444.676815D0, 624.033993D0,
 : 787.163481D0, 7774.939698D0,
 : 274.20D0, 2600.009459D0, 1251.136893D0,
 : 1251.296102D0, -2219.534038D0,
 : 241.45D0, 2755.175630D0, -1102.212834D0,
 : -1257.950837D0, -2523.969396D0,
 : 2309.00D0, -167.659835D0, -2660.664980D0,
 : -2966.799730D0, 247.850422D0,
 : 492.20D0, 871.855056D0, 699.291817D0,
 : 639.744522D0, -846.485643D0,
 : 396.10D0, 44.769698D0, 153.167220D0,
 : 131.600209D0, -1393.124055D0,
 : 288.90D0, -512.313065D0, -950.865637D0,
 : -445.040117D0, 368.526116D0,
 : 231.10D0, -819.415595D0, 499.754645D0,
 : 584.522874D0, 749.045012D0,
 : 1610.00D0, -538.071099D0, -145.188210D0,
 : -89.756563D0, 444.704518D0,
 : 620.00D0, -189.793622D0, 558.116553D0,
 : 524.429630D0, 235.934465D0,
 : 157.87D0, -402.922932D0, -23.923029D0,
 : -13.549067D0, 374.049623D0,
 : 220.30D0, 179.516345D0, -165.405086D0,
 : -210.157124D0, -171.330180D0,
 : 1200.00D0, -9.814756D0, 9.344131D0,
 : -44.919798D0, -22.899655D0 /


* Centuries since J2000.
 T = (EPJ-2000D0)/100D0

* Initialize X and Y accumulators.
 X = 0D0
 Y = 0D0

* Periodic terms.
 DO I=1,NPER
 W = D2PI*T
 A = W/XYPER(1,I)
 S = SIN(A)
 C = COS(A)
 X = X + C*XYPER(2,I) + S*XYPER(4,I)
 Y = Y + C*XYPER(3,I) + S*XYPER(5,I)
 END DO

* Polynomial terms.
 W = 1D0
 DO I=1,NPOL
 X = X + XYPOL(I,1)*W
 Y = Y + XYPOL(I,2)*W
 W = W*T
 END DO

* X and Y (direction cosines).
 X = X*AS2R
 Y = Y*AS2R

* Form the equator pole vector.
 VEQ(1) = X
 VEQ(2) = Y
 W = X*X + Y*Y
 IF ( W.LT.1D0 ) THEN
 VEQ(3) = SQRT(1D0-W)
 ELSE
 VEQ(3) = 0D0
 END IF

 END

A.3. Precession matrix, mean J2000.0

The Fortran subroutine ltp_PMAT generates the 3 × 3 rotation matrix P, constructed using Fabri parameterization (i.e. directly from the unit vectors for the ecliptic and equator poles – see Sect. 5.4). As well as calling the two previous subroutines, ltp_PMAT calls subroutines from the IAU SOFA library.1 The resulting matrix transforms vectors with respect to the mean equator and equinox of epoch 2000.0 into mean place of date.

 SUBROUTINE ltp_PMAT ( EPJ, RP )
*+
* - - - - -
* P M A T
* - - - - -
*
* Long-term precession matrix.
*
* Given:
* EPJ d Julian epoch (TT)
*
* Returned:
* RP d precession matrix, J2000.0 to date
*
* The matrix is in the sense
*
* P_date = RP x P_J2000,
*
* where P_J2000 is a vector with respect to the J2000.0 mean
* equator and equinox and P_date is the same vector with respect to
* the equator and equinox of epoch EPJ.
*
* Called:
* ltp_PEQU equator pole
* ltp_ECL ecliptic pole
* iau_PXP vector product (SOFA)
* iau_PN normalize vector (SOFA)
*
* Reference: Vondrak et al., A&A (2011), Eq.23
*
* Date: 2011 April 30
*
* Authors: J.Vondrak, N.Capitaine, P.Wallace
*
*-----------------------------------------------------------------------

 IMPLICIT NONE
 DOUBLE PRECISION EPJ, RP(3,3)

 DOUBLE PRECISION PEQR(3), PECL(3), V(3), W, EQX(3)

* Equator pole (bottom row of matrix).
 CALL ltp_PEQU ( EPJ, PEQR )

* Ecliptic pole.
 CALL ltp_PECL ( EPJ, PECL )

* Equinox (top row of matrix).
 CALL iau_PXP ( PEQR, PECL, V )
 CALL iau_PN ( V, W, EQX )

* Middle row of matrix.
 CALL iau_PXP ( PEQR, EQX, V )

* The matrix elements.
 RP(1,1) = EQX(1)
 RP(1,2) = EQX(2)
 RP(1,3) = EQX(3)
 RP(2,1) = V(1)
 RP(2,2) = V(2)
 RP(2,3) = V(3)
 RP(3,1) = PEQR(1)
 RP(3,2) = PEQR(2)
 RP(3,3) = PEQR(3)

 END

A.4. Precession matrix, GCRS

The Fortran subroutine ltp_PBMAT generates the 3 × 3 rotation matrix P × B, where B is the “frame bias matrix” that expresses the relative orientation of the GCRS and mean J2000.0 reference systems. A simple first-order implementation of the frame bias is used, the departure from rigor being well under 1   μas.

 SUBROUTINE ltp_PBMAT ( EPJ, RPB )
*+
* - - - - - -
* P B M A T
* - - - - - -
*
* Long-term precession matrix, including GCRS frame bias.
*
* Given:
* EPJ d Julian epoch (TT)
*
* Returned:
* RPB d precession-bias matrix, J2000.0 to date
*
* The matrix is in the sense
*
* P_date = RPB x P_J2000,
*
* where P_J2000 is a vector in the Geocentric Celestial Reference
* System, and P_date is the vector with respect to the Celestial
* Intermediate Reference System at that date but with nutation
* neglected.
*
* A first order bias formulation is used, of sub-microarcsecond
* accuracy compared with a full 3D rotation.
*
* Called:
* ltp_PMAT precession matrix
*
* Reference: Vondrak et al., A&A (2011), Section A.4.
*
* Date: 2011 April 30
*
* Authors: J.Vondrak, N.Capitaine, P.Wallace
*
*-----------------------------------------------------------------------

 IMPLICIT NONE
 DOUBLE PRECISION EPJ, RPB(3,3)

* Arcseconds to radians
 DOUBLE PRECISION AS2R
 PARAMETER ( AS2R = 4.848136811095359935899141D-6 )

* Frame bias (IERS Conventions 2010, Eqs. 5.21 and 5.33)
 DOUBLE PRECISION DX, DE, DR
 PARAMETER ( DX = -0.016617D0 * AS2R,
 : DE = -0.0068192D0 * AS2R,
 : DR = -0.0146D0 * AS2R )

 DOUBLE PRECISION RP(3,3)


* Precession matrix.
 CALL ltp_PMAT ( EPJ, RP )

* Apply the bias.
 RPB(1,1) = RP(1,1) - RP(1,2)*DR + RP(1,3)*DX
 RPB(1,2) = RP(1,1)*DR + RP(1,2) + RP(1,3)*DE
 RPB(1,3) = - RP(1,1)*DX - RP(1,2)*DE + RP(1,3)
 RPB(2,1) = RP(2,1) - RP(2,2)*DR + RP(2,3)*DX
 RPB(2,2) = RP(2,1)*DR + RP(2,2) + RP(2,3)*DE
 RPB(2,3) = - RP(2,1)*DX - RP(2,2)*DE + RP(2,3)
 RPB(3,1) = RP(3,1) - RP(3,2)*DR + RP(3,3)*DX
 RPB(3,2) = RP(3,1)*DR + RP(3,2) + RP(3,3)*DE
 RPB(3,3) = - RP(3,1)*DX - RP(3,2)*DE + RP(3,3)

 END

A.5. Test case

To assist in checking for correct use of the above subroutines, we present below the results of calling each of them for the following test date: –1374 (i.e. 1375 BCE) May 3 (Gregorian calendar) at 13:52:19.2 TT. The equivalent Julian date and Julian epoch are 1219339.078000 (TT) and  − 1373.5959534565 (TT) respectively.  

Calling ltp_PECL gives the following ecliptic pole (with respect to the J2000.0 mean equator and equinox): (A.1)Calling ltp_PEQU gives the following equator pole (with respect to the J2000.0 mean equator and equinox): (A.2)Calling ltp_PMAT gives the following precession matrix (mean equator and equinox of J2000.0 to mean equator and equinox of date): (A.3)Calling ltp_PECL gives the following bias + precession matrix (Geocentric Celestial Reference System to a precession-only variant of the Celestial Intermediate Reference System of date): (A.4)The above computations were performed using quadruple precision (REAL*128) so that all the reported decimals (except perhaps for the least significant digit in rare critical cases) are correct. Note that on typical computers ordinary double precision has only 15–16 decimal places of precision, and this must be taken into account when comparing results.

All Tables

Table 1

Periodic terms in PA,QA.

Table 2

Periodic terms in XA, YA.

Table 3

Periodic terms in pA,εA.

Table 4

Periodic terms in ψA, ωA.

Table 5

Periodic terms in VA, WA.

Table 6

Periodic terms in χA.

Table 7

Periodic terms in ϕ,γ.

Table 8

Periodic terms in ψ.

Table 9

Periodic terms in sA.

All Figures

thumbnail Fig. 1

Precession parameters.

Open with DEXTER
In the text
thumbnail Fig. 2

Long-term model of precession parameters PA,QA – new model and integrated values (solid), IAU2006 (dashed).

Open with DEXTER
In the text
thumbnail Fig. 3

Long-term model of precession parameters XA, YA – new model and integrated values (solid), IAU2006 (dashed).

Open with DEXTER
In the text
thumbnail Fig. 4

Long-term model of precession parameters pA,εA – new model and integrated values (solid), IAU2006 (dashed).

Open with DEXTER
In the text
thumbnail Fig. 5

Long-term model of precession parameters ψA,ωA – new model and integrated values (solid), IAU2006 (dashed).

Open with DEXTER
In the text
thumbnail Fig. 6

Long-term model of precession parameters VA,WA – new model and integrated values (solid), IAU2006 (dashed).

Open with DEXTER
In the text
thumbnail Fig. 7

Long-term model of precession parameter χA – new model and integrated values (solid), IAU2006 (dashed).

Open with DEXTER
In the text
thumbnail Fig. 8

Long-term model of precession parameters ϕ,γ – new model and integrated values (solid), IAU2006 (dashed).

Open with DEXTER
In the text
thumbnail Fig. 9

Long-term model of precession parameter ψ – new model and integrated values (solid), IAU2006 (dashed).

Open with DEXTER
In the text
thumbnail Fig. 10

Locator sA due to precession – new model and integrated values (solid), IAU2006 (dashed).

Open with DEXTER
In the text
thumbnail Fig. 11

Overall agreement between the new model and the numerical integration.

Open with DEXTER
In the text
thumbnail Fig. 12

New-model poles compared, using Eq. (9) as the reference model.

Open with DEXTER
In the text
thumbnail Fig. 13

Comparison of Eq. (9) pole with other models, medium-term.

Open with DEXTER
In the text
thumbnail Fig. 14

Comparison of Eq. (9) pole with IAU 2006, short-term.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.