source: trunk/MagicSoft/slalib/el2ue.c@ 9022

Last change on this file since 9022 was 731, checked in by tbretz, 25 years ago
*** empty log message ***
  • Property svn:executable set to *
File size: 9.5 KB
Line 
1#include "slalib.h"
2#include "slamac.h"
3void slaEl2ue ( double date, int jform, double epoch, double orbinc,
4 double anode, double perih, double aorq, double e,
5 double aorl, double dm, double u[], int *jstat )
6/*
7** - - - - - - - - -
8** s l a E l 2 u e
9** - - - - - - - - -
10**
11** Transform conventional osculating orbital elements into "universal" form.
12**
13** Given:
14** date double epoch (TT MJD) of osculation (Note 3)
15** jform int choice of element set (1-3, Note 6)
16** epoch double epoch (TT MJD) of the elements
17** orbinc double inclination (radians)
18** anode double longitude of the ascending node (radians)
19** perih double longitude or argument of perihelion (radians)
20** aorq double mean distance or perihelion distance (AU)
21** e double eccentricity
22** aorl double mean anomaly or longitude (radians, jform=1,2 only)
23** dm double daily motion (radians, jform=1 only)
24**
25** Returned:
26** u double[13] universal orbital elements (Note 1)
27**
28** [0] combined mass (M+m)
29** [1] total energy of the orbit (alpha)
30** [2] reference (osculating) epoch (t0)
31** [3-5] position at reference epoch (r0)
32** [6-8] velocity at reference epoch (v0)
33** [9] heliocentric distance at reference epoch
34** [10] r0.v0
35** [11] date (t)
36** [12] universal eccentric anomaly (psi) of date, approx
37**
38** jstat int* status: 0 = OK
39** -1 = illegal jform
40** -2 = illegal e
41** -3 = illegal aorq
42** -4 = illegal dm
43** -5 = numerical error
44**
45** Called: slaUe2pv, slaPv2ue
46**
47** Notes
48**
49** 1 The "universal" elements are those which define the orbit for the
50** purposes of the method of universal variables (see reference).
51** They consist of the combined mass of the two bodies, an epoch,
52** and the position and velocity vectors (arbitrary reference frame)
53** at that epoch. The parameter set used here includes also various
54** quantities that can, in fact, be derived from the other
55** information. This approach is taken to avoiding unnecessary
56** computation and loss of accuracy. The supplementary quantities
57** are (i) alpha, which is proportional to the total energy of the
58** orbit, (ii) the heliocentric distance at epoch, (iii) the
59** outwards component of the velocity at the given epoch, (iv) an
60** estimate of psi, the "universal eccentric anomaly" at a given
61** date and (v) that date.
62**
63** 2 The companion routine is slaUe2pv. This takes the set of numbers
64** that the present routine outputs and uses them to derive the
65** object's position and velocity. A single prediction requires one
66** call to the present routine followed by one call to slaUe2pv;
67** for convenience, the two calls are packaged as the routine
68** slaPlanel. Multiple predictions may be made by again calling the
69** present routine once, but then calling slaUe2pv multiple times,
70** which is faster than multiple calls to slaPlanel.
71**
72** 3 date is the epoch of osculation. It is in the TT timescale
73** (formerly Ephemeris Time, ET) and is a Modified Julian Date
74** (JD-2400000.5).
75**
76** 4 The supplied orbital elements are with respect to the J2000
77** ecliptic and equinox. The position and velocity parameters
78** returned in the array u are with respect to the mean equator and
79** equinox of epoch J2000, and are for the perihelion prior to the
80** specified epoch.
81**
82** 5 The universal elements returned in the array u are in canonical
83** units (solar masses, AU and canonical days).
84**
85** 6 Three different element-format options are available:
86**
87** Option jform=1, suitable for the major planets:
88**
89** epoch = epoch of elements (TT MJD)
90** orbinc = inclination i (radians)
91** anode = longitude of the ascending node, big omega (radians)
92** perih = longitude of perihelion, curly pi (radians)
93** aorq = mean distance, a (AU)
94** e = eccentricity, e (range 0 to <1)
95** aorl = mean longitude L (radians)
96** dm = daily motion (radians)
97**
98** Option jform=2, suitable for minor planets:
99**
100** epoch = epoch of elements (TT MJD)
101** orbinc = inclination i (radians)
102** anode = longitude of the ascending node, big omega (radians)
103** perih = argument of perihelion, little omega (radians)
104** aorq = mean distance, a (AU)
105** e = eccentricity, e (range 0 to <1)
106** aorl = mean anomaly M (radians)
107**
108** Option jform=3, suitable for comets:
109**
110** epoch = epoch of perihelion (TT MJD)
111** orbinc = inclination i (radians)
112** anode = longitude of the ascending node, big omega (radians)
113** perih = argument of perihelion, little omega (radians)
114** aorq = perihelion distance, q (AU)
115** e = eccentricity, e (range 0 to 10)
116**
117** 7 Unused elements (dm for jform=2, aorl and dm for jform=3) are
118** not accessed.
119**
120** 8 The algorithm was originally adapted from the EPHSLA program of
121** D.H.P.Jones (private communication, 1996). The method is based on
122** Stumpff's Universal Variables.
123**
124** Reference: Everhart, E. & Pitkin, E.T., Am.J.Phys. 51, 712, 1983.
125**
126** Last revision: 18 March 1999
127**
128** Copyright P.T.Wallace. All rights reserved.
129*/
130
131/* Gaussian gravitational constant (exact) */
132#define GCON 0.01720209895
133
134/* Sin and cos of J2000 mean obliquity (IAU 1976) */
135#define SE 0.3977771559319137
136#define CE 0.9174820620691818
137
138{
139 int j;
140 double pht, argph, q, w, cm, alpha, phs, sw, cw, si, ci, so, co,
141 x, y, z, px, py, pz, vx, vy, vz, dt, fc, fp, psi, ul[13], pv[6];
142
143
144
145/* Validate arguments. */
146 if ( jform < 1 || jform > 3 ) {
147 *jstat = -1;
148 return;
149 }
150 if ( e < 0.0 || e > 10.0 || ( e >= 1.0 && jform != 3 ) ) {
151 *jstat = -2;
152 return;
153 }
154 if ( aorq <= 0.0 ) {
155 *jstat = -3;
156 return;
157 }
158 if ( jform == 1 && dm <= 0.0 ) {
159 *jstat = -4;
160 return;
161 }
162
163/*
164** Transform elements into standard form:
165**
166** pht = epoch of perihelion passage
167** argph = argument of perihelion (little omega)
168** q = perihelion distance (q)
169** cm = combined mass, M+m (mu)
170*/
171
172 if ( jform == 1 ) {
173 pht = epoch - ( aorl - perih ) / dm;
174 argph = perih - anode;
175 q = aorq * ( 1.0 - e );
176 w = dm / GCON;
177 cm = w * w * aorq * aorq * aorq;
178 } else if ( jform == 2 ) {
179 pht = epoch - aorl * sqrt ( aorq * aorq * aorq ) / GCON;
180 argph = perih;
181 q = aorq * ( 1.0 - e );
182 cm = 1.0;
183 } else if ( jform == 3 ) {
184 pht = epoch;
185 argph = perih;
186 q = aorq;
187 cm = 1.0;
188 }
189
190/*
191** The universal variable alpha. This is proportional to the total
192** energy of the orbit: -ve for an ellipse, zero for a parabola,
193** +ve for a hyperbola.
194*/
195
196 alpha = cm * ( e - 1.0 ) / q;
197
198/* Speed at perihelion. */
199
200 phs = sqrt ( alpha + 2.0 * cm / q );
201
202/*
203** In a Cartesian coordinate system which has the x-axis pointing
204** to perihelion and the z-axis normal to the orbit (such that the
205** object orbits counter-clockwise as seen from +ve z), the
206** perihelion position and velocity vectors are:
207**
208** position [Q,0,0]
209** velocity [0,phs,0]
210**
211** To express the results in J2000 equatorial coordinates we make a
212** series of four rotations of the Cartesian axes:
213**
214** axis Euler angle
215**
216** 1 z argument of perihelion (little omega)
217** 2 x inclination (i)
218** 3 z longitude of the ascending node (big omega)
219** 4 x J2000 obliquity (epsilon)
220**
221** In each case the rotation is clockwise as seen from the +ve end
222** of the axis concerned.
223*/
224
225/* Functions of the Euler angles. */
226 sw = sin ( argph );
227 cw = cos ( argph );
228 si = sin ( orbinc );
229 ci = cos ( orbinc );
230 so = sin ( anode );
231 co = cos ( anode );
232
233/* Position at perihelion (AU). */
234 x = q * cw;
235 y = q * sw;
236 z = y * si;
237 y = y * ci;
238 px = x * co - y * so;
239 y = x * so + y * co;
240 py = y * CE - z * SE;
241 pz = y * SE + z * CE;
242
243/* Velocity at perihelion (AU per canonical day). */
244 x = - phs * sw;
245 y = phs * cw;
246 z = y * si;
247 y = y * ci;
248 vx = x * co - y * so;
249 y = x * so + y * co;
250 vy = y * CE - z * SE;
251 vz = y * SE + z * CE;
252
253/* Time from perihelion to date (in Canonical Days: a canonical */
254/* day is 58.1324409... days, defined as 1/GCON). */
255
256 dt = ( date - pht ) * GCON;
257
258/* First Approximation to the Universal Eccentric Anomaly, psi, */
259/* based on the circle (fc) and parabola (fp) values. */
260 fc = dt / q;
261 w = pow ( 3.0 * dt + sqrt ( 9.0 * dt * dt + 8.0 * q * q * q ),
262 1.0 / 3.0 );
263 fp = w - 2.0 * q / w;
264 psi = ( 1.0 - e ) * fc + e * fp;
265
266/* Assemble local copy of element set. */
267 ul[0] = cm;
268 ul[1] = alpha;
269 ul[2] = pht;
270 ul[3] = px;
271 ul[4] = py;
272 ul[5] = pz;
273 ul[6] = vx;
274 ul[7] = vy;
275 ul[8] = vz;
276 ul[9] = q;
277 ul[10] = 0.0;
278 ul[11] = date;
279 ul[12] = psi;
280
281/* Predict position+velocity at epoch of osculation. */
282 slaUe2pv ( date, ul, pv, &j );
283 if ( j ) {
284 *jstat = -5;
285 return;
286 }
287
288/* Convert back to universal elements. */
289 slaPv2ue ( pv, date, cm - 1.0, u, &j );
290 if ( j ) {
291 *jstat = -5;
292 return;
293 }
294
295/* OK exit. */
296 *jstat = 0;
297
298}
Note: See TracBrowser for help on using the repository browser.