New precession expressions, valid for long time intervals^{⋆}
^{1} Astronomical Institute, Acad. Sci. Czech Rep., Boční II, 141 31 Prague 4, Czech Republic
email: vondrak@ig.cas.cz
^{2}
SYRTE, Observatoire de Paris, CNRS, UPMC, 61 avenue de l’Observatoire, 75014 Paris, France
email: n.capitaine@obspm.fr
^{3}
STFC/Rutherford Appleton Laboratory, Harwell Oxford, Oxon OX11 0QX, UK
email: patrick.wallace@stfc.ac.uk
Received: 17 May 2011
Accepted: 28 July 2011
Context. The present IAU model of precession, like its predecessors, is given as a set of polynomial approximations of various precession parameters intended for highaccuracy applications over a limited time span. Earlier comparisons with numerical integrations have shown that this model is valid only for a few centuries around the basic epoch, J2000.0, while for more distant epochs it rapidly diverges from the numerical solution. In our preceding studies we also obtained preliminary developments for the precessional contribution to the motion of the equator: coordinates X,Y of the precessing pole and precession parameters ψ_{A},ω_{A}, suitable for use over long time intervals.
Aims. The goal of the present paper is to obtain upgraded developments for various sets of precession angles that would fit modern observations near J2000.0 and at the same time fit numerical integration of the motions of solar system bodies on scales of several thousand centuries.
Methods. We used the IAU 2006 solutions to represent the precession of the ecliptic and of the equator close to J2000.0 and, for more distant epochs, a numerical integration using the Mercury 6 package and solutions by Laskar et al. (1993, A&A, 270, 522) with upgraded initial conditions and constants to represent the ecliptic, and general precession and obliquity, respectively. From them, different precession parameters were calculated in the interval ± 200 millennia from J2000.0, and analytical expressions are found that provide a good fit for the whole interval.
Results. Series for the various precessional parameters, comprising a cubic polynomial plus from 8 to 14 periodic terms, are derived that allow precession to be computed with an accuracy comparable to IAU 2006 around the central epoch J2000.0, a few arcseconds throughout the historical period, and a few tenths of a degree at the ends of the ± 200 millennia time span. Computer algorithms are provided that compute the ecliptic and mean equator poles and the precession matrix.
Key words: astrometry / ephemerides / reference systems
The Appendix containing the computer code is available in electronic form at http://www.aanda.org
© 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 lunisolar 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 lunisolar 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 EarthMoon barycenter (EMB) orbit are quasiperiodic 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 secularlymoving 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 highaccuracy 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 longperiodic 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 secularlymoving 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 1day step. Mercury 6 is a generalpurpose 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 EarthMoon barycenter in the solar system. The numerical integrator is based on a secondorder mixedvariable symplectic (MVS) algorithm incorporating simple symplectic correctors. To compute the orbital evolution of the EarthMoon 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, EarthMoon 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 P_{A} and Q_{A}. The integration was performed for the precession parameters P_{A} = sinπ_{A}sinΠ_{A}, Q_{A} = sinπ_{A}cosΠ_{A}, which represent the orientation parameters of the secularlymoving 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 100year step in order to retain only variations with periods longer than 3 millennia. The values of P_{A}, Q_{A} were then checked against the results of Laskar (1993), who lists the values p = sini/2sinΩ, q = sini/2cosΩ in 1000year steps from − 20 to + 10 Myr. We used the obvious relations , and interpolated the values, using a cubic spline function, at 100year 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 P_{A} and Q_{A} 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 semianalytical 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, p_{A}, 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) ofp_{A} and ε_{A}, available in the interval ± 1 million years with a 1000year step, was interpolated to give values in 100year 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, p_{A}, 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 p_{A} to account for the slightly different La93 value for the dynamical ellipticity of the Earth (H_{D} = 0.00328005) compared with the IAU 2006 value (H_{D} = 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′′/cy^{2} to p_{A} 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, dJ_{2}/dt = −3.001 × 10^{9} cy^{1}, in the dynamical form factor, (ii) the J_{2} 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).


b)
Inside the interval ± 1000 years around J2000.0, the integrated values for p_{A} 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 P_{A}, Q_{A}, p_{A} 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 ΥΥ_{o}N:
Fig. 1 Precession parameters. 

