source: trunk/MagicSoft/slalib/pv2el.c

Last change on this file was 731, checked in by tbretz, 24 years ago
*** empty log message ***
  • Property svn:executable set to *
File size: 10.4 KB
Line 
1#include "slalib.h"
2#include "slamac.h"
3void slaPv2el ( double pv[], double date, double pmass, int jformr,
4 int *jform, double *epoch, double *orbinc,
5 double *anode, double *perih, double *aorq, double *e,
6 double *aorl, double *dm, int *jstat )
7/*
8** - - - - - - - - -
9** s l a P v 2 e l
10** - - - - - - - - -
11**
12** Heliocentric osculating elements obtained from instantaneous position
13** and velocity.
14**
15** Given:
16** pv double[6] heliocentric x,y,z,xdot,ydot,zdot of date,
17** J2000 equatorial triad (AU,AU/s; Note 1)
18** date double date (TT Modified Julian Date = JD-2400000.5)
19** pmass double mass of the planet (Sun=1; Note 2)
20** jformr int requested element set (1-3; Note 3)
21**
22** Returned:
23** jform double* element set actually returned (1-3; Note 4)
24** epoch double* epoch of elements (TT MJD)
25** orbinc double* inclination (radians)
26** anode double* longitude of the ascending node (radians)
27** perih double* longitude or argument of perihelion (radians)
28** aorq double* mean distance or perihelion distance (AU)
29** e double* eccentricity
30** aorl double* mean anomaly or longitude (radians, jform=1,2 only)
31** dm double* daily motion (radians, jform=1 only)
32** jstat int* status: 0 = OK
33** -1 = illegal pmass
34** -2 = illegal jformr
35** -3 = position/velocity out of range
36**
37** Notes
38**
39** 1 The pv 6-vector is with respect to the mean equator and equinox of
40** epoch J2000. The orbital elements produced are with respect to
41** the J2000 ecliptic and mean equinox.
42**
43** 2 The mass, pmass, is important only for the larger planets. For
44** most purposes (e.g. asteroids) use 0.0. Values less than zero
45** are illegal.
46**
47** 3 Three different element-format options are supported:
48**
49** Option jformr=1, suitable for the major planets:
50**
51** epoch = epoch of elements (TT MJD)
52** orbinc = inclination i (radians)
53** anode = longitude of the ascending node, big omega (radians)
54** perih = longitude of perihelion, curly pi (radians)
55** aorq = mean distance, a (AU)
56** e = eccentricity, e
57** aorl = mean longitude L (radians)
58** dm = daily motion (radians)
59**
60** Option jformr=2, suitable for minor planets:
61**
62** epoch = epoch of elements (TT MJD)
63** orbinc = inclination i (radians)
64** anode = longitude of the ascending node, big omega (radians)
65** perih = argument of perihelion, little omega (radians)
66** aorq = mean distance, a (AU)
67** e = eccentricity, e
68** aorl = mean anomaly M (radians)
69**
70** Option jformr=3, suitable for comets:
71**
72** epoch = epoch of perihelion (TT MJD)
73** orbinc = inclination i (radians)
74** anode = longitude of the ascending node, big omega (radians)
75** perih = argument of perihelion, little omega (radians)
76** aorq = perihelion distance, q (AU)
77** e = eccentricity, e
78**
79** 4 It may not be possible to generate elements in the form
80** requested through jformr. The caller is notified of the form
81** of elements actually returned by means of the jform argument:
82**
83** jformr jform meaning
84**
85** 1 1 OK - elements are in the requested format
86** 1 2 never happens
87** 1 3 orbit not elliptical
88**
89** 2 1 never happens
90** 2 2 OK - elements are in the requested format
91** 2 3 orbit not elliptical
92**
93** 3 1 never happens
94** 3 2 never happens
95** 3 3 OK - elements are in the requested format
96**
97** 5 The arguments returned for each value of jform (cf Note 5: jform
98** may not be the same as jformr) are as follows:
99**
100** jform 1 2 3
101** epoch t0 t0 T
102** orbinc i i i
103** anode Omega Omega Omega
104** perih curly pi omega omega
105** aorq a a q
106** e e e e
107** aorl L M -
108** dm n - -
109**
110** where:
111**
112** t0 is the epoch of the elements (MJD, TT)
113** T " epoch of perihelion (MJD, TT)
114** i " inclination (radians)
115** Omega " longitude of the ascending node (radians)
116** curly pi " longitude of perihelion (radians)
117** omega " argument of perihelion (radians)
118** a " mean distance (AU)
119** q " perihelion distance (AU)
120** e " eccentricity
121** L " longitude (radians, 0-2pi)
122** M " mean anomaly (radians, 0-2pi)
123** n " daily motion (radians)
124** - means no value is set
125**
126** 6 At very small inclinations, the longitude of the ascending node
127** anode becomes indeterminate and under some circumstances may be
128** set arbitrarily to zero. Similarly, if the orbit is close to
129** circular, the true anomaly becomes indeterminate and under some
130** circumstances may be set arbitrarily to zero. In such cases,
131** the other elements are automatically adjusted to compensate,
132** and so the elements remain a valid description of the orbit.
133**
134** Reference: Sterne, Theodore E., "An Introduction to Celestial
135** Mechanics", Interscience Publishers, 1960
136**
137** Called: slaDranrm
138**
139** Last revision: 21 February 1999
140**
141** Copyright P.T.Wallace. All rights reserved.
142*/
143
144/* Seconds to days */
145#define DAY 86400.0
146
147/* Gaussian gravitational constant (exact) */
148#define GCON 0.01720209895
149
150/* Sin and cos of J2000 mean obliquity (IAU 1976) */
151#define SE 0.3977771559319137
152#define CE 0.9174820620691818
153
154/* Minimum allowed distance (AU) and speed (AU/day) */
155#define RMIN 1e-3
156#define VMIN 1e-8
157
158/* How close to unity the eccentricity has to be to call it a parabola */
159#define PARAB 1e-8
160
161{
162 double x, y, z, xd, yd, zd, r, v2, v, rdv, gmu, hx, hy, hz,
163 hx2py2, h2, h, oi, bigom, ar, e2, ecc, s, c, at, u, om,
164 gar3, em1, ep1, hat, shat, chat, ae, am, dn, pl,
165 el, q, tp, that, thhf, f;
166 int jf;
167
168
169/* Validate arguments pmass and jformr. */
170 if ( pmass < 0.0 ) {
171 *jstat = -1;
172 return;
173 }
174 if ( jformr < 1 || jformr > 3 ) {
175 *jstat = -2;
176 return;
177 }
178
179/* Provisionally assume the elements will be in the chosen form. */
180 jf = jformr;
181
182/* Rotate the position from equatorial to ecliptic coordinates. */
183 x = pv [ 0 ];
184 y = pv [ 1 ] * CE + pv [ 2 ] * SE;
185 z = - pv [ 1 ] * SE + pv [ 2 ] * CE;
186
187/* Rotate the velocity similarly, scaling to AU/day. */
188 xd = DAY * pv [ 3 ];
189 yd = DAY * ( pv [ 4 ] * CE + pv [ 5 ] * SE );
190 zd = DAY * ( - pv [ 4 ] * SE + pv [ 5 ] * CE );
191
192/* Distance and speed. */
193 r = sqrt ( x * x + y * y + z * z );
194 v2 = xd * xd + yd * yd + zd * zd;
195 v = sqrt ( v2 );
196
197/* Reject unreasonably small values. */
198 if ( r < RMIN || v < VMIN ) {
199 *jstat = -3;
200 return;
201 }
202
203/* R dot V. */
204 rdv = x * xd + y * yd + z * zd;
205
206/* Mu. */
207 gmu = ( 1.0 + pmass ) * GCON * GCON;
208
209/* Vector angular momentum per unit reduced mass. */
210 hx = y * zd - z * yd;
211 hy = z * xd - x * zd;
212 hz = x * yd - y * xd;
213
214/* Areal constant. */
215 hx2py2 = hx * hx + hy * hy;
216 h2 = hx2py2 + hz * hz;
217 h = sqrt ( h2 );
218
219/* Inclination. */
220 oi = atan2 ( sqrt ( hx2py2 ), hz );
221
222/* Longitude of ascending node. */
223 if ( hx != 0.0 || hy != 0.0 ) {
224 bigom = atan2 ( hx, -hy );
225 } else {
226 bigom = 0.0;
227 }
228
229/* Reciprocal of mean distance etc. */
230 ar = 2.0 / r - v2 / gmu;
231
232/* Eccentricity. */
233 e2 = 1.0 - ar * h2 / gmu;
234 ecc = ( e2 >= 0.0 ) ? sqrt ( e2 ) : 0.0;
235
236/* True anomaly. */
237 s = h * rdv;
238 c = h2 - r * gmu;
239 if ( s != 0.0 && c != 0.0 ) {
240 at = atan2 ( s, c );
241 } else {
242 at = 0.0;
243 }
244
245/* Argument of the latitude. */
246 s = sin ( bigom );
247 c = cos ( bigom );
248 u = atan2 ( ( - x * s + y * c ) * cos ( oi ) + z * sin ( oi ),
249 x * c + y * s );
250
251/* Argument of perihelion. */
252 om = u - at;
253
254/* Capture near-parabolic cases. */
255 if ( fabs ( ecc - 1.0 ) < PARAB ) ecc = 1.0;
256
257/* Comply with jformr = 1 or 2 only if orbit is elliptical. */
258 if ( ecc >= 1.0 ) jf = 3;
259
260/* Functions. */
261 gar3 = gmu * ar * ar * ar;
262 em1 = ecc - 1.0;
263 ep1 = ecc + 1.0;
264 hat = at / 2.0;
265 shat = sin ( hat );
266 chat = cos ( hat );
267
268/* Ellipse? */
269 if ( ecc < 1.0 ) {
270
271 /* Eccentric anomaly. */
272 ae = 2.0 * atan2 ( sqrt ( -em1 ) * shat, sqrt ( ep1 ) * chat );
273
274 /* Mean anomaly. */
275 am = ae - ecc * sin ( ae );
276
277 /* Daily motion. */
278 dn = sqrt ( gar3 );
279 }
280
281/* "Major planet" element set? */
282 if ( jf == 1 ) {
283
284 /* Longitude of perihelion. */
285 pl = bigom + om;
286
287 /* Longitude at epoch. */
288 el = pl + am;
289 }
290
291/* "Comet" element set? */
292 if ( jf == 3 ) {
293
294 /* Perihelion distance. */
295 q = h2 / ( gmu * ep1 );
296
297 /* Ellipse, parabola, hyperbola? */
298 if ( ecc < 1.0 ) {
299
300 /* Ellipse: epoch of perihelion. */
301 tp = date - am / dn;
302 } else {
303
304 /* Parabola or hyperbola: evaluate tan ( ( true anomaly ) / 2 ) */
305 that = shat / chat;
306 if ( ecc == 1.0 ) {
307
308 /* Parabola: epoch of perihelion. */
309 tp = date - that * ( 1.0 + that * that / 3.0 ) * h * h2 /
310 ( 2.0 * gmu * gmu );
311 } else {
312
313 /* Hyperbola: epoch of perihelion. */
314 thhf = sqrt ( em1 / ep1 ) * that;
315 f = log ( 1.0 + thhf ) - log ( 1.0 - thhf );
316 tp = date - ( ecc * sinh ( f ) - f ) / sqrt ( - gar3 );
317 }
318 }
319 }
320
321/* Return the appropriate set of elements. */
322 *jform = jf;
323 *orbinc = oi;
324 *anode = slaDranrm ( bigom );
325 *e = ecc;
326 if ( jf == 1 ) {
327 *perih = slaDranrm ( pl );
328 *aorl = slaDranrm ( el );
329 *dm = dn;
330 } else {
331 *perih = slaDranrm ( om );
332 if ( jf == 2 ) *aorl = slaDranrm ( am );
333 }
334 if ( jf != 3 ) {
335 *epoch = date;
336 *aorq = 1.0 / ar;
337 } else {
338 *epoch = tp;
339 *aorq = q;
340 }
341 *jstat = 0;
342
343}
Note: See TracBrowser for help on using the repository browser.