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
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
|
SUBROUTINE sla_PRECL (EP0, EP1, RMATP)
*+
* - - - - - -
* P R E C L
* - - - - - -
*
* Form the matrix of precession between two epochs, using the
* model of Simon et al (1994), which is suitable for long
* periods of time.
*
* (double precision)
*
* Given:
* EP0 dp beginning epoch
* EP1 dp ending epoch
*
* Returned:
* RMATP dp(3,3) precession matrix
*
* Notes:
*
* 1) The epochs are TDB Julian epochs.
*
* 2) The matrix is in the sense V(EP1) = RMATP * V(EP0)
*
* 3) The absolute accuracy of the model is limited by the
* uncertainty in the general precession, about 0.3 arcsec per
* 1000 years. The remainder of the formulation provides a
* precision of 1 mas over the interval from 1000AD to 3000AD,
* 0.1 arcsec from 1000BC to 5000AD and 1 arcsec from
* 4000BC to 8000AD.
*
* Reference:
* Simon, J.L. et al., 1994. Astron.Astrophys., 282, 663-683.
*
* 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 (1000JY)
T0 = (EP0-2000D0)/1000D0
* Interval over which precession required (1000JY)
T = (EP1-EP0)/1000D0
* Euler angles
TAS2R = T*AS2R
W = 23060.9097D0+
: (139.7459D0+
: (-0.0038D0+
: (-0.5918D0+
: (-0.0037D0+
: 0.0007D0*T0)*T0)*T0)*T0)*T0
ZETA = (W+(30.2226D0+
: (-0.2523D0+
: (-0.3840D0+
: (-0.0014D0+
: 0.0007D0*T0)*T0)*T0)*T0+
: (18.0183D0+
: (-0.1326D0+
: (0.0006D0+
: 0.0005D0*T0)*T0)*T0+
: (-0.0583D0+
: (-0.0001D0+
: 0.0007D0*T0)*T0+
: (-0.0285D0+
: (-0.0002D0)*T)*T)*T)*T)*T)*TAS2R
Z = (W+(109.5270D0+
: (0.2446D0+
: (-1.3913D0+
: (-0.0134D0+
: 0.0026D0*T0)*T0)*T0)*T0+
: (18.2667D0+
: (-1.1400D0+
: (-0.0173D0+
: 0.0044D0*T0)*T0)*T0+
: (-0.2821D0+
: (-0.0093D0+
: 0.0032D0*T0)*T0+
: (-0.0301D0+
: 0.0006D0*T0
: -0.0001D0*T)*T)*T)*T)*T)*TAS2R
THETA = (20042.0207D0+
: (-85.3131D0+
: (-0.2111D0+
: (0.3642D0+
: (0.0008D0+
: (-0.0005D0)*T0)*T0)*T0)*T0)*T0+
: (-42.6566D0+
: (-0.2111D0+
: (0.5463D0+
: (0.0017D0+
: (-0.0012D0)*T0)*T0)*T0)*T0+
: (-41.8238D0+
: (0.0359D0+
: (0.0027D0+
: (-0.0001D0)*T0)*T0)*T0+
: (-0.0731D0+
: (0.0019D0+
: 0.0009D0*T0)*T0+
: (-0.0127D0+
: 0.0011D0*T0+0.0004D0*T)*T)*T)*T)*T)*TAS2R
* Rotation matrix
CALL sla_DEULER('ZYZ',-ZETA,THETA,-Z,RMATP)
END
|