Open with DEXTER 
(2)then the angles η, δ by solving the triangle ΥΥ_{o}P_{t}:
and, from triangle Υ_{o}P_{t}P_{o}, we get the precession parameters θ_{A},ζ_{A}:
From the triangle P_{o}P_{t}C_{o} the precession parameters ω_{A},ψ_{A} then follow:
and from the triangles P_{t}CC_{o}, P_{o}P_{t}C_{o} the parameters χ_{A},z_{A}:
Instead of parameters θ_{A},z_{A},ζ_{A} we use their combinations X_{A} = sinθ_{A}cosζ_{A}, Y_{A} = −sinθ_{A}sinζ_{A}, V_{A} = sinθ_{A}sinz_{A}, W_{A} = sinθ_{A}cosz_{A} that, unlike the former, are continuous functions of time. θ_{A}, ζ_{A} and z_{A} exhibit rather large discontinuities (of about 94° for θ_{A}, 180° for ζ_{A} and z_{A}) at intervals that are not quite regular: there is a change of sign approximately each 26 000 years. This would make the longterm analytical approximation of these parameters extremely difficult.
Note that X_{A} and Y_{A} are the precession part of the coordinates X,Y of the CIP unit vector in the GCRS, or equivalently the coordinates of the secularlymoving CIP with respect to the mean equator and equinox at J2000.0, while W_{A} and V_{A} 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 P_{o}C_{o}C and P_{o}CP_{t}: (7)
We used Eqs. (2) through (7) to calculate time series of all the above precession parameters in the interval ± 2000 cy with 1cy 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 slightlyoffset 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 P_{A},Q_{A} for the ecliptic and X_{A},Y_{A} 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 longterm 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 onleastsquares 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 leastsquares 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 longterm character of the numerical integration: very high close to J2000.0 (equal to 10^{4} inside the interval ± 100 yr), decreasing quadratically with time (equal to 1/T^{2} 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 longterm expressions will be called “longterm model” or “new model”, with the understanding that the model is empirical, not physical.
Below we present longterm 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 ∑ (C_{i}cos2πT/P_{i} + S_{i}sin2πT/P_{i}). In all these tables, according to Laskar et al. (1993), g_{i} and s_{i} (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.
Fig. 2 Longterm model of precession parameters P_{A},Q_{A} – new model and integrated values (solid), IAU2006 (dashed). 

Open with DEXTER 
4.2. Primary precession parameters
Periodic terms in P_{A},Q_{A}.
The longterm expressions for the precession of the ecliptic, P_{A}, Q_{A}, 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 longterm models of the precession of the ecliptic, P_{A} (top), Q_{A} (bottom) with integrated values and the IAU 2006 model are shown in Fig. 2.
Fig. 3 Longterm model of precession parameters X_{A}, Y_{A} – new model and integrated values (solid), IAU2006 (dashed). 

Open with DEXTER 
Periodic terms in X_{A}, Y_{A}.
The longterm expressions for the precession angles X_{A} = sinθ_{A}cosζ_{A}, Y_{A} = −sinθ_{A}sinζ_{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 longterm models of precession angles X_{A}(top) and Y_{A} (bottom) are shown in Fig. 3.
As P_{A},Q_{A} and X_{A},Y_{A} represent direction cosines, it is easy to understand that Eqs. (8) and (9), having small coefficients in T and T^{2}, 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 T^{2}.
Fig. 4 Longterm model of precession parameters p_{A},ε_{A} – new model and integrated values (solid), IAU2006 (dashed). 

Open with DEXTER 
4.3. Other precession parameters
Periodic terms in p_{A},ε_{A}.
Similarly, the longterm expressions for the general precession and obliquity, p_{A}, ε_{A}, are given as (10)where the cosine/sine amplitudes of the periodic parts ∑ _{p}, ∑ _{ε} are given in Table 3. The comparisons of the longterm models of general precession p_{A} (reduced by a conventional rate 5045′′/cy for clarity) and obliquity ε_{A} (bottom) are shown in Fig. 4.
Fig. 5 Longterm model of precession parameters ψ_{A},ω_{A} – new model and integrated values (solid), IAU2006 (dashed). 

Open with DEXTER 
Longterm 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 longterm 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.
Periodic terms in ψ_{A}, ω_{A}.
Longterm expressions for the precession angles V_{A} = sinθ_{A}sinz_{A}, W_{A} = sinθ_{A}cosz_{A}, are given as (12)where the cosine/sine amplitudes of the periodic parts ∑ _{V}, ∑ _{W} are given in Table 5. The comparisons of the longterm models of precession angles V_{A} (top) and W_{A} (bottom) are shown in Fig. 6.
Periodic terms in V_{A}, W_{A}.
Fig. 6 Longterm model of precession parameters V_{A},W_{A} – 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 longterm model of precession angle χ_{A} is shown in Fig. 7.
Periodic terms in χ_{A}.
Fig. 7 Longterm model of precession parameter χ_{A} – new model and integrated values (solid), IAU2006 (dashed). 

Open with DEXTER 
The longterm expressions for precession parameters ϕ, γ are (14)where the cosine/sine amplitudes of the periodic parts ∑ _{ϕ}, ∑ _{γ} are given in Table 7. The comparisons of the longterm models of precession parameters ϕ and γ are illustrated in Fig. 8.
Periodic terms in ϕ,γ.
Fig. 8 Longterm 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 longterm model of precession angle ψ, with a nominal reduction by 5045′′/cy, is shown in Fig. 9.
Fig. 9 Longterm model of precession parameter ψ – new model and integrated values (solid), IAU2006 (dashed). 

Open with DEXTER 
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 zaxis), θ_{A} (around the new position of the yaxis), and − z_{A} (around the new position of the zaxis): (16)where R_{i} 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 longterm series problematical. These difficulties are avoided by instead reexpressing the matrix elements in terms of the angles V_{A}, W_{A}, X_{A} and Y_{A}: (18)where the developments of X_{A}, Y_{A} and V_{A}, W_{A} are given by (9) and (12), respectively, and Z_{A} = (1 − r^{2})^{1/2}.
An obvious disadvantage of this parameterization is the singularity at the epoch J2000.0, when X_{A} = Y_{A} = V_{A} = W_{A} = r^{2} = 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 xaxis by ε_{0}, around the zaxis by − ψ_{A}, around the xaxis by − ω_{A}, and around the zaxis 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 zaxis, by ϕ around the xaxis, by − ψ around the zaxis, and by − ε_{A} around the xaxis (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 precessionnutation 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 precessionnutation 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 X_{A},Y_{A} 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 s_{A} 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 s_{A} derived from the expressions for X_{A} and Y_{A} (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 longterm model of locator s_{A} 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 shortperiodic part, the amplitude of which does not exceed a few arcseconds.
Periodic terms in s_{A}.
Fig. 10 Locator s_{A} 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 Υ_{o}N = Σ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 yaxis. 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 X_{A},Y_{A} 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 unitweight 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.
Fig. 11 Overall agreement between the new model and the numerical integration. 

Open with DEXTER 
Fig. 12 Newmodel 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 X_{A},Y_{A} 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},z_{A} method, these angles being deduced from V_{A},W_{A},X_{A},Y_{A}. 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 X_{A},Y_{A} methods agree less well. The sporadic improvements in the ζ_{A},θ_{A},z_{A} case are probably because it and the reference method both use the X_{A},Y_{A} 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},z_{A} 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 X_{A},Y_{A}, 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.
Fig. 13 Comparison of Eq. (9) pole with other models, mediumterm. 

Open with DEXTER 
Fig. 14 Comparison of Eq. (9) pole with IAU 2006, shortterm. 

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 multipurpose 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 longterm 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 p_{A}, 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 J_{2}) and in the evolution of the dissipative effects, especially the tidal dissipation in the EarthMoon system and the corresponding evolution of the orbit of the Moon. Note in particular that, even if the IAU 2006 values for dJ_{2}/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 multipurpose 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” (P_{A},Q_{A} for the ecliptic and X_{A},Y_{A} 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.
See www.iausofa.org.
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
 Bretagnon, P., & Francou, G. 1988, A&A, 202, 309 [NASA ADS] [Google Scholar]
 Capitaine, N., & Wallace, P. T. 2006, A&A, 450, 855 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Capitaine, N., Wallace, P. T., & Chapront, J. 2003, A&A, 412, 567 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chambers, J. E. 1999, MNRAS, 304, 793 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Fabri, E. 1980, A&A, 82, 123 [NASA ADS] [Google Scholar]
 Fienga, A., Manche, M., Laskar, H., & Gastineau, J. 2008, A&A, 477, 315 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fukushima, T. 2003, AJ, 126, 494 [NASA ADS] [CrossRef] [Google Scholar]
 Hilton, J., Capitaine, N., Chapront, J., et al. 2006, Cel. Mech. Dyn. Astr., 94, 351 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J., Joutel, F., & Boudin, F. 1993, A&A, 270, 522 [NASA ADS] [Google Scholar]
 Laskar, J., Robutel, P., Joutel, F., et al. 2004, A&A, 428, 261 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lieske, J. H., Lederle, T., Fricke, W., & Morando, B. 1977, A&A, 58, 1 [NASA ADS] [Google Scholar]
 Morrison, L., & Stephenson, R. 1997, Contemp. Phys., 38, 13 [NASA ADS] [CrossRef] [Google Scholar]
 Murray, C. A. 1983, Vectorial Astrometry, Adam Hilger Ltd, Bristol, ISBN 0852743726 [Google Scholar]
 Standish, E. M. 1998, in JPL IOM 312, F98048 [Google Scholar]
 Vaníček, P. 1969, Ap&SS, 4, 387 [NASA ADS] [CrossRef] [Google Scholar]
 Vondrák, J. 1977, Studia Geopys. Geod., 21, 107 [NASA ADS] [CrossRef] [Google Scholar]
 Vondrák, J., Capitaine, N., & Wallace, P. T. 2009, in Journées 2008 Systèmes de référence spatiotemporels, ed. M. Soffel, & N. Capitaine, Lohrmann Observatorium Dresden and Observatoire de Paris, 23 [Google Scholar]
 Vondrák, J., Capitaine, N., & Wallace, P. T. 2011, in Journées 2010 Systèmes de référence spatiotemporels, ed. N. Capitaine, Observatoire de Paris, in press [Google Scholar]
 Williams, J. G. 1994, AJ, 108, 711 [NASA ADS] [CrossRef] [Google Scholar]
Online material
Appendix: Implementation in software
Here we present four computer subroutines that implement the longterm precession model. The chosen parameterization is P_{A},Q_{A} for the ecliptic pole and X_{A},Y_{A} 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 P_{A},Q_{A} (Eq. (8) and Table 1):
SUBROUTINE ltp_PECL ( EPJ, VEC ) *+ *      * P E C L *      * * Longterm 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.848136811095359935899141D6 ) * 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 = (EPJ2000D0)/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(1D0P*PQ*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 X_{A},Y_{A} (Eq. (9) and Table 2):
SUBROUTINE ltp_PEQU ( EPJ, VEQ ) *+ *      * P E Q U *      * * Longterm 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.848136811095359935899141D6 ) * 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 = (EPJ2000D0)/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(1D0W) 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 *      * * Longterm 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 firstorder 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 *       * * Longterm precession matrix, including GCRS frame bias. * * Given: * EPJ d Julian epoch (TT) * * Returned: * RPB d precessionbias 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 submicroarcsecond * 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.848136811095359935899141D6 ) * 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 precessiononly 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
All Figures
Fig. 1 Precession parameters. 

