1 /*
2 *+
3 * Name:
4 * palPv2el
5
6 * Purpose:
7 * Position velocity to heliocentirc osculating elements
8
9 * Language:
10 * Starlink ANSI C
11
12 * Type of Module:
13 * Library routine
14
15 * Invocation:
16 * void palPv2el ( const double pv[6], double date, double pmass, 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 * pv = const double [6] (Given)
23 * Heliocentric x,y,z,xdot,ydot,zdot of date,
24 * J2000 equatorial triad (AU,AU/s; Note 1)
25 * date = double (Given)
26 * Date (TT Modified Julian Date = JD-2400000.5)
27 * pmass = double (Given)
28 * Mass of the planet (Sun=1; Note 2)
29 * jformr = int (Given)
30 * Requested element set (1-3; Note 3)
31 * jform = int * (Returned)
32 * Element set actually returned (1-3; Note 4)
33 * epoch = double * (Returned)
34 * Epoch of elements (TT MJD)
35 * orbinc = double * (Returned)
36 * inclination (radians)
37 * anode = double * (Returned)
38 * longitude of the ascending node (radians)
39 * perih = double * (Returned)
40 * longitude or argument of perihelion (radians)
41 * aorq = double * (Returned)
42 * mean distance or perihelion distance (AU)
43 * e = double * (Returned)
44 * eccentricity
45 * aorl = double * (Returned)
46 * mean anomaly or longitude (radians, JFORM=1,2 only)
47 * dm = double * (Returned)
48 * daily motion (radians, JFORM=1 only)
49 * jstat = int * (Returned)
50 * status: 0 = OK
51 * - -1 = illegal PMASS
52 * - -2 = illegal JFORMR
53 * - -3 = position/velocity out of range
54
55 * Description:
56 * Heliocentric osculating elements obtained from instantaneous position
57 * and velocity.
58
59 * Authors:
60 * PTW: Pat Wallace (STFC)
61 * TIMJ: Tim Jenness (JAC, Hawaii)
62 * {enter_new_authors_here}
63
64 * Notes:
65 * - The PV 6-vector is with respect to the mean equator and equinox of
66 * epoch J2000. The orbital elements produced are with respect to
67 * the J2000 ecliptic and mean equinox.
68 * - The mass, PMASS, is important only for the larger planets. For
69 * most purposes (e.g. asteroids) use 0D0. Values less than zero
70 * are illegal.
71 * - Three different element-format options are supported:
72 *
73 * Option JFORM=1, suitable for the major planets:
74 *
75 * EPOCH = epoch of elements (TT MJD)
76 * ORBINC = inclination i (radians)
77 * ANODE = longitude of the ascending node, big omega (radians)
78 * PERIH = longitude of perihelion, curly pi (radians)
79 * AORQ = mean distance, a (AU)
80 * E = eccentricity, e
81 * AORL = mean longitude L (radians)
82 * DM = daily motion (radians)
83 *
84 * Option JFORM=2, suitable for minor planets:
85 *
86 * EPOCH = epoch of elements (TT MJD)
87 * ORBINC = inclination i (radians)
88 * ANODE = longitude of the ascending node, big omega (radians)
89 * PERIH = argument of perihelion, little omega (radians)
90 * AORQ = mean distance, a (AU)
91 * E = eccentricity, e
92 * AORL = mean anomaly M (radians)
93 *
94 * Option JFORM=3, suitable for comets:
95 *
96 * EPOCH = epoch of perihelion (TT MJD)
97 * ORBINC = inclination i (radians)
98 * ANODE = longitude of the ascending node, big omega (radians)
99 * PERIH = argument of perihelion, little omega (radians)
100 * AORQ = perihelion distance, q (AU)
101 * E = eccentricity, e
102 *
103 * - It may not be possible to generate elements in the form
104 * requested through JFORMR. The caller is notified of the form
105 * of elements actually returned by means of the JFORM argument:
106
107 * JFORMR JFORM meaning
108 *
109 * 1 1 OK - elements are in the requested format
110 * 1 2 never happens
111 * 1 3 orbit not elliptical
112 *
113 * 2 1 never happens
114 * 2 2 OK - elements are in the requested format
115 * 2 3 orbit not elliptical
116 *
117 * 3 1 never happens
118 * 3 2 never happens
119 * 3 3 OK - elements are in the requested format
120 *
121 * - The arguments returned for each value of JFORM (cf Note 5: JFORM
122 * may not be the same as JFORMR) are as follows:
123 *
124 * JFORM 1 2 3
125 * EPOCH t0 t0 T
126 * ORBINC i i i
127 * ANODE Omega Omega Omega
128 * PERIH curly pi omega omega
129 * AORQ a a q
130 * E e e e
131 * AORL L M -
132 * DM n - -
133 *
134 * where:
135 *
136 * t0 is the epoch of the elements (MJD, TT)
137 * T " epoch of perihelion (MJD, TT)
138 * i " inclination (radians)
139 * Omega " longitude of the ascending node (radians)
140 * curly pi " longitude of perihelion (radians)
141 * omega " argument of perihelion (radians)
142 * a " mean distance (AU)
143 * q " perihelion distance (AU)
144 * e " eccentricity
145 * L " longitude (radians, 0-2pi)
146 * M " mean anomaly (radians, 0-2pi)
147 * n " daily motion (radians)
148 * - means no value is set
149 *
150 * - At very small inclinations, the longitude of the ascending node
151 * ANODE becomes indeterminate and under some circumstances may be
152 * set arbitrarily to zero. Similarly, if the orbit is close to
153 * circular, the true anomaly becomes indeterminate and under some
154 * circumstances may be set arbitrarily to zero. In such cases,
155 * the other elements are automatically adjusted to compensate,
156 * and so the elements remain a valid description of the orbit.
157 * - The osculating epoch for the returned elements is the argument
158 * DATE.
159 *
160 * - Reference: Sterne, Theodore E., "An Introduction to Celestial
161 * Mechanics", Interscience Publishers, 1960
162
163 * History:
164 * 2012-03-09 (TIMJ):
165 * Initial version converted from SLA/F.
166 * Adapted with permission from the Fortran SLALIB library.
167 * {enter_further_changes_here}
168
169 * Copyright:
170 * Copyright (C) 2005 Patrick T. Wallace
171 * Copyright (C) 2012 Science and Technology Facilities Council.
172 * All Rights Reserved.
173
174 * Licence:
175 * This program is free software; you can redistribute it and/or
176 * modify it under the terms of the GNU General Public License as
177 * published by the Free Software Foundation; either version 3 of
178 * the License, or (at your option) any later version.
179 *
180 * This program is distributed in the hope that it will be
181 * useful, but WITHOUT ANY WARRANTY; without even the implied
182 * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
183 * PURPOSE. See the GNU General Public License for more details.
184 *
185 * You should have received a copy of the GNU General Public License
186 * along with this program; if not, write to the Free Software
187 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
188 * MA 02110-1301, USA.
189
190 * Bugs:
191 * {note_any_bugs_here}
192 *-
193 */
194
195 #include <math.h>
196
197 #include "pal1sofa.h"
198 #include "pal.h"
199 #include "palmac.h"
200
palPv2el(const double pv[6],double date,double pmass,int jformr,int * jform,double * epoch,double * orbinc,double * anode,double * perih,double * aorq,double * e,double * aorl,double * dm,int * jstat)201 void palPv2el ( const double pv[6], double date, double pmass, int jformr,
202 int *jform, double *epoch, double *orbinc,
203 double *anode, double *perih, double *aorq, double *e,
204 double *aorl, double *dm, int *jstat ) {
205
206 /* Sin and cos of J2000 mean obliquity (IAU 1976) */
207 const double SE = 0.3977771559319137;
208 const double CE = 0.9174820620691818;
209
210 /* Minimum allowed distance (AU) and speed (AU/day) */
211 const double RMIN = 1e-3;
212 const double VMIN = 1e-8;
213
214 /* How close to unity the eccentricity has to be to call it a parabola */
215 const double PARAB = 1.0e-8;
216
217 double X,Y,Z,XD,YD,ZD,R,V2,V,RDV,GMU,HX,HY,HZ,
218 HX2PY2,H2,H,OI,BIGOM,AR,ECC,S,C,AT,U,OM,
219 GAR3,EM1,EP1,HAT,SHAT,CHAT,AE,AM,DN,PL,
220 EL,Q,TP,THAT,THHF,F;
221
222 int JF;
223
224 /* Validate arguments PMASS and JFORMR.*/
225 if (pmass < 0.0) {
226 *jstat = -1;
227 return;
228 }
229 if (jformr < 1 || jformr > 3) {
230 *jstat = -2;
231 return;
232 }
233
234 /* Provisionally assume the elements will be in the chosen form. */
235 JF = jformr;
236
237 /* Rotate the position from equatorial to ecliptic coordinates. */
238 X = pv[0];
239 Y = pv[1]*CE+pv[2]*SE;
240 Z = -pv[1]*SE+pv[2]*CE;
241
242 /* Rotate the velocity similarly, scaling to AU/day. */
243 XD = PAL__SPD*pv[3];
244 YD = PAL__SPD*(pv[4]*CE+pv[5]*SE);
245 ZD = PAL__SPD*(-pv[4]*SE+pv[5]*CE);
246
247 /* Distance and speed. */
248 R = sqrt(X*X+Y*Y+Z*Z);
249 V2 = XD*XD+YD*YD+ZD*ZD;
250 V = sqrt(V2);
251
252 /* Reject unreasonably small values. */
253 if (R < RMIN || V < VMIN) {
254 *jstat = -3;
255 return;
256 }
257
258 /* R dot V. */
259 RDV = X*XD+Y*YD+Z*ZD;
260
261 /* Mu. */
262 GMU = (1.0+pmass)*PAL__GCON*PAL__GCON;
263
264 /* Vector angular momentum per unit reduced mass. */
265 HX = Y*ZD-Z*YD;
266 HY = Z*XD-X*ZD;
267 HZ = X*YD-Y*XD;
268
269 /* Areal constant. */
270 HX2PY2 = HX*HX+HY*HY;
271 H2 = HX2PY2+HZ*HZ;
272 H = sqrt(H2);
273
274 /* Inclination. */
275 OI = atan2(sqrt(HX2PY2),HZ);
276
277 /* Longitude of ascending node. */
278 if (HX != 0.0 || HY != 0.0) {
279 BIGOM = atan2(HX,-HY);
280 } else {
281 BIGOM=0.0;
282 }
283
284 /* Reciprocal of mean distance etc. */
285 AR = 2.0/R-V2/GMU;
286
287 /* Eccentricity. */
288 ECC = sqrt(DMAX(1.0-AR*H2/GMU,0.0));
289
290 /* True anomaly. */
291 S = H*RDV;
292 C = H2-R*GMU;
293 if (S != 0.0 || C != 0.0) {
294 AT = atan2(S,C);
295 } else {
296 AT = 0.0;
297 }
298
299 /* Argument of the latitude. */
300 S = sin(BIGOM);
301 C = cos(BIGOM);
302 U = atan2((-X*S+Y*C)*cos(OI)+Z*sin(OI),X*C+Y*S);
303
304 /* Argument of perihelion. */
305 OM = U-AT;
306
307 /* Capture near-parabolic cases. */
308 if (fabs(ECC-1.0) < PARAB) ECC=1.0;
309
310 /* Comply with JFORMR = 1 or 2 only if orbit is elliptical. */
311 if (ECC > 1.0) JF=3;
312
313 /* Functions. */
314 GAR3 = GMU*AR*AR*AR;
315 EM1 = ECC-1.0;
316 EP1 = ECC+1.0;
317 HAT = AT/2.0;
318 SHAT = sin(HAT);
319 CHAT = cos(HAT);
320
321 /* Variable initializations to avoid compiler warnings. */
322 AM = 0.0;
323 DN = 0.0;
324 PL = 0.0;
325 EL = 0.0;
326 Q = 0.0;
327 TP = 0.0;
328
329 /* Ellipse? */
330 if (ECC < 1.0 ) {
331
332 /* Eccentric anomaly. */
333 AE = 2.0*atan2(sqrt(-EM1)*SHAT,sqrt(EP1)*CHAT);
334
335 /* Mean anomaly. */
336 AM = AE-ECC*sin(AE);
337
338 /* Daily motion. */
339 DN = sqrt(GAR3);
340 }
341
342 /* "Major planet" element set? */
343 if (JF == 1) {
344
345 /* Longitude of perihelion. */
346 PL = BIGOM+OM;
347
348 /* Longitude at epoch. */
349 EL = PL+AM;
350 }
351
352 /* "Comet" element set? */
353 if (JF == 3) {
354
355 /* Perihelion distance. */
356 Q = H2/(GMU*EP1);
357
358 /* Ellipse, parabola, hyperbola? */
359 if (ECC < 1.0) {
360
361 /* Ellipse: epoch of perihelion. */
362 TP = date-AM/DN;
363
364 } else {
365
366 /* Parabola or hyperbola: evaluate tan ( ( true anomaly ) / 2 ) */
367 THAT = SHAT/CHAT;
368 if (ECC == 1.0) {
369
370 /* Parabola: epoch of perihelion. */
371 TP = date-THAT*(1.0+THAT*THAT/3.0)*H*H2/(2.0*GMU*GMU);
372
373 } else {
374
375 /* Hyperbola: epoch of perihelion. */
376 THHF = sqrt(EM1/EP1)*THAT;
377 F = log(1.0+THHF)-log(1.0-THHF);
378 TP = date-(ECC*sinh(F)-F)/sqrt(-GAR3);
379 }
380 }
381 }
382
383 /* Return the appropriate set of elements. */
384 *jform = JF;
385 *orbinc = OI;
386 *anode = eraAnp(BIGOM);
387 *e = ECC;
388 if (JF == 1) {
389 *perih = eraAnp(PL);
390 *aorl = eraAnp(EL);
391 *dm = DN;
392 } else {
393 *perih = eraAnp(OM);
394 if (JF == 2) *aorl = eraAnp(AM);
395 }
396 if (JF != 3) {
397 *epoch = date;
398 *aorq = 1.0/AR;
399 } else {
400 *epoch = TP;
401 *aorq = Q;
402 }
403 *jstat = 0;
404
405 }
406