source: trunk/MagicSoft/slalib/ue2pv.c

Last change on this file was 732, checked in by tbretz, 24 years ago
*** empty log message ***
  • Property svn:executable set to *
File size: 6.8 KB
Line 
1#include "slalib.h"
2#include "slamac.h"
3void slaUe2pv ( double date, double u[], double pv[], int *jstat )
4/*
5** - - - - - - - - -
6** s l a U e 2 p v
7** - - - - - - - - -
8**
9** Heliocentric position and velocity of a planet, asteroid or comet,
10** starting from orbital elements in the "universal variables" form.
11**
12** Given:
13** date double date, Modified Julian Date (JD-2400000.5)
14**
15** Given and returned:
16** u double[13] universal orbital elements (updated; Note 1)
17**
18** given [0] combined mass (M+m)
19** " [1] total energy of the orbit (alpha)
20** " [2] reference (osculating) epoch (t0)
21** " [3-5] position at reference epoch (r0)
22** " [6-8] velocity at reference epoch (v0)
23** " [9] heliocentric distance at reference epoch
24** " [10] r0.v0
25** returned [11] date (t)
26** " [12] universal eccentric anomaly (psi) of date, approx
27**
28** Returned:
29** pv double[6] position (AU) and velocity (AU/s)
30** jstat int* status: 0 = OK
31** -1 = radius vector zero
32** -2 = failed to converge
33**
34** Notes
35**
36** 1 The "universal" elements are those which define the orbit for the
37** purposes of the method of universal variables (see reference).
38** They consist of the combined mass of the two bodies, an epoch,
39** and the position and velocity vectors (arbitrary reference frame)
40** at that epoch. The parameter set used here includes also various
41** quantities that can, in fact, be derived from the other
42** information. This approach is taken to avoiding unnecessary
43** computation and loss of accuracy. The supplementary quantities
44** are (i) alpha, which is proportional to the total energy of the
45** orbit, (ii) the heliocentric distance at epoch, (iii) the
46** outwards component of the velocity at the given epoch, (iv) an
47** estimate of psi, the "universal eccentric anomaly" at a given
48** date and (v) that date.
49**
50** 2 The companion routine is slaEl2ue. This takes the conventional
51** orbital elements and transforms them into the set of numbers
52** needed by the present routine. A single prediction requires one
53** one call to slaEl2ue followed by one call to the present routine;
54** for convenience, the two calls are packaged as the routine
55** slaPlanel. Multiple predictions may be made by again calling
56** slaEl2ue once, but then calling the present routine multiple times,
57** which is faster than multiple calls to slaPlanel.
58**
59** It is not obligatory to use slaEl2ue to obtain the parameters.
60** However, it should be noted that because slaEl2ue performs its
61** own validation, no checks on the contents of the array U are made
62** by the present routine.
63**
64** 3 date is the instant for which the prediction is required. It is
65** in the TT timescale (formerly Ephemeris Time, ET) and is a
66** Modified Julian Date (JD-2400000.5).
67**
68** 4 The universal elements supplied in the array u are in canonical
69** units (solar masses, AU and canonical days). The position and
70** velocity are not sensitive to the choice of reference frame. The
71** slaEl2ue routine in fact produces coordinates with respect to the
72** J2000 equator and equinox.
73**
74** 5 The algorithm was originally adapted from the EPHSLA program of
75** D.H.P.Jones (private communication, 1996). The method is based
76** on Stumpff's Universal Variables.
77**
78** Reference: Everhart, E. & Pitkin, E.T., Am.J.Phys. 51, 712, 1983.
79**
80** Last revision: 19 March 1999
81**
82** Copyright P.T.Wallace. All rights reserved.
83*/
84
85/* Gaussian gravitational constant (exact) */
86#define GCON 0.01720209895
87
88/* Canonical days to seconds */
89#define CD2S ( GCON / 86400.0 )
90
91/* Test value for solution and maximum number of iterations */
92#define TEST 1e-13
93#define NITMAX 20
94
95{
96 int i, nit, n;
97 double cm, alpha, t0, p0[3], v0[3], r0, sigma0, t, psi, dt, w,
98 tol, psj, psj2, beta, s0, s1, s2, s3, ff, r, f, g, fd, gd;
99
100
101
102/* Unpack the parameters. */
103 cm = u[0];
104 alpha = u[1];
105 t0 = u[2];
106 for ( i = 0; i < 3; i++ ) {
107 p0[i] = u[i+3];
108 v0[i] = u[i+6];
109 }
110 r0 = u[9];
111 sigma0 = u[10];
112 t = u[11];
113 psi = u[12];
114
115/* Approximately update the universal eccentric anomaly. */
116 psi = psi + ( date - t ) * GCON / r0;
117
118/* Time from reference epoch to date (in Canonical Days: a canonical */
119/* day is 58.1324409... days, defined as 1/GCON). */
120 dt = ( date - t0 ) * GCON;
121
122/* Refine the universal eccentric anomaly. */
123 nit = 1;
124 w = 1.0;
125 tol = 0.0;
126 while ( fabs ( w ) >= tol ) {
127
128 /* Form half angles until BETA small enough. */
129 n = 0;
130 psj = psi;
131 psj2 = psj * psj;
132 beta = alpha * psj2;
133 while ( fabs ( beta ) > 0.7 ) {
134 n++;
135 beta /= 4.0;
136 psj /= 2.0;
137 psj2 /= 4.0;
138 }
139
140 /* Calculate Universal Variables S0,S1,S2,S3 by nested series. */
141 s3 = psj * psj2 * ( ( ( ( ( ( beta / 210.0 + 1.0 )
142 * beta / 156.0 + 1.0 )
143 * beta / 110.0 + 1.0 )
144 * beta / 72.0 + 1.0 )
145 * beta / 42.0 + 1.0 )
146 * beta / 20.0 + 1.0 ) / 6.0;
147 s2 = psj2 * ( ( ( ( ( ( beta / 182.0 + 1.0 )
148 * beta / 132.0 + 1.0 )
149 * beta / 90.0 + 1.0 )
150 * beta / 56.0 + 1.0 )
151 * beta / 30.0 + 1.0 )
152 * beta / 12.0 + 1.0 ) / 2.0;
153 s1 = psj + alpha * s3;
154 s0 = 1.0 + alpha * s2;
155
156 /* Undo the angle-halving. */
157 tol = TEST;
158 while ( n > 0 ) {
159 s3 = 2.0 * ( s0 * s3 + psj * s2 );
160 s2 = 2.0 * s1 * s1;
161 s1 = 2.0 * s0 * s1;
162 s0 = 2.0 * s0 * s0 - 1.0;
163 psj += psj;
164 tol += tol;
165 n--;
166 }
167
168 /* Improve the approximation to psi. */
169 ff = r0 * s1 + sigma0 * s2 + cm * s3 - dt;
170 r = r0 * s0 + sigma0 * s1 + cm * s2;
171 if ( r == 0.0 ) {
172 *jstat = -1;
173 return;
174 }
175 w = ff / r;
176 psi -= w;
177
178 /* Next iteration, unless too many already. */
179 if ( nit >= NITMAX ) {
180 *jstat = -2;
181 return;
182 }
183 nit++;
184 }
185
186/* Project the position and velocity vectors (scaling velocity to AU/s). */
187 w = cm * s2;
188 f = 1.0 - w / r0;
189 g = dt - cm * s3;
190 fd = - cm * s1 / ( r0 * r );
191 gd = 1.0 - w / r;
192 for ( i = 0; i < 3; i++ ) {
193 pv[i] = p0[i] * f + v0[i] * g;
194 pv[i+3] = CD2S * ( p0[i] * fd + v0[i] * gd );
195 }
196
197/* Update the parameters to allow speedy prediction of psi next time. */
198 u[11] = date;
199 u[12] = psi;
200
201/* OK exit. */
202 *jstat = 0;
203
204}
Note: See TracBrowser for help on using the repository browser.