1 /*
2 * - - - - - - - - - - -
3 * g a l _ p l a n 9 4
4 * - - - - - - - - - - -
5 *
6 * This routine is part of the General Astrodynamics Library
7 *
8 * Description:
9 *
10 * This routine is an independent translation of a FORTRAN routine
11 * that is part of IAU's SOFA software collection.
12 *
13 * Status:
14 *
15 * support routine.
16 *
17 * Approximate heliocentric position and velocity of a nominated major
18 * planet: Mercury, Venus, EMB, Mars, Jupiter, Saturn, Uranus or
19 * Neptune (but not the Earth itself).
20 *
21 * Given:
22 *
23 * date1 d TDB date part A (Note 1)
24 * date2 d TDB date part B (Note 1)
25 * np i planet (1=Mercury, 2=Venus, 3=EMB ... 8=Neptune)
26 *
27 * Returned:
28 *
29 * pv d[2][3] planet pos,vel (heliocentric, J2000, AU, AU/d)
30 * gal_plan94 i status: -1 = illegal NP (outside 1-8)
31 * 0 = OK
32 * +1 = warning: date outside 1000-3000 AD
33 * +2 = warning: solution failed to converge
34 *
35 * Notes:
36 *
37 * 1) The TT date date1+date2 is a Julian Date, apportioned in any
38 * convenient way between the two arguments. For example,
39 * JD(TT)=2450123.7 could be expressed in any of these ways,
40 * among others:
41 *
42 * date1 date2
43 *
44 * 2450123.7 0.0 (JD method)
45 * 2451545.0 -1421.3 (J2000 method)
46 * 2400000.5 50123.2 (MJD method)
47 * 2450123.5 0.2 (date & time method)
48 *
49 * The JD method is the most natural and convenient to use in
50 * cases where the loss of several decimal digits of resolution
51 * is acceptable. The J2000 method is best matched to the way
52 * the argument is handled internally and will deliver the
53 * optimum resolution. The MJD method and the date & time methods
54 * are both good compromises between resolution and convenience.
55 * The limited accuracy of the present algorithm is such that any
56 * of the methods is satisfactory.
57 *
58 * 2) If an np value outside the range 1-8 is supplied, an error
59 * status (j = -1) is returned and the pv vector set to zeroes.
60 *
61 * 3) For np=3 the result is for the Earth-Moon Barycenter. To
62 * obtain the heliocentric position and velocity of the Earth,
63 * use instead the SOFA routine gal_epv00.
64 *
65 * 4) On successful return, the array PV contains the following:
66 *
67 * pv[0][0] x }
68 * pv[0][1] y } heliocentric position, AU
69 * pv[0][2] z }
70 *
71 * pv[1][0] xdot }
72 * pv[1][1] ydot } heliocentric velocity, AU/d
73 * pv[1][2] zdot }
74 *
75 * The reference frame is equatorial and is with respect to the
76 * mean equator and equinox of epoch J2000.
77 *
78 * 5) The algorithm is due to J.L. Simon, P. Bretagnon, J. Chapront,
79 * M. Chapront-Touze, G. Francou and J. Laskar (Bureau des
80 * Longitudes, Paris, France). From comparisons with JPL
81 * ephemeris DE102, they quote the following maximum errors
82 * over the interval 1800-2050:
83 *
84 * L (arcsec) B (arcsec) R (km)
85 *
86 * Mercury 4 1 300
87 * Venus 5 1 800
88 * EMB 6 1 1000
89 * Mars 17 1 7700
90 * Jupiter 71 5 76000
91 * Saturn 81 13 267000
92 * Uranus 86 7 712000
93 * Neptune 11 1 253000
94 *
95 * Over the interval 1000-3000, they report that the accuracy is no
96 * worse than 1.5 times that over 1800-2050. Outside 1000-3000 the
97 * accuracy declines.
98 *
99 * Comparisons of the present routine with the JPL DE200 ephemeris
100 * give the following RMS errors over the interval 1960-2025:
101 *
102 * position (km) velocity (m/s)
103 *
104 * Mercury 334 0.437
105 * Venus 1060 0.855
106 * EMB 2010 0.815
107 * Mars 7690 1.98
108 * Jupiter 71700 7.70
109 * Saturn 199000 19.4
110 * Uranus 564000 16.4
111 * Neptune 158000 14.4
112 *
113 * Comparisons against DE200 over the interval 1800-2100 gave the
114 * following maximum absolute differences. (The results using
115 * DE406 were essentially the same.)
116 *
117 * L (arcsec) B (arcsec) R (km) Rdot (m/s)
118 *
119 * Mercury 7 1 500 0.7
120 * Venus 7 1 1100 0.9
121 * EMB 9 1 1300 1.0
122 * Mars 26 1 9000 2.5
123 * Jupiter 78 6 82000 8.2
124 * Saturn 87 14 263000 24.6
125 * Uranus 86 7 661000 27.4
126 * Neptune 11 2 248000 21.4
127 *
128 * 6) The present SOFA re-implementation of the original Simon et al.
129 * Fortran code differs from the original in the following respects:
130 *
131 * * The date is supplied in two parts.
132 *
133 * * The result is returned only in equatorial Cartesian form;
134 * the ecliptic longitude, latitude and radius vector are not
135 * returned.
136 *
137 * * The result is in the J2000 equatorial frame, not ecliptic.
138 *
139 * * More is done in-line: there are fewer calls to other
140 * routines.
141 *
142 * * Different error/warning status values are used.
143 *
144 * * A different Kepler's-equation-solver is used (avoiding
145 * use of COMPLEX*16).
146 *
147 * * Polynomials in T are nested to minimize rounding errors.
148 *
149 * * Explicit double-precision constants are used to avoid
150 * mixed-mode expressions.
151 *
152 * * There are other, cosmetic, changes to comply with SOFA
153 * style conventions.
154 *
155 * None of the above changes affects the result significantly.
156 *
157 * 7) The returned status, j, indicates the most serious condition
158 * encountered during execution of the routine. Illegal np is
159 * considered the most serious, overriding failure to converge,
160 * which in turn takes precedence over the remote epoch warning.
161 *
162 * Called:
163 *
164 * gal_anp normalize angle into range 0 to 2pi
165 *
166 * References: Simon, J.L, Bretagnon, P., Chapront, J.,
167 * Chapront-Touze, M., Francou, G., and Laskar, J.,
168 * Astron. Astrophys. 282, 663 (1994).
169 *
170 * This revision:
171 *
172 * 2006 November 13 ( c version 2008 June 8 )
173 *
174 *
175 * Copyright (C) 2008 Paul C. L. Willmott. See notes at end.
176 *
177 *-----------------------------------------------------------------------
178 */
179
180 #include <math.h>
181 #include "gal_const.h"
182 #include "gal_plan94.h"
183 #include "gal_anpm.h"
184
185 int
gal_plan94(double date1,double date2,int np,double pv[2][3])186 gal_plan94
187 (
188 double date1,
189 double date2,
190 int np,
191 double pv[2][3] )
192 {
193
194 /*
195 * Maximum number of iterations allowed to solve Kepler's equation
196 */
197
198 const int KMAX = 10 ;
199
200 /*
201 * Sin and cos of J2000 mean obliquity (IAU 1976)
202 */
203
204 const double SINEPS = 0.3977771559319137 ;
205 const double COSEPS = 0.9174820620691818 ;
206
207 /*
208 * Gaussian constant
209 */
210
211 const double GK = 0.017202098950 ;
212
213 int jstat, k, i ;
214
215 double t, da, dl, de, dp, di, dom, dmu, arga, argl, am,
216 ae, dae, ae2, at, r, v, si2, xq, xp, tl, xsw,
217 xcw, xm2, xf, ci2, xms, xmc, xpxq2, x, y, z ;
218
219 /*
220 * Planetary inverse masses
221 */
222
223 static const double AMAS[8] = {
224 6023600.0,
225 408523.5,
226 328900.5,
227 3098710.0,
228 1047.355,
229 3498.5,
230 22869.0,
231 19314.0
232 } ;
233
234 /*
235 * Tables giving the mean Keplerian elements, limited to T**2 terms:
236 *
237 * A semi-major axis (AU)
238 * DLM mean longitude (degree and arcsecond)
239 * E eccentricity
240 * PI longitude of the perihelion (degree and arcsecond)
241 * DINC inclination (degree and arcsecond)
242 * OMEGA longitude of the ascending node (degree and arcsecond)
243 */
244
245 static const double A[8][3] = {
246
247 { 0.3870983098e0, 0e0, 0e0, } ,
248 { 0.7233298200e0, 0e0, 0e0, } ,
249 { 1.0000010178e0, 0e0, 0e0, } ,
250 { 1.5236793419e0, 3e-10, 0e0, } ,
251 { 5.2026032092e0, 19132e-10, -39e-10, } ,
252 { 9.5549091915e0, -0.0000213896e0, 444e-10, } ,
253 { 19.2184460618e0, -3716e-10, 979e-10, } ,
254 { 30.1103868694e0, -16635e-10, 686e-10, } ,
255
256 } ;
257
258 static const double DLM[8][3] = {
259
260 { 252.25090552e0, 5381016286.88982e0, -1.92789e0, } ,
261 { 181.97980085e0, 2106641364.33548e0, 0.59381e0, } ,
262 { 100.46645683e0, 1295977422.83429e0, -2.04411e0, } ,
263 { 355.43299958e0, 689050774.93988e0, 0.94264e0, } ,
264 { 34.35151874e0, 109256603.77991e0, -30.60378e0, } ,
265 { 50.07744430e0, 43996098.55732e0, 75.61614e0, } ,
266 { 314.05500511e0, 15424811.93933e0, -1.75083e0, } ,
267 { 304.34866548e0, 7865503.20744e0, 0.21103e0, } ,
268
269 } ;
270
271 static const double E[8][3] = {
272
273 { 0.2056317526e0, 0.0002040653e0, -28349e-10, } ,
274 { 0.0067719164e0, -0.0004776521e0, 98127e-10, } ,
275 { 0.0167086342e0, -0.0004203654e0, -0.0000126734e0, } ,
276 { 0.0934006477e0, 0.0009048438e0, -80641e-10, } ,
277 { 0.0484979255e0, 0.0016322542e0, -0.0000471366e0, } ,
278 { 0.0555481426e0, -0.0034664062e0, -0.0000643639e0, } ,
279 { 0.0463812221e0, -0.0002729293e0, 0.0000078913e0, } ,
280 { 0.0094557470e0, 0.0000603263e0, 0e0, } ,
281
282 } ;
283
284 static const double PI[8][3] = {
285
286 { 77.45611904e0, 5719.11590e0, -4.83016e0, } ,
287 { 131.56370300e0, 175.48640e0, -498.48184e0, } ,
288 { 102.93734808e0, 11612.35290e0, 53.27577e0, } ,
289 { 336.06023395e0, 15980.45908e0, -62.32800e0, } ,
290 { 14.33120687e0, 7758.75163e0, 259.95938e0, } ,
291 { 93.05723748e0, 20395.49439e0, 190.25952e0, } ,
292 { 173.00529106e0, 3215.56238e0, -34.09288e0, } ,
293 { 48.12027554e0, 1050.71912e0, 27.39717e0, } ,
294
295 } ;
296
297 static const double DINC[8][3] = {
298
299 { 7.00498625e0, -214.25629e0, 0.28977e0, } ,
300 { 3.39466189e0, -30.84437e0, -11.67836e0, } ,
301 { 0e0, 469.97289e0, -3.35053e0, } ,
302 { 1.84972648e0, -293.31722e0, -8.11830e0, } ,
303 { 1.30326698e0, -71.55890e0, 11.95297e0, } ,
304 { 2.48887878e0, 91.85195e0, -17.66225e0, } ,
305 { 0.77319689e0, -60.72723e0, 1.25759e0, } ,
306 { 1.76995259e0, 8.12333e0, 0.08135e0, } ,
307
308 } ;
309
310 static const double OMEGA[8][3] = {
311
312 { 48.33089304e0, -4515.21727e0, -31.79892e0, } ,
313 { 76.67992019e0, -10008.48154e0, -51.32614e0, } ,
314 { 174.87317577e0, -8679.27034e0, 15.34191e0, } ,
315 { 49.55809321e0, -10620.90088e0, -230.57416e0, } ,
316 { 100.46440702e0, 6362.03561e0, 326.52178e0, } ,
317 { 113.66550252e0, -9240.19942e0, -66.23743e0, } ,
318 { 74.00595701e0, 2669.15033e0, 145.93964e0, } ,
319 { 131.78405702e0, -221.94322e0, -0.78728e0, } ,
320
321 } ;
322
323 /*
324 * Tables for trigonometric terms to be added to the mean elements
325 * of the semi-major axes.
326 */
327
328 static const double KP[8][9] = {
329
330 { 69613.0, 75645.0, 88306.0, 59899.0, 15746.0, 71087.0, 142173.0, 3086.0, 0.0, } ,
331 { 21863.0, 32794.0, 26934.0, 10931.0, 26250.0, 43725.0, 53867.0, 28939.0, 0.0, } ,
332 { 16002.0, 21863.0, 32004.0, 10931.0, 14529.0, 16368.0, 15318.0, 32794.0, 0.0, } ,
333 { 6345.0, 7818.0, 15636.0, 7077.0, 8184.0, 14163.0, 1107.0, 4872.0, 0.0, } ,
334 { 1760.0, 1454.0, 1167.0, 880.0, 287.0, 2640.0, 19.0, 2047.0, 1454.0, } ,
335 { 574.0, 0.0, 880.0, 287.0, 19.0, 1760.0, 1167.0, 306.0, 574.0, } ,
336 { 204.0, 0.0, 177.0, 1265.0, 4.0, 385.0, 200.0, 208.0, 204.0, } ,
337 { 0.0, 102.0, 106.0, 4.0, 98.0, 1367.0, 487.0, 204.0, 0.0, } ,
338
339 } ;
340
341 static const double CA[8][9] = {
342
343 { 4.0, -13.0, 11.0, -9.0, -9.0, -3.0, -1.0, 4.0, 0.0, } ,
344 { -156.0, 59.0, -42.0, 6.0, 19.0, -20.0, -10.0, -12.0, 0.0, } ,
345 { 64.0, -152.0, 62.0, -8.0, 32.0, -41.0, 19.0, -11.0, 0.0, } ,
346 { 124.0, 621.0, -145.0, 208.0, 54.0, -57.0, 30.0, 15.0, 0.0, } ,
347 { -23437.0, -2634.0, 6601.0, 6259.0, -1507.0, -1821.0, 2620.0, -2115.0, -1489.0, } ,
348 { 62911.0, -119919.0, 79336.0, 17814.0, -24241.0, 12068.0, 8306.0, -4893.0, 8902.0, } ,
349 { 389061.0, -262125.0, -44088.0, 8387.0, -22976.0, -2093.0, -615.0, -9720.0, 6633.0, } ,
350 { -412235.0, -157046.0, -31430.0, 37817.0, -9740.0, -13.0, -7449.0, 9644.0, 0.0, } ,
351
352 } ;
353
354 static const double SA[8][9] = {
355
356 { -29.0, -1.0, 9.0, 6.0, -6.0, 5.0, 4.0, 0.0, 0.0, } ,
357 { -48.0, -125.0, -26.0, -37.0, 18.0, -13.0, -20.0, -2.0, 0.0, } ,
358 { -150.0, -46.0, 68.0, 54.0, 14.0, 24.0, -28.0, 22.0, 0.0, } ,
359 { -621.0, 532.0, -694.0, -20.0, 192.0, -94.0, 71.0, -73.0, 0.0, } ,
360 { -14614.0, -19828.0, -5869.0, 1881.0, -4372.0, -2255.0, 782.0, 930.0, 913.0, } ,
361 { 139737.0, 0.0, 24667.0, 51123.0, -5102.0, 7429.0, -4095.0, -1976.0, -9566.0, } ,
362 { -138081.0, 0.0, 37205.0, -49039.0, -41901.0, -33872.0, -27037.0, -12474.0, 18797.0, } ,
363 { 0.0, 28492.0, 133236.0, 69654.0, 52322.0, -49577.0, -26430.0, -3593.0, 0.0, } ,
364
365 } ;
366
367 /*
368 * Tables giving the trigonometric terms to be added to the mean
369 * elements of the mean longitudes.
370 */
371
372 static const double KQ[8][10] = {
373
374 { 3086.0, 15746.0, 69613.0, 59899.0, 75645.0, 88306.0, 12661.0, 2658.0, 0.0, 0.0, } ,
375 { 21863.0, 32794.0, 10931.0, 73.0, 4387.0, 26934.0, 1473.0, 2157.0, 0.0, 0.0, } ,
376 { 10.0, 16002.0, 21863.0, 10931.0, 1473.0, 32004.0, 4387.0, 73.0, 0.0, 0.0, } ,
377 { 10.0, 6345.0, 7818.0, 1107.0, 15636.0, 7077.0, 8184.0, 532.0, 10.0, 0.0, } ,
378 { 19.0, 1760.0, 1454.0, 287.0, 1167.0, 880.0, 574.0, 2640.0, 19.0, 1454.0, } ,
379 { 19.0, 574.0, 287.0, 306.0, 1760.0, 12.0, 31.0, 38.0, 19.0, 574.0, } ,
380 { 4.0, 204.0, 177.0, 8.0, 31.0, 200.0, 1265.0, 102.0, 4.0, 204.0, } ,
381 { 4.0, 102.0, 106.0, 8.0, 98.0, 1367.0, 487.0, 204.0, 4.0, 102.0, } ,
382
383 } ;
384
385 static const double CL[8][10] = {
386
387 { 21.0, -95.0, -157.0, 41.0, -5.0, 42.0, 23.0, 30.0, 0.0, 0.0, } ,
388 { -160.0, -313.0, -235.0, 60.0, -74.0, -76.0, -27.0, 34.0, 0.0, 0.0, } ,
389 { -325.0, -322.0, -79.0, 232.0, -52.0, 97.0, 55.0, -41.0, 0.0, 0.0, } ,
390 { 2268.0, -979.0, 802.0, 602.0, -668.0, -33.0, 345.0, 201.0, -55.0, 0.0, } ,
391 { 7610.0, -4997.0, -7689.0, -5841.0, -2617.0, 1115.0, -748.0, -607.0, 6074.0, 354.0, } ,
392 { -18549.0, 30125.0, 20012.0, -730.0, 824.0, 23.0, 1289.0, -352.0, -14767.0, -2062.0, } ,
393 { -135245.0, -14594.0, 4197.0, -4030.0, -5630.0, -2898.0, 2540.0, -306.0, 2939.0, 1986.0, } ,
394 { 89948.0, 2103.0, 8963.0, 2695.0, 3682.0, 1648.0, 866.0, -154.0, -1963.0, -283.0, } ,
395
396 } ;
397
398 static const double SL[8][10] = {
399
400 { -342.0, 136.0, -23.0, 62.0, 66.0, -52.0, -33.0, 17.0, 0.0, 0.0, } ,
401 { 524.0, -149.0, -35.0, 117.0, 151.0, 122.0, -71.0, -62.0, 0.0, 0.0, } ,
402 { -105.0, -137.0, 258.0, 35.0, -116.0, -88.0, -112.0, -80.0, 0.0, 0.0, } ,
403 { 854.0, -205.0, -936.0, -240.0, 140.0, -341.0, -97.0, -232.0, 536.0, 0.0, } ,
404 { -56980.0, 8016.0, 1012.0, 1448.0, -3024.0, -3710.0, 318.0, 503.0, 3767.0, 577.0, } ,
405 { 138606.0, -13478.0, -4964.0, 1441.0, -1319.0, -1482.0, 427.0, 1236.0, -9167.0, -1918.0, } ,
406 { 71234.0, -41116.0, 5334.0, -4935.0, -1848.0, 66.0, 434.0, -1748.0, 3780.0, -701.0, } ,
407 { -47645.0, 11647.0, 2166.0, 3194.0, 679.0, 0.0, -244.0, -419.0, -2531.0, 48.0, } ,
408
409 } ;
410
411 /*
412 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
413 */
414
415 /*
416 * Validate the planet number.
417 */
418
419 if ( np < 1 || np > 8 ) {
420
421 jstat = -1 ;
422
423 /*
424 * Reset the result in case of failure.
425 */
426
427 for ( k = 0; k < 2; k++ ) {
428 for ( i = 0; i < 3; i++ ) {
429 pv[k][i] = 0.0 ;
430 }
431 }
432
433 }
434 else {
435
436 /*
437 * Time: Julian millennia since J2000.
438 */
439
440 t = ( ( date1 - GAL_J2000 ) + date2 ) / GAL_DJM ;
441
442 /*
443 * OK status unless remote epoch.
444 */
445
446 if ( fabs ( t ) <= 1.0 ) {
447 jstat = 0 ;
448 }
449 else {
450 jstat = 1 ;
451 }
452
453 /*
454 * Adjust Planet number for C indicies starting at 0 not 1
455 */
456
457 np-- ;
458
459 /*
460 * Compute the mean elements.
461 */
462
463 da = A[np][0] +
464 ( A[np][1] +
465 A[np][2] * t ) * t ;
466
467 dl = ( 3600.0 * DLM[np][0] +
468 ( DLM[np][1] +
469 DLM[np][2] * t ) * t ) * GAL_AS2R ;
470
471 de = E[np][0] +
472 ( E[np][1] +
473 E[np][2] * t ) * t ;
474
475 dp = gal_anpm ( ( 3600.0 * PI[np][0] +
476 ( PI[np][1] +
477 PI[np][2] * t ) * t ) * GAL_AS2R ) ;
478
479 di = ( 3600.0 * DINC[np][0] +
480 ( DINC[np][1] +
481 DINC[np][2] * t ) * t ) * GAL_AS2R ;
482
483 dom = gal_anpm ( ( 3600.0 * OMEGA[np][0]
484 + ( OMEGA[np][1]
485 + OMEGA[np][2] * t ) * t ) * GAL_AS2R ) ;
486
487 /*
488 * Apply the trigonometric terms.
489 */
490
491 dmu = 0.35953620 * t ;
492
493 for (k=0; k<8; k++) {
494
495 arga = KP[np][k] * dmu ;
496
497 argl = KQ[np][k] * dmu ;
498
499 da += ( CA[np][k] * cos ( arga ) +
500 SA[np][k] * sin ( arga ) ) * 1e-7 ;
501
502 dl += ( CL[np][k] * cos ( argl ) +
503 SL[np][k] * sin ( argl ) ) * 1e-7 ;
504
505 }
506
507 arga = KP[np][8] * dmu ;
508
509 da += t * ( CA[np][8] * cos ( arga ) +
510 SA[np][8] * sin ( arga ) ) * 1e-7 ;
511
512 for (k=8; k<10; k++) {
513
514 argl = KQ[np][k] * dmu ;
515
516 dl += t * ( CL[np][k] * cos ( argl ) +
517 SL[np][k] * sin ( argl ) ) * 1e-7 ;
518
519 }
520
521 dl = fmod ( dl, GAL_2PI ) ;
522
523 /*
524 * Iterative solution of Kepler's equation to get eccentric anomaly.
525 */
526
527 am = dl - dp ;
528 ae = am + de * sin ( am ) ;
529
530 k = 0 ;
531 do {
532 dae = ( am - ae + de * sin ( ae ) ) / ( 1.0 - de * cos ( ae ) ) ;
533 ae += dae ;
534 k++ ;
535 if ( k >= KMAX ) {
536 jstat = 2 ;
537 }
538 } while ( k < KMAX && fabs ( dae ) > 1e-12 ) ;
539
540 /*
541 * True anomaly.
542 */
543
544 ae2 = ae / 2.0 ;
545 at = 2.0 * atan2 ( sqrt ( ( 1.0 + de ) / ( 1.0 - de ) ) * sin ( ae2 ), cos ( ae2 ) ) ;
546
547 /*
548 * Distance (AU) and speed (radians per day).
549 */
550
551 r = da * ( 1.0 - de * cos ( ae ) ) ;
552 v = GK * sqrt ( ( 1.0 + 1.0 / AMAS[np] ) / ( da * da * da ) ) ;
553
554 si2 = sin ( di / 2.0 ) ;
555 xq = si2 * cos ( dom ) ;
556 xp = si2 * sin ( dom ) ;
557 tl = at + dp ;
558 xsw = sin ( tl ) ;
559 xcw = cos ( tl ) ;
560 xm2 = 2.0 * ( xp * xcw - xq * xsw ) ;
561 xf = da / sqrt ( 1.0 - de * de ) ;
562 ci2 = cos( di / 2.0 ) ;
563 xms = ( de * sin ( dp ) + xsw ) * xf ;
564 xmc = ( de * cos ( dp ) + xcw ) * xf ;
565 xpxq2 = 2.0 * xp * xq ;
566
567 /*
568 * Position (J2000 ecliptic x,y,z in AU).
569 */
570
571 x = r * ( xcw - xm2 * xp ) ;
572 y = r * ( xsw + xm2 * xq ) ;
573 z = r * ( -xm2 * ci2 ) ;
574
575 /*
576 * Rotate to equatorial.
577 */
578
579 pv[0][0] = x ;
580 pv[0][1] = y * COSEPS - z * SINEPS ;
581 pv[0][2] = y * SINEPS + z * COSEPS ;
582
583 /*
584 * Velocity (J2000 ecliptic xdot,ydot,zdot in AU/d).
585 */
586
587 x = v * ( ( -1.0 + 2.0 * xp * xp ) * xms + xpxq2 * xmc ) ;
588 y = v * ( ( 1.0 - 2.0 * xq * xq ) * xmc - xpxq2 * xms ) ;
589 z = v * ( 2.0 * ci2 * ( xp * xms + xq * xmc ) ) ;
590
591 /*
592 * Rotate to equatorial.
593 */
594
595 pv[1][0] = x ;
596 pv[1][1] = y * COSEPS - z * SINEPS ;
597 pv[1][2] = y * SINEPS + z * COSEPS ;
598
599 }
600
601 /*
602 * Return the status.
603 */
604
605 return jstat ;
606
607 /*
608 * Finished.
609 */
610
611 }
612
613 /*
614 * gal - General Astrodynamics Library
615 * Copyright (C) 2008 Paul C. L. Willmott
616 *
617 * This program is free software; you can redistribute it and/or modify
618 * it under the terms of the GNU General Public License as published by
619 * the Free Software Foundation; either version 2 of the License, or
620 * (at your option) any later version.
621 *
622 * This program is distributed in the hope that it will be useful,
623 * but WITHOUT ANY WARRANTY; without even the implied warranty of
624 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
625 * GNU General Public License for more details.
626 *
627 * You should have received a copy of the GNU General Public License along
628 * with this program; if not, write to the Free Software Foundation, Inc.,
629 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
630 *
631 * Contact:
632 *
633 * Paul Willmott
634 * vp9mu@amsat.org
635 */
636