Open with DEXTER  
In the text 
Fig. 2 Longterm model of precession parameters P_{A},Q_{A} – new model and integrated values (solid), IAU2006 (dashed). 

Open with DEXTER  
In the text 
Fig. 3 Longterm model of precession parameters X_{A}, Y_{A} – new model and integrated values (solid), IAU2006 (dashed). 

Open with DEXTER  
In the text 
Fig. 4 Longterm model of precession parameters p_{A},ε_{A} – new model and integrated values (solid), IAU2006 (dashed). 

Open with DEXTER  
In the text 
Fig. 5 Longterm model of precession parameters ψ_{A},ω_{A} – new model and integrated values (solid), IAU2006 (dashed). 

Open with DEXTER  
In the text 
Fig. 6 Longterm model of precession parameters V_{A},W_{A} – new model and integrated values (solid), IAU2006 (dashed). 

Open with DEXTER  
In the text 
Fig. 7 Longterm model of precession parameter χ_{A} – new model and integrated values (solid), IAU2006 (dashed). 

Open with DEXTER  
In the text 
Fig. 8 Longterm model of precession parameters ϕ,γ – new model and integrated values (solid), IAU2006 (dashed). 

Open with DEXTER  
In the text 
Fig. 9 Longterm model of precession parameter ψ – new model and integrated values (solid), IAU2006 (dashed). 

Open with DEXTER  
In the text 
Fig. 10 Locator s_{A} due to precession – new model and integrated values (solid), IAU2006 (dashed). 

Open with DEXTER  
In the text 
Fig. 11 Overall agreement between the new model and the numerical integration. 

Open with DEXTER  
In the text 
Fig. 12 Newmodel poles compared, using Eq. (9) as the reference model. 

Open with DEXTER  
In the text 
Fig. 13 Comparison of Eq. (9) pole with other models, mediumterm. 

Open with DEXTER  
In the text 
Fig. 14 Comparison of Eq. (9) pole with IAU 2006, shortterm. 

Open with DEXTER  
In the text 