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
|
SUBROUTINE sla_PLANEL (DATE, JFORM, EPOCH, ORBINC, ANODE, PERIH,
: AORQ, E, AORL, DM, PV, JSTAT)
*+
* - - - - - - -
* P L A N E 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)
* JFORM i choice of element set (1-3; Note 3)
* EPOCH d epoch of elements (TT MJD)
* 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: sla_EL2UE, sla_UE2PV
*
* 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 Three different element-format options are 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 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)
*
* 4 Unused elements (DM for JFORM=2, AORL and DM for JFORM=3) are
* not accessed.
*
* 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 18 March 1999
*
* Copyright (C) 1999 Rutherford Appleton Laboratory
*-
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 sla_EL2UE(DATE,JFORM,
: EPOCH,ORBINC,ANODE,PERIH,AORQ,E,AORL,DM,U,J)
* Determine the position and velocity.
IF (J.EQ.0) THEN
CALL sla_UE2PV(DATE,U,PV,J)
IF (J.NE.0) J=-5
END IF
* Wrap up.
JSTAT = J
END
|
