1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
|
SUBROUTINE sla_PREC (EP0, EP1, RMATP)
*+
* - - - - -
* P R E C
* - - - - -
*
* Form the matrix of precession between two epochs (IAU 1976, FK5)
* (double precision)
*
* Given:
* EP0 dp beginning epoch
* EP1 dp ending epoch
*
* Returned:
* RMATP dp(3,3) precession matrix
*
* Notes:
*
* 1) The epochs are TDB (loosely ET) Julian epochs.
*
* 2) The matrix is in the sense V(EP1) = RMATP * V(EP0)
*
* 3) Though the matrix method itself is rigorous, the precession
* angles are expressed through canonical polynomials which are
* valid only for a limited time span. There are also known
* errors in the IAU precession rate. The absolute accuracy
* of the present formulation is better than 0.1 arcsec from
* 1960AD to 2040AD, better than 1 arcsec from 1640AD to 2360AD,
* and remains below 3 arcsec for the whole of the period
* 500BC to 3000AD. The errors exceed 10 arcsec outside the
* range 1200BC to 3900AD, exceed 100 arcsec outside 4200BC to
* 5600AD and exceed 1000 arcsec outside 6800BC to 8200AD.
* The SLALIB routine sla_PRECL implements a more elaborate
* model which is suitable for problems spanning several
* thousand years.
*
* References:
* Lieske,J.H., 1979. Astron.Astrophys.,73,282.
* equations (6) & (7), p283.
* Kaplan,G.H., 1981. USNO circular no. 163, pA2.
*
* Called: sla_DEULER
*
* P.T.Wallace Starlink 23 August 1996
*
* Copyright (C) 1996 Rutherford Appleton Laboratory
*-
IMPLICIT NONE
DOUBLE PRECISION EP0,EP1,RMATP(3,3)
* Arc seconds to radians
DOUBLE PRECISION AS2R
PARAMETER (AS2R=0.484813681109535994D-5)
DOUBLE PRECISION T0,T,TAS2R,W,ZETA,Z,THETA
* Interval between basic epoch J2000.0 and beginning epoch (JC)
T0 = (EP0-2000D0)/100D0
* Interval over which precession required (JC)
T = (EP1-EP0)/100D0
* Euler angles
TAS2R = T*AS2R
W = 2306.2181D0+(1.39656D0-0.000139D0*T0)*T0
ZETA = (W+((0.30188D0-0.000344D0*T0)+0.017998D0*T)*T)*TAS2R
Z = (W+((1.09468D0+0.000066D0*T0)+0.018203D0*T)*T)*TAS2R
THETA = ((2004.3109D0+(-0.85330D0-0.000217D0*T0)*T0)
: +((-0.42665D0-0.000217D0*T0)-0.041833D0*T)*T)*TAS2R
* Rotation matrix
CALL sla_DEULER('ZYZ',-ZETA,THETA,-Z,RMATP)
END
|
