source: trunk/FACT++/pal/palEl2ue.c

Last change on this file was 18347, checked in by tbretz, 9 years ago
File size: 10.7 KB
Line 
1/*
2*+
3* Name:
4* palEl2ue
5
6* Purpose:
7* Transform conventional elements into "universal" form
8
9* Language:
10* Starlink ANSI C
11
12* Type of Module:
13* Library routine
14
15* Invocation:
16* void palEl2ue ( double date, int jform, double epoch, double orbinc,
17* double anode, double perih, double aorq, double e,
18* double aorl, double dm, double u[13], int *jstat );
19
20* Arguments:
21* date = double (Given)
22* Epoch (TT MJD) of osculation (Note 3)
23* jform = int (Given)
24* Element set actually returned (1-3; Note 6)
25* epoch = double (Given)
26* Epoch of elements (TT MJD)
27* orbinc = double (Given)
28* inclination (radians)
29* anode = double (Given)
30* longitude of the ascending node (radians)
31* perih = double (Given)
32* longitude or argument of perihelion (radians)
33* aorq = double (Given)
34* mean distance or perihelion distance (AU)
35* e = double (Given)
36* eccentricity
37* aorl = double (Given)
38* mean anomaly or longitude (radians, JFORM=1,2 only)
39* dm = double (Given)
40* daily motion (radians, JFORM=1 only)
41* u = double [13] (Returned)
42* Universal orbital elements (Note 1)
43* - (0) combined mass (M+m)
44* - (1) total energy of the orbit (alpha)
45* - (2) reference (osculating) epoch (t0)
46* - (3-5) position at reference epoch (r0)
47* - (6-8) velocity at reference epoch (v0)
48* - (9) heliocentric distance at reference epoch
49* - (10) r0.v0
50* - (11) date (t)
51* - (12) universal eccentric anomaly (psi) of date, approx
52* jstat = int * (Returned)
53* status: 0 = OK
54* - -1 = illegal JFORM
55* - -2 = illegal E
56* - -3 = illegal AORQ
57* - -4 = illegal DM
58* - -5 = numerical error
59
60* Description:
61* Transform conventional osculating elements into "universal" form.
62
63* Authors:
64* PTW: Pat Wallace (STFC)
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).
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 companion routine is palUe2pv. This takes the set of numbers
83* that the present routine outputs and uses them to derive the
84* object's position and velocity. A single prediction requires one
85* call to the present routine followed by one call to palUe2pv;
86* for convenience, the two calls are packaged as the routine
87* palPlanel. Multiple predictions may be made by again calling the
88* present routine once, but then calling palUe2pv multiple times,
89* which is faster than multiple calls to palPlanel.
90* - DATE is the epoch of osculation. It is in the TT timescale
91* (formerly Ephemeris Time, ET) and is a Modified Julian Date
92* (JD-2400000.5).
93* - The supplied orbital elements are with respect to the J2000
94* ecliptic and equinox. The position and velocity parameters
95* returned in the array U are with respect to the mean equator and
96* equinox of epoch J2000, and are for the perihelion prior to the
97* specified epoch.
98* - The universal elements returned in the array U are in canonical
99* units (solar masses, AU and canonical days).
100* - Three different element-format options are available:
101*
102* Option JFORM=1, suitable for the major planets:
103*
104* EPOCH = epoch of elements (TT MJD)
105* ORBINC = inclination i (radians)
106* ANODE = longitude of the ascending node, big omega (radians)
107* PERIH = longitude of perihelion, curly pi (radians)
108* AORQ = mean distance, a (AU)
109* E = eccentricity, e (range 0 to <1)
110* AORL = mean longitude L (radians)
111* DM = daily motion (radians)
112*
113* Option JFORM=2, suitable for minor planets:
114*
115* EPOCH = epoch of elements (TT MJD)
116* ORBINC = inclination i (radians)
117* ANODE = longitude of the ascending node, big omega (radians)
118* PERIH = argument of perihelion, little omega (radians)
119* AORQ = mean distance, a (AU)
120* E = eccentricity, e (range 0 to <1)
121* AORL = mean anomaly M (radians)
122*
123* Option JFORM=3, suitable for comets:
124*
125* EPOCH = epoch of perihelion (TT MJD)
126* ORBINC = inclination i (radians)
127* ANODE = longitude of the ascending node, big omega (radians)
128* PERIH = argument of perihelion, little omega (radians)
129* AORQ = perihelion distance, q (AU)
130* E = eccentricity, e (range 0 to 10)
131*
132* - Unused elements (DM for JFORM=2, AORL and DM for JFORM=3) are
133* not accessed.
134* - The algorithm was originally adapted from the EPHSLA program of
135* D.H.P.Jones (private communication, 1996). The method is based
136* on Stumpff's Universal Variables.
137*
138* See Also:
139* Everhart & Pitkin, Am.J.Phys. 51, 712 (1983).
140
141* History:
142* 2012-03-12 (TIMJ):
143* Initial version taken directly from SLA/F.
144* Adapted with permission from the Fortran SLALIB library.
145* {enter_further_changes_here}
146
147* Copyright:
148* Copyright (C) 2005 Patrick T. Wallace
149* Copyright (C) 2012 Science and Technology Facilities Council.
150* All Rights Reserved.
151
152* Licence:
153* This program is free software; you can redistribute it and/or
154* modify it under the terms of the GNU General Public License as
155* published by the Free Software Foundation; either version 3 of
156* the License, or (at your option) any later version.
157*
158* This program is distributed in the hope that it will be
159* useful, but WITHOUT ANY WARRANTY; without even the implied
160* warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
161* PURPOSE. See the GNU General Public License for more details.
162*
163* You should have received a copy of the GNU General Public License
164* along with this program; if not, write to the Free Software
165* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
166* MA 02110-1301, USA.
167
168* Bugs:
169* {note_any_bugs_here}
170*-
171*/
172
173#include <math.h>
174
175#include "pal.h"
176#include "palmac.h"
177
178void palEl2ue ( double date, int jform, double epoch, double orbinc,
179 double anode, double perih, double aorq, double e,
180 double aorl, double dm, double u[13], int *jstat ) {
181
182 /* Sin and cos of J2000 mean obliquity (IAU 1976) */
183 const double SE=0.3977771559319137;
184 const double CE=0.9174820620691818;
185
186 int J;
187
188 double PHT,ARGPH,Q,W,CM,ALPHA,PHS,SW,CW,SI,CI,SO,CO,
189 X,Y,Z,PX,PY,PZ,VX,VY,VZ,DT,FC,FP,PSI,
190 UL[13],PV[6];
191
192 /* Validate arguments. */
193 if (jform < 1 || jform > 3) {
194 *jstat = -1;
195 return;
196 }
197 if (e < 0.0 || e > 10.0 || (e >= 1.0 && jform != 3)) {
198 *jstat = -2;
199 return;
200 }
201 if (aorq <= 0.0) {
202 *jstat = -3;
203 return;
204 }
205 if (jform == 1 && dm <= 0.0) {
206 *jstat = -4;
207 return;
208 }
209
210 /*
211 * Transform elements into standard form:
212 *
213 * PHT = epoch of perihelion passage
214 * ARGPH = argument of perihelion (little omega)
215 * Q = perihelion distance (q)
216 * CM = combined mass, M+m (mu)
217 */
218
219 if (jform == 1) {
220
221 /* Major planet. */
222 PHT = epoch-(aorl-perih)/dm;
223 ARGPH = perih-anode;
224 Q = aorq*(1.0-e);
225 W = dm/PAL__GCON;
226 CM = W*W*aorq*aorq*aorq;
227
228 } else if (jform == 2) {
229
230 /* Minor planet. */
231 PHT = epoch-aorl*sqrt(aorq*aorq*aorq)/PAL__GCON;
232 ARGPH = perih;
233 Q = aorq*(1.0-e);
234 CM = 1.0;
235
236 } else {
237
238 /* Comet. */
239 PHT = epoch;
240 ARGPH = perih;
241 Q = aorq;
242 CM = 1.0;
243
244 }
245
246 /* The universal variable alpha. This is proportional to the total
247 * energy of the orbit: -ve for an ellipse, zero for a parabola,
248 * +ve for a hyperbola. */
249
250 ALPHA = CM*(e-1.0)/Q;
251
252 /* Speed at perihelion. */
253
254 PHS = sqrt(ALPHA+2.0*CM/Q);
255
256 /* In a Cartesian coordinate system which has the x-axis pointing
257 * to perihelion and the z-axis normal to the orbit (such that the
258 * object orbits counter-clockwise as seen from +ve z), the
259 * perihelion position and velocity vectors are:
260 *
261 * position [Q,0,0]
262 * velocity [0,PHS,0]
263 *
264 * To express the results in J2000 equatorial coordinates we make a
265 * series of four rotations of the Cartesian axes:
266 *
267 * axis Euler angle
268 *
269 * 1 z argument of perihelion (little omega)
270 * 2 x inclination (i)
271 * 3 z longitude of the ascending node (big omega)
272 * 4 x J2000 obliquity (epsilon)
273 *
274 * In each case the rotation is clockwise as seen from the +ve end of
275 * the axis concerned.
276 */
277
278 /* Functions of the Euler angles. */
279 SW = sin(ARGPH);
280 CW = cos(ARGPH);
281 SI = sin(orbinc);
282 CI = cos(orbinc);
283 SO = sin(anode);
284 CO = cos(anode);
285
286 /* Position at perihelion (AU). */
287 X = Q*CW;
288 Y = Q*SW;
289 Z = Y*SI;
290 Y = Y*CI;
291 PX = X*CO-Y*SO;
292 Y = X*SO+Y*CO;
293 PY = Y*CE-Z*SE;
294 PZ = Y*SE+Z*CE;
295
296 /* Velocity at perihelion (AU per canonical day). */
297 X = -PHS*SW;
298 Y = PHS*CW;
299 Z = Y*SI;
300 Y = Y*CI;
301 VX = X*CO-Y*SO;
302 Y = X*SO+Y*CO;
303 VY = Y*CE-Z*SE;
304 VZ = Y*SE+Z*CE;
305
306 /* Time from perihelion to date (in Canonical Days: a canonical day
307 * is 58.1324409... days, defined as 1/PAL__GCON). */
308
309 DT = (date-PHT)*PAL__GCON;
310
311 /* First approximation to the Universal Eccentric Anomaly, PSI,
312 * based on the circle (FC) and parabola (FP) values. */
313
314 FC = DT/Q;
315 W = pow(3.0*DT+sqrt(9.0*DT*DT+8.0*Q*Q*Q), 1.0/3.0);
316 FP = W-2.0*Q/W;
317 PSI = (1.0-e)*FC+e*FP;
318
319 /* Assemble local copy of element set. */
320 UL[0] = CM;
321 UL[1] = ALPHA;
322 UL[2] = PHT;
323 UL[3] = PX;
324 UL[4] = PY;
325 UL[5] = PZ;
326 UL[6] = VX;
327 UL[7] = VY;
328 UL[8] = VZ;
329 UL[9] = Q;
330 UL[10] = 0.0;
331 UL[11] = date;
332 UL[12] = PSI;
333
334 /* Predict position+velocity at epoch of osculation. */
335 palUe2pv( date, UL, PV, &J );
336 if (J != 0) {
337 *jstat = -5;
338 return;
339 }
340
341 /* Convert back to universal elements. */
342 palPv2ue( PV, date, CM-1.0, u, &J );
343 if (J != 0) {
344 *jstat = -5;
345 return;
346 }
347
348 /* OK exit. */
349 *jstat = 0;
350
351}
Note: See TracBrowser for help on using the repository browser.