FACT++  1.0
palPv2el.c
Go to the documentation of this file.
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 
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 }
static const double PAL__GCON
Definition: palmac.h:117
in C
Definition: README_v19.txt:379
void 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)
Definition: palPv2el.c:201
#define DMAX(A, B)
Definition: palmac.h:126
double eraAnp(double a)
Definition: anp.c:3
static const double PAL__SPD
Definition: palmac.h:102