1 /*
2 *+
3 * Name:
4 * palPv2ue
5
6 * Purpose:
7 * Universal elements to position and velocity.
8
9 * Language:
10 * Starlink ANSI C
11
12 * Type of Module:
13 * Library routine
14
15 * Invocation:
16 * void palPv2ue( const double pv[6], double date, double pmass,
17 * double u[13], int * jstat );
18
19 * Arguments:
20 * pv = double [6] (Given)
21 * Heliocentric x,y,z,xdot,ydot,zdot of date, (AU,AU/s; Note 1)
22 * date = double (Given)
23 * Date (TT modified Julian Date = JD-2400000.5)
24 * pmass = double (Given)
25 * Mass of the planet (Sun=1; note 2)
26 * u = double [13] (Returned)
27 * Universal orbital elements (Note 3)
28 *
29 * - (0) combined mass (M+m)
30 * - (1) total energy of the orbit (alpha)
31 * - (2) reference (osculating) epoch (t0)
32 * - (3-5) position at reference epoch (r0)
33 * - (6-8) velocity at reference epoch (v0)
34 * - (9) heliocentric distance at reference epoch
35 * - (10) r0.v0
36 * - (11) date (t)
37 * - (12) universal eccentric anomaly (psi) of date, approx
38 * jstat = int * (Returned)
39 * status: 0 = OK
40 * - -1 = illegal PMASS
41 * - -2 = too close to Sun
42 * - -3 = too slow
43
44 * Description:
45 * Construct a universal element set based on an instantaneous position
46 * and velocity.
47
48
49 * Authors:
50 * PTW: Pat Wallace (STFC)
51 * TIMJ: Tim Jenness (JAC, Hawaii)
52 * {enter_new_authors_here}
53
54 * Notes:
55 * - The PV 6-vector can be with respect to any chosen inertial frame,
56 * and the resulting universal-element set will be with respect to
57 * the same frame. A common choice will be mean equator and ecliptic
58 * of epoch J2000.
59 * - The mass, PMASS, is important only for the larger planets. For
60 * most purposes (e.g. asteroids) use 0D0. Values less than zero
61 * are illegal.
62 * - The "universal" elements are those which define the orbit for the
63 * purposes of the method of universal variables (see reference).
64 * They consist of the combined mass of the two bodies, an epoch,
65 * and the position and velocity vectors (arbitrary reference frame)
66 * at that epoch. The parameter set used here includes also various
67 * quantities that can, in fact, be derived from the other
68 * information. This approach is taken to avoiding unnecessary
69 * computation and loss of accuracy. The supplementary quantities
70 * are (i) alpha, which is proportional to the total energy of the
71 * orbit, (ii) the heliocentric distance at epoch, (iii) the
72 * outwards component of the velocity at the given epoch, (iv) an
73 * estimate of psi, the "universal eccentric anomaly" at a given
74 * date and (v) that date.
75 * - Reference: Everhart, E. & Pitkin, E.T., Am.J.Phys. 51, 712, 1983.
76
77 * History:
78 * 2012-03-09 (TIMJ):
79 * Initial version from the SLA/F implementation.
80 * Adapted with permission from the Fortran SLALIB library.
81 * {enter_further_changes_here}
82
83 * Copyright:
84 * Copyright (C) 1999 Rutherford Appleton Laboratory
85 * Copyright (C) 2012 Science and Technology Facilities Council.
86 * All Rights Reserved.
87
88 * Licence:
89 * This program is free software; you can redistribute it and/or
90 * modify it under the terms of the GNU General Public License as
91 * published by the Free Software Foundation; either version 3 of
92 * the License, or (at your option) any later version.
93 *
94 * This program is distributed in the hope that it will be
95 * useful, but WITHOUT ANY WARRANTY; without even the implied
96 * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
97 * PURPOSE. See the GNU General Public License for more details.
98 *
99 * You should have received a copy of the GNU General Public License
100 * along with this program; if not, write to the Free Software
101 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
102 * MA 02110-1301, USA.
103
104 * Bugs:
105 * {note_any_bugs_here}
106 *-
107 */
108
109 #include <math.h>
110
111 #include "pal.h"
112 #include "palmac.h"
113
palPv2ue(const double pv[6],double date,double pmass,double u[13],int * jstat)114 void palPv2ue( const double pv[6], double date, double pmass,
115 double u[13], int * jstat ) {
116
117 /* Canonical days to seconds */
118 const double CD2S = PAL__GCON / PAL__SPD;
119
120 /* Minimum allowed distance (AU) and speed (AU per canonical day) */
121 const double RMIN = 1e-3;
122 const double VMIN = 1e-3;
123
124 double T0,CM,X,Y,Z,XD,YD,ZD,R,V2,V,ALPHA,RDV;
125
126 /* Reference epoch. */
127 T0 = date;
128
129 /* Combined mass (mu=M+m). */
130 if (pmass < 0.0 ) { /* Negative planet mass */
131 *jstat = -1;
132 return;
133 }
134 CM = 1.0+pmass;
135
136 /* Unpack the state vector, expressing velocity in AU per canonical day. */
137 X = pv[0];
138 Y = pv[1];
139 Z = pv[2];
140 XD = pv[3]/CD2S;
141 YD = pv[4]/CD2S;
142 ZD = pv[5]/CD2S;
143
144 /* Heliocentric distance, and speed. */
145 R = sqrt(X*X+Y*Y+Z*Z);
146 V2 = XD*XD+YD*YD+ZD*ZD;
147 V = sqrt(V2);
148
149 /* Reject unreasonably small values. */
150 if (R < RMIN) { /* Too close */
151 *jstat = -2;
152 return;
153 }
154 if (V < VMIN) { /* Too slow */
155 *jstat = -3;
156 return;
157 }
158
159 /* Total energy of the orbit. */
160 ALPHA = V2-2.0*CM/R;
161
162 /* Outward component of velocity. */
163 RDV = X*XD+Y*YD+Z*ZD;
164
165 /* Construct the universal-element set. */
166 u[0] = CM;
167 u[1] = ALPHA;
168 u[2] = T0;
169 u[3] = X;
170 u[4] = Y;
171 u[5] = Z;
172 u[6] = XD;
173 u[7] = YD;
174 u[8] = ZD;
175 u[9] = R;
176 u[10] = RDV;
177 u[11] = T0;
178 u[12] = 0.0;
179
180 *jstat = 0;
181 return;
182 }
183