Long-term precession of the equator.
Status: support routine.
The returned vector is with respect to the J2000.0 mean equator and equinox.
The Vondrak et al. (2011, 2012) 400 millennia precession model agrees with the IAU 2006 precession at J2000.0 and stays within 100 microarcseconds during the 20th and 21st centuries. It is accurate to a few arcseconds throughout the historical period, worsening to a few tenths of a degree at the end of the +/- 200,000 year time span.
Vondrak, J., Capitaine, N. and Wallace, P., 2011, New precession expressions, valid for long time intervals, Astron.Astrophys. 534, A22
Vondrak, J., Capitaine, N. and Wallace, P., 2012, New precession expressions, valid for long time intervals (Corrigendum), Astron.Astrophys. 541, C1
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=wp), | intent(in) | :: | epj | Julian epoch (TT) |
||
| real(kind=wp), | intent(out), | dimension(3) | :: | veq | equator pole unit vector |
subroutine LTPEQU ( epj, veq )
implicit none
real(wp),intent(in) :: epj !! Julian epoch (TT)
real(wp),dimension(3),intent(out) :: veq !! equator pole unit vector
! Number of polynomial terms
integer,parameter :: npol = 4
! Number of periodic terms
integer,parameter :: nper = 14
! Miscellaneous
integer :: i, j
real(wp) :: t, x, y, w, a, s, c
! Polynomial and periodic coefficients
! Polynomials
real(wp),dimension(npol,2),parameter :: xypol = reshape([&
+5453.282155_wp, &
+0.4252841_wp, &
-0.00037173_wp, &
-0.000000152_wp, &
-73750.930350_wp, &
-0.7675452_wp, &
-0.00018725_wp, &
+0.000000231_wp ], [npol,2])
! Periodics
real(wp),dimension(5,nper),parameter :: xyper = reshape([&
256.75_wp, -819.940624_wp, 75004.344875_wp, &
81491.287984_wp, 1558.515853_wp, &
708.15_wp, -8444.676815_wp, 624.033993_wp, &
787.163481_wp, 7774.939698_wp, &
274.20_wp, 2600.009459_wp, 1251.136893_wp, &
1251.296102_wp, -2219.534038_wp, &
241.45_wp, 2755.175630_wp, -1102.212834_wp, &
-1257.950837_wp, -2523.969396_wp, &
2309.00_wp, -167.659835_wp, -2660.664980_wp, &
-2966.799730_wp, 247.850422_wp, &
492.20_wp, 871.855056_wp, 699.291817_wp, &
639.744522_wp, -846.485643_wp, &
396.10_wp, 44.769698_wp, 153.167220_wp, &
131.600209_wp, -1393.124055_wp, &
288.90_wp, -512.313065_wp, -950.865637_wp, &
-445.040117_wp, 368.526116_wp, &
231.10_wp, -819.415595_wp, 499.754645_wp, &
584.522874_wp, 749.045012_wp, &
1610.00_wp, -538.071099_wp, -145.188210_wp, &
-89.756563_wp, 444.704518_wp, &
620.00_wp, -189.793622_wp, 558.116553_wp, &
524.429630_wp, 235.934465_wp, &
157.87_wp, -402.922932_wp, -23.923029_wp, &
-13.549067_wp, 374.049623_wp, &
220.30_wp, 179.516345_wp, -165.405086_wp, &
-210.157124_wp, -171.330180_wp, &
1200.00_wp, -9.814756_wp, 9.344131_wp, &
-44.919798_wp, -22.899655_wp ], [5,nper])
! Centuries since J2000.
t = (epj-2000.0_wp)/100.0_wp
! Initialize X and Y accumulators.
x = 0.0_wp
y = 0.0_wp
! Periodic terms.
w = d2pi*t
do i=1,nper
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 = 1.0_wp
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*das2r
y = y*das2r
! Form the equator pole vector.
veq(1) = x
veq(2) = y
veq(3) = sqrt(max(1.0_wp-x*x-y*y,0.0_wp))
end subroutine LTPEQU