| 1 | #include "slalib.h"
|
|---|
| 2 | #include "slamac.h"
|
|---|
| 3 | void 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 | }
|
|---|
