source: trunk/FACT++/pal/palPv2el.c@ 19901

Last change on this file since 19901 was 18347, checked in by tbretz, 9 years ago
File size: 11.8 KB
Line 
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
201void 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}
Note: See TracBrowser for help on using the repository browser.