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
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
|
SUBROUTINE slPLNE (DATE, JFORM, EPOCH, ORBINC, ANODE, PERIH,
: AORQ, E, AORL, DM, PV, JSTAT)
*+
* - - - - - - -
* P L N L
* - - - - - - -
*
* Heliocentric position and velocity of a planet, asteroid or comet,
* starting from orbital elements.
*
* Given:
* DATE d date, Modified Julian Date (JD - 2400000.5, Note 1)
* JFORM i choice of element set (1-3; Note 3)
* EPOCH d epoch of elements (TT MJD, Note 4)
* ORBINC d inclination (radians)
* ANODE d longitude of the ascending node (radians)
* PERIH d longitude or argument of perihelion (radians)
* AORQ d mean distance or perihelion distance (AU)
* E d eccentricity
* AORL d mean anomaly or longitude (radians, JFORM=1,2 only)
* DM d daily motion (radians, JFORM=1 only)
*
* Returned:
* PV d(6) heliocentric x,y,z,xdot,ydot,zdot of date,
* J2000 equatorial triad (AU,AU/s)
* JSTAT i status: 0 = OK
* -1 = illegal JFORM
* -2 = illegal E
* -3 = illegal AORQ
* -4 = illegal DM
* -5 = numerical error
*
* Called: slELUE, slUEPV
*
* Notes
*
* 1 DATE is the instant for which the prediction is required. It is
* in the TT timescale (formerly Ephemeris Time, ET) and is a
* Modified Julian Date (JD-2400000.5).
*
* 2 The elements are with respect to the J2000 ecliptic and equinox.
*
* 3 A choice of three different element-set options is available:
*
* Option JFORM = 1, suitable for the major planets:
*
* EPOCH = epoch of elements (TT MJD)
* ORBINC = inclination i (radians)
* ANODE = longitude of the ascending node, big omega (radians)
* PERIH = longitude of perihelion, curly pi (radians)
* AORQ = mean distance, a (AU)
* E = eccentricity, e (range 0 to <1)
* AORL = mean longitude L (radians)
* DM = daily motion (radians)
*
* Option JFORM = 2, suitable for minor planets:
*
* EPOCH = epoch of elements (TT MJD)
* ORBINC = inclination i (radians)
* ANODE = longitude of the ascending node, big omega (radians)
* PERIH = argument of perihelion, little omega (radians)
* AORQ = mean distance, a (AU)
* E = eccentricity, e (range 0 to <1)
* AORL = mean anomaly M (radians)
*
* Option JFORM = 3, suitable for comets:
*
* EPOCH = epoch of elements and perihelion (TT MJD)
* ORBINC = inclination i (radians)
* ANODE = longitude of the ascending node, big omega (radians)
* PERIH = argument of perihelion, little omega (radians)
* AORQ = perihelion distance, q (AU)
* E = eccentricity, e (range 0 to 10)
*
* Unused arguments (DM for JFORM=2, AORL and DM for JFORM=3) are not
* accessed.
*
* 4 Each of the three element sets defines an unperturbed heliocentric
* orbit. For a given epoch of observation, the position of the body
* in its orbit can be predicted from these elements, which are
* called "osculating elements", using standard two-body analytical
* solutions. However, due to planetary perturbations, a given set
* of osculating elements remains usable for only as long as the
* unperturbed orbit that it describes is an adequate approximation
* to reality. Attached to such a set of elements is a date called
* the "osculating epoch", at which the elements are, momentarily,
* a perfect representation of the instantaneous position and
* velocity of the body.
*
* Therefore, for any given problem there are up to three different
* epochs in play, and it is vital to distinguish clearly between
* them:
*
* . The epoch of observation: the moment in time for which the
* position of the body is to be predicted.
*
* . The epoch defining the position of the body: the moment in time
* at which, in the absence of purturbations, the specified
* position (mean longitude, mean anomaly, or perihelion) is
* reached.
*
* . The osculating epoch: the moment in time at which the given
* elements are correct.
*
* For the major-planet and minor-planet cases it is usual to make
* the epoch that defines the position of the body the same as the
* epoch of osculation. Thus, only two different epochs are
* involved: the epoch of the elements and the epoch of observation.
*
* For comets, the epoch of perihelion fixes the position in the
* orbit and in general a different epoch of osculation will be
* chosen. Thus, all three types of epoch are involved.
*
* For the present routine:
*
* . The epoch of observation is the argument DATE.
*
* . The epoch defining the position of the body is the argument
* EPOCH.
*
* . The osculating epoch is not used and is assumed to be close
* enough to the epoch of observation to deliver adequate accuracy.
* If not, a preliminary call to slPRTL may be used to update
* the element-set (and its associated osculating epoch) by
* applying planetary perturbations.
*
* 5 The reference frame for the result is with respect to the mean
* equator and equinox of epoch J2000.
*
* 6 The algorithm was originally adapted from the EPHSLA program of
* D.H.P.Jones (private communication, 1996). The method is based
* on Stumpff's Universal Variables.
*
* Reference: Everhart, E. & Pitkin, E.T., Am.J.Phys. 51, 712, 1983.
*
* P.T.Wallace Starlink 31 December 2002
*
* Copyright (C) 2002 Rutherford Appleton Laboratory
*
* License:
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program (see SLA_CONDITIONS); if not, write to the
* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
* Boston, MA 02110-1301 USA
*
* Copyright (C) 1995 Association of Universities for Research in Astronomy Inc.
*-
IMPLICIT NONE
DOUBLE PRECISION DATE
INTEGER JFORM
DOUBLE PRECISION EPOCH,ORBINC,ANODE,PERIH,AORQ,E,AORL,DM,PV(6)
INTEGER JSTAT
DOUBLE PRECISION U(13)
INTEGER J
* Validate elements and convert to "universal variables" parameters.
CALL slELUE(DATE,JFORM,
: EPOCH,ORBINC,ANODE,PERIH,AORQ,E,AORL,DM,U,J)
* Determine the position and velocity.
IF (J.EQ.0) THEN
CALL slUEPV(DATE,U,PV,J)
IF (J.NE.0) J=-5
END IF
* Wrap up.
JSTAT = J
END
|
