FACT++  1.0
palEl2ue.c
Go to the documentation of this file.
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 
178 void 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 }
void palPv2ue(const double pv[6], double date, double pmass, double u[13], int *jstat)
Definition: palPv2ue.c:114
void palEl2ue(double date, int jform, double epoch, double orbinc, double anode, double perih, double aorq, double e, double aorl, double dm, double u[13], int *jstat)
Definition: palEl2ue.c:178
static const double PAL__GCON
Definition: palmac.h:117
void palUe2pv(double date, double u[13], double pv[], int *jstat)