FACT++  1.0
palPv2ue.c
Go to the documentation of this file.
1 /*
2 *+
3 * Name:
4 * palPv2ue
5 
6 * Purpose:
7 * Universal elements to position and velocity.
8 
9 * Language:
10 * Starlink ANSI C
11 
12 * Type of Module:
13 * Library routine
14 
15 * Invocation:
16 * void palPv2ue( const double pv[6], double date, double pmass,
17 * double u[13], int * jstat );
18 
19 * Arguments:
20 * pv = double [6] (Given)
21 * Heliocentric x,y,z,xdot,ydot,zdot of date, (AU,AU/s; Note 1)
22 * date = double (Given)
23 * Date (TT modified Julian Date = JD-2400000.5)
24 * pmass = double (Given)
25 * Mass of the planet (Sun=1; note 2)
26 * u = double [13] (Returned)
27 * Universal orbital elements (Note 3)
28 *
29 * - (0) combined mass (M+m)
30 * - (1) total energy of the orbit (alpha)
31 * - (2) reference (osculating) epoch (t0)
32 * - (3-5) position at reference epoch (r0)
33 * - (6-8) velocity at reference epoch (v0)
34 * - (9) heliocentric distance at reference epoch
35 * - (10) r0.v0
36 * - (11) date (t)
37 * - (12) universal eccentric anomaly (psi) of date, approx
38 * jstat = int * (Returned)
39 * status: 0 = OK
40 * - -1 = illegal PMASS
41 * - -2 = too close to Sun
42 * - -3 = too slow
43 
44 * Description:
45 * Construct a universal element set based on an instantaneous position
46 * and velocity.
47 
48 
49 * Authors:
50 * PTW: Pat Wallace (STFC)
51 * TIMJ: Tim Jenness (JAC, Hawaii)
52 * {enter_new_authors_here}
53 
54 * Notes:
55 * - The PV 6-vector can be with respect to any chosen inertial frame,
56 * and the resulting universal-element set will be with respect to
57 * the same frame. A common choice will be mean equator and ecliptic
58 * of epoch J2000.
59 * - The mass, PMASS, is important only for the larger planets. For
60 * most purposes (e.g. asteroids) use 0D0. Values less than zero
61 * are illegal.
62 * - The "universal" elements are those which define the orbit for the
63 * purposes of the method of universal variables (see reference).
64 * They consist of the combined mass of the two bodies, an epoch,
65 * and the position and velocity vectors (arbitrary reference frame)
66 * at that epoch. The parameter set used here includes also various
67 * quantities that can, in fact, be derived from the other
68 * information. This approach is taken to avoiding unnecessary
69 * computation and loss of accuracy. The supplementary quantities
70 * are (i) alpha, which is proportional to the total energy of the
71 * orbit, (ii) the heliocentric distance at epoch, (iii) the
72 * outwards component of the velocity at the given epoch, (iv) an
73 * estimate of psi, the "universal eccentric anomaly" at a given
74 * date and (v) that date.
75 * - Reference: Everhart, E. & Pitkin, E.T., Am.J.Phys. 51, 712, 1983.
76 
77 * History:
78 * 2012-03-09 (TIMJ):
79 * Initial version from the SLA/F implementation.
80 * Adapted with permission from the Fortran SLALIB library.
81 * {enter_further_changes_here}
82 
83 * Copyright:
84 * Copyright (C) 1999 Rutherford Appleton Laboratory
85 * Copyright (C) 2012 Science and Technology Facilities Council.
86 * All Rights Reserved.
87 
88 * Licence:
89 * This program is free software; you can redistribute it and/or
90 * modify it under the terms of the GNU General Public License as
91 * published by the Free Software Foundation; either version 3 of
92 * the License, or (at your option) any later version.
93 *
94 * This program is distributed in the hope that it will be
95 * useful, but WITHOUT ANY WARRANTY; without even the implied
96 * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
97 * PURPOSE. See the GNU General Public License for more details.
98 *
99 * You should have received a copy of the GNU General Public License
100 * along with this program; if not, write to the Free Software
101 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
102 * MA 02110-1301, USA.
103 
104 * Bugs:
105 * {note_any_bugs_here}
106 *-
107 */
108 
109 #include <math.h>
110 
111 #include "pal.h"
112 #include "palmac.h"
113 
114 void palPv2ue( const double pv[6], double date, double pmass,
115  double u[13], int * jstat ) {
116 
117  /* Canonical days to seconds */
118  const double CD2S = PAL__GCON / PAL__SPD;
119 
120  /* Minimum allowed distance (AU) and speed (AU per canonical day) */
121  const double RMIN = 1e-3;
122  const double VMIN = 1e-3;
123 
124  double T0,CM,X,Y,Z,XD,YD,ZD,R,V2,V,ALPHA,RDV;
125 
126  /* Reference epoch. */
127  T0 = date;
128 
129  /* Combined mass (mu=M+m). */
130  if (pmass < 0.0 ) { /* Negative planet mass */
131  *jstat = -1;
132  return;
133  }
134  CM = 1.0+pmass;
135 
136  /* Unpack the state vector, expressing velocity in AU per canonical day. */
137  X = pv[0];
138  Y = pv[1];
139  Z = pv[2];
140  XD = pv[3]/CD2S;
141  YD = pv[4]/CD2S;
142  ZD = pv[5]/CD2S;
143 
144  /* Heliocentric distance, and speed. */
145  R = sqrt(X*X+Y*Y+Z*Z);
146  V2 = XD*XD+YD*YD+ZD*ZD;
147  V = sqrt(V2);
148 
149  /* Reject unreasonably small values. */
150  if (R < RMIN) { /* Too close */
151  *jstat = -2;
152  return;
153  }
154  if (V < VMIN) { /* Too slow */
155  *jstat = -3;
156  return;
157  }
158 
159  /* Total energy of the orbit. */
160  ALPHA = V2-2.0*CM/R;
161 
162  /* Outward component of velocity. */
163  RDV = X*XD+Y*YD+Z*ZD;
164 
165  /* Construct the universal-element set. */
166  u[0] = CM;
167  u[1] = ALPHA;
168  u[2] = T0;
169  u[3] = X;
170  u[4] = Y;
171  u[5] = Z;
172  u[6] = XD;
173  u[7] = YD;
174  u[8] = ZD;
175  u[9] = R;
176  u[10] = RDV;
177  u[11] = T0;
178  u[12] = 0.0;
179 
180  *jstat = 0;
181  return;
182 }
static const double PAL__GCON
Definition: palmac.h:117
void palPv2ue(const double pv[6], double date, double pmass, double u[13], int *jstat)
Definition: palPv2ue.c:114
static const double PAL__SPD
Definition: palmac.h:102