1 /*
2 *+
3 * Name:
4 * palUe2el
5
6 * Purpose:
7 * Universal elements to heliocentric osculating elements
8
9 * Language:
10 * Starlink ANSI C
11
12 * Type of Module:
13 * Library routine
14
15 * Invocation:
16 * void palUe2el ( const double u[13], int jformr,
17 * int *jform, double *epoch, double *orbinc,
18 * double *anode, double *perih, double *aorq, double *e,
19 * double *aorl, double *dm, int *jstat );
20
21 * Arguments:
22 * u = const double [13] (Given)
23 * Universal orbital elements (Note 1)
24 * (0) combined mass (M+m)
25 * (1) total energy of the orbit (alpha)
26 * (2) reference (osculating) epoch (t0)
27 * (3-5) position at reference epoch (r0)
28 * (6-8) velocity at reference epoch (v0)
29 * (9) heliocentric distance at reference epoch
30 * (10) r0.v0
31 * (11) date (t)
32 * (12) universal eccentric anomaly (psi) of date, approx
33 * jformr = int (Given)
34 * Requested element set (1-3; Note 3)
35 * jform = int * (Returned)
36 * Element set actually returned (1-3; Note 4)
37 * epoch = double * (Returned)
38 * Epoch of elements (TT MJD)
39 * orbinc = double * (Returned)
40 * inclination (radians)
41 * anode = double * (Returned)
42 * longitude of the ascending node (radians)
43 * perih = double * (Returned)
44 * longitude or argument of perihelion (radians)
45 * aorq = double * (Returned)
46 * mean distance or perihelion distance (AU)
47 * e = double * (Returned)
48 * eccentricity
49 * aorl = double * (Returned)
50 * mean anomaly or longitude (radians, JFORM=1,2 only)
51 * dm = double * (Returned)
52 * daily motion (radians, JFORM=1 only)
53 * jstat = int * (Returned)
54 * status: 0 = OK
55 * -1 = illegal combined mass
56 * -2 = illegal JFORMR
57 * -3 = position/velocity out of range
58
59 * Description:
60 * Transform universal elements into conventional heliocentric
61 * osculating elements.
62
63 * Authors:
64 * PTW: Patrick T. Wallace
65 * TIMJ: Tim Jenness (JAC, Hawaii)
66 * {enter_new_authors_here}
67
68 * Notes:
69 * - The "universal" elements are those which define the orbit for the
70 * purposes of the method of universal variables (see reference 2).
71 * They consist of the combined mass of the two bodies, an epoch,
72 * and the position and velocity vectors (arbitrary reference frame)
73 * at that epoch. The parameter set used here includes also various
74 * quantities that can, in fact, be derived from the other
75 * information. This approach is taken to avoiding unnecessary
76 * computation and loss of accuracy. The supplementary quantities
77 * are (i) alpha, which is proportional to the total energy of the
78 * orbit, (ii) the heliocentric distance at epoch, (iii) the
79 * outwards component of the velocity at the given epoch, (iv) an
80 * estimate of psi, the "universal eccentric anomaly" at a given
81 * date and (v) that date.
82 * - The universal elements are with respect to the mean equator and
83 * equinox of epoch J2000. The orbital elements produced are with
84 * respect to the J2000 ecliptic and mean equinox.
85 * - Three different element-format options are supported:
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
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
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
116 *
117 * - It may not be possible to generate elements in the form
118 * requested through JFORMR. The caller is notified of the form
119 * of elements actually returned by means of the JFORM argument:
120 *
121 * JFORMR JFORM meaning
122 *
123 * 1 1 OK - elements are in the requested format
124 * 1 2 never happens
125 * 1 3 orbit not elliptical
126 *
127 * 2 1 never happens
128 * 2 2 OK - elements are in the requested format
129 * 2 3 orbit not elliptical
130 *
131 * 3 1 never happens
132 * 3 2 never happens
133 * 3 3 OK - elements are in the requested format
134 *
135 * - The arguments returned for each value of JFORM (cf Note 6: JFORM
136 * may not be the same as JFORMR) are as follows:
137 *
138 * JFORM 1 2 3
139 * EPOCH t0 t0 T
140 * ORBINC i i i
141 * ANODE Omega Omega Omega
142 * PERIH curly pi omega omega
143 * AORQ a a q
144 * E e e e
145 * AORL L M -
146 * DM n - -
147 *
148 * where:
149 *
150 * t0 is the epoch of the elements (MJD, TT)
151 * T " epoch of perihelion (MJD, TT)
152 * i " inclination (radians)
153 * Omega " longitude of the ascending node (radians)
154 * curly pi " longitude of perihelion (radians)
155 * omega " argument of perihelion (radians)
156 * a " mean distance (AU)
157 * q " perihelion distance (AU)
158 * e " eccentricity
159 * L " longitude (radians, 0-2pi)
160 * M " mean anomaly (radians, 0-2pi)
161 * n " daily motion (radians)
162 * - means no value is set
163 *
164 * - At very small inclinations, the longitude of the ascending node
165 * ANODE becomes indeterminate and under some circumstances may be
166 * set arbitrarily to zero. Similarly, if the orbit is close to
167 * circular, the true anomaly becomes indeterminate and under some
168 * circumstances may be set arbitrarily to zero. In such cases,
169 * the other elements are automatically adjusted to compensate,
170 * and so the elements remain a valid description of the orbit.
171
172 * See Also:
173 * - Sterne, Theodore E., "An Introduction to Celestial Mechanics",
174 * Interscience Publishers Inc., 1960. Section 6.7, p199.
175 * - Everhart, E. & Pitkin, E.T., Am.J.Phys. 51, 712, 1983.
176
177 * History:
178 * 2012-03-09 (TIMJ):
179 * Initial version
180 * Adapted with permission from the Fortran SLALIB library.
181 * {enter_further_changes_here}
182
183 * Copyright:
184 * Copyright (C) 1999 Rutherford Appleton Laboratory
185 * Copyright (C) 2012 Science and Technology Facilities Council.
186 * All Rights Reserved.
187
188 * Licence:
189 * This program is free software; you can redistribute it and/or
190 * modify it under the terms of the GNU General Public License as
191 * published by the Free Software Foundation; either version 3 of
192 * the License, or (at your option) any later version.
193 *
194 * This program is distributed in the hope that it will be
195 * useful, but WITHOUT ANY WARRANTY; without even the implied
196 * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
197 * PURPOSE. See the GNU General Public License for more details.
198 *
199 * You should have received a copy of the GNU General Public License
200 * along with this program; if not, write to the Free Software
201 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
202 * MA 02110-1301, USA.
203
204 * Bugs:
205 * {note_any_bugs_here}
206 *-
207 */
208
209 #include "pal.h"
210 #include "palmac.h"
211
palUe2el(const double u[],int jformr,int * jform,double * epoch,double * orbinc,double * anode,double * perih,double * aorq,double * e,double * aorl,double * dm,int * jstat)212 void palUe2el ( const double u[], int jformr,
213 int *jform, double *epoch, double *orbinc,
214 double *anode, double *perih, double *aorq, double *e,
215 double *aorl, double *dm, int *jstat ) {
216
217 /* Canonical days to seconds */
218 const double CD2S = PAL__GCON / PAL__SPD;
219
220 int i;
221 double pmass, date, pv[6];
222
223 /* Unpack the universal elements */
224 pmass = u[0] - 1.0;
225 date = u[2];
226 for (i=0; i<3; i++) {
227 pv[i] = u[i+3];
228 pv[i+3] = u[i+6] * CD2S;
229 }
230
231 /* Convert the position and velocity etc into conventional elements */
232 palPv2el( pv, date, pmass, jformr, jform, epoch, orbinc, anode,
233 perih, aorq, e, aorl, dm, jstat );
234 }
235