source: trunk/MagicSoft/slalib/pv2ue.c@ 2514

Last change on this file since 2514 was 731, checked in by tbretz, 24 years ago
*** empty log message ***
  • Property svn:executable set to *
File size: 4.0 KB
Line 
1#include "slalib.h"
2#include "slamac.h"
3void slaPv2ue ( double pv[], double date, double pmass,
4 double u[], int *jstat )
5/*
6** - - - - - - - - -
7** s l a P v 2 u e
8** - - - - - - - - -
9**
10** Construct a universal element set based on an instantaneous position
11** and velocity.
12**
13** Given:
14** pv double[6] heliocentric x,y,z,xdot,ydot,zdot of date,
15** (au,au/s; Note 1)
16** date double date (TT Modified Julian Date = JD-2400000.5)
17** pmass double mass of the planet (Sun=1; Note 2)
18**
19** Returned:
20**
21** u double[13] universal orbital elements (Note 3)
22**
23** [0] combined mass (M+m)
24** [1] total energy of the orbit (alpha)
25** [2] reference (osculating) epoch (t0)
26** [3-5] position at reference epoch (r0)
27** [6-8] velocity at reference epoch (v0)
28** [9] heliocentric distance at reference epoch
29** [10] r0.v0
30** [11] date (t)
31** [12] universal eccentric anomaly (psi) of date
32**
33** jstat int* status: 0 = OK
34** -1 = illegal pmass
35** -2 = too close to Sun
36** -3 = too slow
37**
38** Notes
39**
40** 1 The pv 6-vector can be with respect to any chosen inertial frame,
41** and the resulting universal-element set will be with respect to
42** the same frame. A common choice will be mean equator and ecliptic
43** of epoch J2000.
44**
45** 2 The mass, pmass, is important only for the larger planets. For
46** most purposes (e.g. asteroids) use 0.0. Values less than zero
47** are illegal.
48**
49** 3 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** Reference: Everhart, E. & Pitkin, E.T., Am.J.Phys. 51, 712, 1983.
64**
65** Last revision: 17 March 1999
66**
67** Copyright P.T.Wallace. All rights reserved.
68*/
69
70/* Gaussian gravitational constant (exact) */
71#define GCON 0.01720209895
72
73/* Canonical days to seconds */
74#define CD2S ( GCON / 86400.0 );
75
76/* Minimum allowed distance (AU) and speed (AU per canonical day) */
77#define RMIN 1e-3
78#define VMIN 1e-3
79
80{
81 double t0, cm, x, y, z, xd, yd, zd, r, v2, v, alpha, rdv;
82
83
84/* Reference epoch. */
85 t0 = date;
86
87/* Combined mass (mu=M+m). */
88 if ( pmass < 0.0 ) {
89 *jstat = -1;
90 return;
91 }
92 cm = 1.0 + pmass;
93
94/* Unpack the state vector, expressing velocity in AU per canonical day. */
95 x = pv[0];
96 y = pv[1];
97 z = pv[2];
98 xd = pv[3] / CD2S;
99 yd = pv[4] / CD2S;
100 zd = pv[5] / CD2S;
101
102/* Heliocentric distance, and speed. */
103 r = sqrt ( x * x + y * y + z * z );
104 v2 = xd * xd + yd * yd + zd * zd;
105 v = sqrt ( v2 );
106
107/* Reject unreasonably small values. */
108 if ( r < RMIN ) {
109 *jstat = -2;
110 return;
111 }
112 if ( v < VMIN ) {
113 *jstat = -3;
114 return;
115 }
116
117/* Total energy of the orbit. */
118 alpha = v2 - 2.0 * cm / r;
119
120/* Outward component of velocity. */
121 rdv = x * xd + y * yd + z * zd;
122
123/* Construct the universal-element set. */
124 u[0] = cm;
125 u[1] = alpha;
126 u[2] = t0;
127 u[3] = x;
128 u[4] = y;
129 u[5] = z;
130 u[6] = xd;
131 u[7] = yd;
132 u[8] = zd;
133 u[9 ] = r;
134 u[10] = rdv;
135 u[11] = t0;
136 u[12] = 0.0;
137
138/* Exit. */
139 *jstat = 0;
140
141}
Note: See TracBrowser for help on using the repository browser.