Online material

Appendix:  Implementation in software

Here we present four computer subroutines that implement the long-term 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
* - - - - -
*
* 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 X_A,Y_A (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.¹ 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.

© ESO, 2011