1 /*
2 * - - - - - - - - - - - -
3 * g a l _ h p l p v 8 7
4 * - - - - - - - - - - - -
5 *
6 * This routine is part of the General Astrodynamics Library
7 *
8 * Description:
9 *
10 * Pluto heliocentric position and velocity, with respect to the
11 * FK5 Reference Frame.
12 *
13 * Status:
14 *
15 * internal support routine.
16 *
17 * Given:
18 *
19 * tdb1 d TDB epoch part 1 (Note 1)
20 * tdb2 d TDB epoch part 2 (Note 1)
21 * ref i Reference frame
22 * 0 = dynamical equinox and ecliptic J2000.
23 * 1 = FK5 (VSOP87)
24 *
25 * Returned:
26 *
27 * pv d[2][3] position/velocity (AU, AU/Day)
28 *
29 * Called:
30 *
31 * gal_zpv Zero pv-vector
32 * gal_sxpv Scale pv-vector
33 * gal_rxpv Product of r-matrix and pv-vector
34 *
35 * Notes:
36 *
37 * 1) The epoch tt1+tt2 is a Julian Date, apportioned in
38 * any convenient way between the two arguments. For example,
39 * JD(TDB)=2450123.7 could be expressed in any of these ways,
40 * among others:
41 *
42 * tdb1 tdb2
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 * However, the accuracy of the result is more likely to be
56 * limited by the algorithm itself than the way the epoch has been
57 * expressed.
58 *
59 * 2) On return, the arrays pvh and pvb contain the following:
60 *
61 * pv[0][0] x }
62 * pv[0][1] y } position, AU
63 * pv[0][2] z }
64 *
65 * pv[1][0] xdot }
66 * pv[1][1] ydot } velocity, AU/day
67 * pv[1][2] zdot }
68 *
69 * The vectors are heliocentric with respect to the FK5 Reference Frame.
70 * The time unit is one day in TT
71 *
72 * 3) Computation of heliocentric rectangular coordinates of Pluto.
73 * Source : DE200 (Jet Propulsion Laboratory).
74 * Frame : dynamical mean equinox and equator J2000 (DE200).
75 * The ephemerides are computed with a method of approximation
76 * using frequency analysis.
77 *
78 * 4) Valid date range is JD2341972.5 (Jan 01 1700 0h) to
79 * JD2488092.5 (Jan 01 2100 0h).
80 *
81 * References:
82 *
83 * J. Chapront, 1995, Astron. Atrophys. Sup. Ser., 109, 181.
84 *
85 * This revision:
86 *
87 * 2008 July 5
88 *
89 * Copyright (C) 2008 Paul C. L. Willmott. See notes at end.
90 *
91 *-----------------------------------------------------------------------
92 */
93
94 #include <math.h>
95 #include "gal_hplpv87.h"
96 #include "gal_zpv.h"
97 #include "gal_sxpv.h"
98 #include "gal_rxpv.h"
99
100 void
gal_hplpv87(double tdb1,double tdb2,int ref,double pv[2][3])101 gal_hplpv87
102 (
103 double tdb1,
104 double tdb2,
105 int ref,
106 double pv[2][3]
107 )
108
109 {
110
111 double v[2][3], w[3] ;
112
113 /*
114 * Initializations
115 */
116
117 static const int nf[3] = { 82, 19, 5 } ;
118
119 static const double tdeb = 2341972.5 ;
120 static const double dt = 146120.0 ;
121
122 /*
123 * Secular Series
124 */
125
126 static const double ax[4] = { 98083308510.0, -1465718392.0, 11528487809.0, 55397965917.0 } ;
127 static const double ay[4] = { 101846243715.0, 57789.0, -5487929294.0, 8520205290.0 } ;
128 static const double az[4] = { 2183700004.0, 433209785.0, -4911803413.0, -14029741184.0 } ;
129
130 /*
131 * Frequencies
132 */
133
134 static double fq[106] = {
135 0.0000645003954767,0.0001083248054773,0.0001302772403167,
136 0.0001647868659960,0.0001935009111902,0.0002223740247147,
137 0.0003032575201026,0.0003259246239385,0.0003564763034914,
138 0.0004265811293132,0.0004503959517513,0.0004638675148284,
139 0.0005009272733421,0.0005163593863414,0.0005578826828210,
140 0.0005882795362847,0.0006450023602974,0.0007097635821639,
141 0.0007630643253588,0.0007740033551209,0.0008385031396726,
142 0.0008950591609720,0.0009545118163938,0.0010255417569600,
143 0.0010826728325744,0.0011680358909203,0.0012405125052369,
144 0.0012931805883876,0.0013460706181008,0.0014190059530383,
145 0.0014394705053002,0.0014502634075377,0.0014992014575181,
146 0.0015434430430867,0.0016000710611098,0.0016562809940875,
147 0.0017275924266291,0.0017454042542465,0.0018215079641428,
148 0.0018694826929211,0.0019274630193251,0.0020276790928706,
149 0.0021822818660433,0.0022885289854970,0.0023167646379420,
150 0.0023445464874575,0.0024069306189938,0.0024473146628449,
151 0.0024778027974419,0.0025244208011161,0.0025682157855485,
152 0.0026028617439482,0.0026544444009919,0.0026987455959123,
153 0.0027308225916697,0.0027735113723168,0.0028728385464030,
154 0.0029001725379479,0.0029379670182566,0.0029750359447782,
155 0.0031326820696785,0.0031822107498712,0.0031931048857857,
156 0.0032268922327691,0.0034657232066225,0.0037838581645670,
157 0.0038055149432355,0.0038631344783149,0.0039129259467328,
158 0.0040311445462510,0.0040607542008930,0.0041490103414206,
159 0.0043500678052272,0.0046321937641054,0.0058000376725240,
160 0.0091460971544658,0.0091560629947357,0.0172021239411871,
161 0.0182919855069063,0.0279624510118796,0.0344040996177640,
162 0.0714245719830324,0.0001083248054773,0.0001647868659960,
163 0.0003032575201026,0.0003259246239385,0.0003564763034914,
164 0.0005009272733421,0.0005882795362847,0.0007097635821639,
165 0.0007630643253588,0.0008950591609720,0.0011680358909203,
166 0.0014502634075377,0.0017454042542465,0.0029001725379479,
167 0.0031822107498712,0.0040607542008930,0.0043500678052272,
168 0.0091460971544658,0.0279624510118796,0.0001083248054773,
169 0.0003032575201026,0.0011680358909203,0.0043500678052272,
170 0.0279624510118796,
171 } ;
172
173 /*
174 * Coefficients X (cosine)
175 */
176
177 static const double cx[106] = {
178 -16338582222.0, -5995086437.0, 23663880362.0, 10304632056.0,
179 -3996936944.0, -4136465568.0, 1188702881.0, -621434363.0,
180 566898160.0, -75880391.0, 576146406.0, -659684298.0,
181 451962774.0, -153724334.0, -603163280.0, 364764379.0,
182 193062130.0, 161493959.0, 1167349082.0, -1417467887.0,
183 15325240.0, -3624391.0, -587306.0, 132022.0,
184 -106501.0, 228373.0, -95106.0, 56299.0,
185 -48339.0, 803937.0, -6172744.0, -18962749.0,
186 133022.0, -25964.0, 7111.0, -4998.0,
187 32034.0, -29666.0, -1983.0, 114.0,
188 191.0, -1063.0, 419.0, 346.0,
189 5059.0, -81.0, 1408.0, 2964.0,
190 -5364.0, 1509.0, -4924.0, 2954.0,
191 2034.0, -5199.0, 604.0, -1247.0,
192 4576.0, -350741.0, -4023.0, 1147.0,
193 -38.0, -99.0, -11686.0, 1129.0,
194 582.0, -83.0, -97.0, 431.0,
195 -134.0, -323.0, -292.0, 195.0,
196 39068.0, 523.0, -1747.0, 3135.0,
197 -619.0, -12095.0, 6.0, 18476.0,
198 -130.0, -438.0, 102345278799.0, -9329130892.0,
199 1484339404.0, 472660593.0, -581239444.0, 1016663241.0,
200 -1054199614.0, 99039105.0, -52190030.0, -3394173.0,
201 -16529.0, 3102430.0, 2286.0, -10955.0,
202 -5293.0, -654.0, 124.0, -85.0,
203 29.0, 418209651.0, -1191875710.0, -823081.0,
204 -558.0, -1091.0,
205 } ;
206
207 /*
208 * Coefficients X (sine)
209 */
210
211 static const double sx[106] = {
212 -308294137468.0,-68820910480.0, 28346466257.0, -1755658975.0,
213 7818660837.0, -1098895702.0, -1192462299.0, -772129982.0,
214 1061702581.0, -639572722.0, 1128327488.0, -423570428.0,
215 -175317704.0, 251601606.0, -869448807.0, 551228298.0,
216 87807522.0, -11540541.0, -103236703.0, 92638954.0,
217 -3624991.0, 1004975.0, 304396.0, -56532.0,
218 55554.0, -799096.0, 56947.0, -48016.0,
219 50599.0, -680660.0, 5858452.0, 38125648.0,
220 -109460.0, 18684.0, -5269.0, 2771.0,
221 -6814.0, 47130.0, 1192.0, -1387.0,
222 379.0, -612.0, -52.0, 813.0,
223 -4354.0, -2275.0, 685.0, -1352.0,
224 4681.0, -1908.0, -6530.0, 8667.0,
225 1675.0, 874.0, 898.0, 965.0,
226 -7124.0, -1145389.0, 2931.0, -618.0,
227 -34.0, -6562.0, 8038.0, -697.0,
228 -8.0, 12.0, -267.0, -131.0,
229 304.0, -756.0, -103.0, -250.0,
230 19816.0, -596.0, 576.0, 4122.0,
231 65.0, -27900.0, 217.0, -137.0,
232 -269.0, 531.0, -24338350765.0, 11210995713.0,
233 2793567155.0, -776019789.0, 1528323591.0, -249354416.0,
234 1127608109.0, -667692329.0, -1570766679.0, -9724425.0,
235 26552.0, 3332520.0, -27607.0, -11696.0,
236 -7297.0, -104.0, -184.0, -455.0,
237 -16.0, 39813894679.0, 3633087275.0, 522728.0,
238 -320.0, -1401.0,
239 } ;
240
241 /*
242 * Coefficients Y (cosine)
243 */
244 static const double cy[106] = {
245 299584895562.0, 75951634908.0, -36135662843.0, 18125610071.0,
246 -20398008415.0, 6125780503.0, -162559485.0, 4352425804.0,
247 -3819676998.0, 1168107376.0, -5041323701.0, 4093828501.0,
248 -1727274544.0, 134214260.0, 5033950069.0, -3071449401.0,
249 -1190419055.0, -775881742.0, -5524713888.0, 6803228005.0,
250 -65675611.0, 15155413.0, 2009509.0, -389682.0,
251 275571.0, 474366.0, 132163.0, -81550.0,
252 69996.0, -706470.0, 4777898.0, -44002785.0,
253 -58735.0, 7624.0, -1922.0, -729.0,
254 -1733.0, -35642.0, -586.0, -258.0,
255 -368.0, 1286.0, -136.0, 883.0,
256 2673.0, 331.0, 50.0, 178.0,
257 2901.0, -654.0, -8972.0, 3034.0,
258 1113.0, 570.0, -72.0, 1950.0,
259 8550.0, 1047593.0, -2348.0, 313.0,
260 432.0, 6765.0, -8240.0, 335.0,
261 140.0, -833.0, 252.0, -210.0,
262 366.0, -920.0, 1215.0, -217.0,
263 -17780.0, 581.0, -560.0, -4131.0,
264 390.0, 25613.0, -206.0, 1850.0,
265 171.0, -471.0, 26437625772.0,-12674907683.0,
266 -1067899665.0, -2082744.0, -43195632.0, 211912497.0,
267 -108307161.0, -63033809.0, -203850703.0, -1672332.0,
268 7136.0, 803655.0, -10985.0, 9126.0,
269 3317.0, -151.0, 160.0, 138.0,
270 -27.0,-36463065062.0, -5816560445.0, 1576292.0,
271 -21.0, -295.0,
272 } ;
273
274 /*
275 * Coefficients Y (sine)
276 */
277
278 static const double sy[106] = {
279 -53545027809.0, -8838029861.0, 23553788174.0, 13775798112.0,
280 -6068121593.0, -2853107588.0, 750355551.0, -82067770.0,
281 230091832.0, -259838942.0, 197944074.0, 27141006.0,
282 -105334544.0, 95175918.0, -139461973.0, 80593104.0,
283 -5126842.0, -21953793.0, -163767784.0, 192436228.0,
284 -2479113.0, 561687.0, 121909.0, -30275.0,
285 16333.0, 68105.0, 24081.0, -11228.0,
286 667.0, -73047.0, 1007089.0, -22814549.0,
287 434.0, 1013.0, 710.0, 1100.0,
288 -4598.0, 1990.0, 564.0, 828.0,
289 -1119.0, -1249.0, -597.0, 227.0,
290 5467.0, 801.0, -2029.0, -1892.0,
291 4713.0, -459.0, 1757.0, -9303.0,
292 -2357.0, 7679.0, -2953.0, 629.0,
293 5011.0, -333905.0, -2388.0, 415.0,
294 139.0, -5726.0, -4583.0, 310.0,
295 681.0, -107.0, 301.0, -525.0,
296 198.0, -379.0, -230.0, -64.0,
297 36069.0, 459.0, -1596.0, 2509.0,
298 -146.0, -11081.0, 4.0, 15764.0,
299 -147.0, -362.0, 117449924600.0, -7691661502.0,
300 -4771148239.0, 3733883366.0, -7081845126.0, 3502526523.0,
301 -8115570206.0, 3607883959.0, 7690328772.0, 37384011.0,
302 -164319.0, -2859257.0, 1593.0, -11997.0,
303 -6476.0, 1419.0, 34.0, 232.0,
304 32.0, 2752753498.0, -672124207.0, 154239.0,
305 -400.0, 372.0,
306 } ;
307
308 /*
309 * Coefficients Z (cosine)
310 */
311
312 static const double cz[106] = {
313 98425296138.0, 25475793908.0, -18424386574.0, 2645968636.0,
314 -5282207967.0, 3278235471.0, -425422632.0, 1526641086.0,
315 -1323182752.0, 235873266.0, -1617466723.0, 1557465867.0,
316 -848586296.0, 218182986.0, 1636044515.0, -1001334243.0,
317 -455739370.0, -348173978.0, -2511254281.0, 3062521470.0,
318 -32079379.0, 7597939.0, 1138566.0, -238849.0,
319 192377.0, 83169.0, 148694.0, -92489.0,
320 87116.0, -1281070.0, 9950106.0, -25105642.0,
321 -171749.0, 31035.0, -8648.0, 5360.0,
322 -30345.0, 11482.0, 1322.0, -467.0,
323 96.0, 894.0, -381.0, -583.0,
324 2525.0, -569.0, 226.0, -2039.0,
325 3728.0, -1540.0, 42.0, -3144.0,
326 658.0, 220.0, 1848.0, 678.0,
327 -7289.0, 463291.0, 3945.0, -1141.0,
328 -26.0, -10607.0, 11458.0, -1005.0,
329 120.0, -301.0, 135.0, -186.0,
330 118.0, 30.0, 197.0, -182.0,
331 -8585.0, 240.0, -226.0, -2049.0,
332 283.0, 11109.0, -100.0, -842.0,
333 71.0, -181.0, -22591501373.0, -1138977908.0,
334 -782718600.0, -141483824.0, 159033355.0, -246222739.0,
335 287284767.0, -48002332.0, -41114335.0, 578004.0,
336 -8420.0, -766779.0, 957.0, 5780.0,
337 4141.0, 417.0, -8.0, 65.0,
338 -22.0,-11656050047.0, -1186276469.0, 1388681.0,
339 201.0, 561.0,
340 } ;
341
342 /*
343 * Coefficients Z (sine)
344 */
345
346 static const double sz[106] = {
347 76159403805.0, 17987340882.0, -1193982379.0, 4828308190.0,
348 -4248985438.0, -559147671.0, 593594960.0, 208799497.0,
349 -249913200.0, 115051024.0, -282588988.0, 135883560.0,
350 23091693.0, -49187976.0, 223956575.0, -137344299.0,
351 -28188872.0, -2636274.0, -14202661.0, 25488216.0,
352 419837.0, -150966.0, -64906.0, 3719.0,
353 -2226.0, 86321.0, -15970.0, 16609.0,
354 -15782.0, 200300.0, -1500491.0, -9161491.0,
355 37481.0, -4616.0, 224.0, -1027.0,
356 5220.0, -6976.0, -267.0, 556.0,
357 -23.0, -711.0, -122.0, -97.0,
358 2440.0, 786.0, -806.0, -167.0,
359 -156.0, 572.0, 2532.0, -4582.0,
360 -1178.0, 875.0, -558.0, 781.0,
361 3230.0, -116132.0, -1440.0, 438.0,
362 176.0, 1072.0, -5850.0, 418.0,
363 267.0, 60.0, 134.0, -85.0,
364 -59.0, 112.0, -168.0, -89.0,
365 14986.0, 190.0, -685.0, 1018.0,
366 -48.0, -4807.0, 0.0, 7066.0,
367 -54.0, -229.0, 44126663549.0, -5626220823.0,
368 -2536450838.0, 1536292657.0, -2916144530.0, 949074586.0,
369 -2842935040.0, 1500396857.0, 3415136438.0, 19702076.0,
370 -46995.0, -5801645.0, 33470.0, 17674.0,
371 7355.0, 199.0, 11.0, 205.0,
372 33.0,-11127973411.0, -1310869292.0, -164753.0,
373 -107.0, 284.0,
374 } ;
375
376 double fk5[2][3][3] = {
377 {
378 { +1.000000000000, +0.000000000000, +0.000000000000, } ,
379 { +0.000000000000, +1.000000000000, +0.000000000000, } ,
380 { 0.000000000000, +0.000000000000, +1.000000000000, } ,
381 } ,
382 {
383 { +1.000000000000, +0.000000440360, -0.000000190919, } ,
384 { -0.000000479966, +0.917482137087, -0.397776982902, } ,
385 { 0.000000000000, +0.397776982902, +0.917482137087, } ,
386 } ,
387 } ;
388
389 double x, fx, wx, f, sf, cf, fw, wt, wy, pvt[2][3] ;
390 int i, iv, imax, imin, m ;
391
392 /*
393 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
394 */
395
396 /*
397 * Results in pv
398 */
399
400 gal_zpv ( pv ) ;
401
402 /*
403 * Change variable
404 */
405
406 x = 2.0 * ( tdb1 - tdeb + tdb2 ) / dt - 1.0 ;
407 fx = x * dt / 2.0 ;
408
409 /*
410 * Compute the positions (secular terms)
411 */
412
413 gal_zpv ( v ) ;
414 wx = 1.0 ;
415 for ( i = 0; i <=3; i++ ) {
416 v[0][0] += ax[i] * wx ;
417 v[0][1] += ay[i] * wx ;
418 v[0][2] += az[i] * wx ;
419 wx *= x ;
420 }
421
422 /*
423 * Compute the positions (periodic and Poisson terms)
424 */
425
426 imax = 0 ;
427 wx = 1.0 ;
428 for ( m = 0; m <=2; m++ ) {
429 gal_zp ( w ) ;
430 imin = imax + 1 ;
431 imax = imax + nf[m] ;
432 for ( i = imin - 1; i < imax; i++ ) {
433 f = fq[i] * fx ;
434 cf = cos ( f ) ;
435 sf = sin ( f ) ;
436 w[0] += cx[i] * cf + sx[i] * sf ;
437 w[1] += cy[i] * cf + sy[i] * sf ;
438 w[2] += cz[i] * cf + sz[i] * sf ;
439 }
440 for ( iv = 0; iv < 3; iv++ ) {
441 v[0][iv] += w[iv] * wx ;
442 }
443 wx *= x ;
444 }
445
446 /*
447 * Compute the velocities (secular terms)
448 */
449
450 wt = 2.0 / dt ;
451 wx = 1.0 ;
452 for ( i = 1; i <= 3; i++ ) {
453 v[1][0] += i * ax[i] * wx ;
454 v[1][1] += i * ay[i] * wx ;
455 v[1][2] += i * az[i] * wx ;
456 wx *= x ;
457 }
458 v[1][0] *= wt ;
459 v[1][1] *= wt ;
460 v[1][2] *= wt ;
461
462 /*
463 * Compute the velocities (periodic and Poisson terms)
464 */
465
466 imax = 0 ;
467 wx = 1 ;
468 for ( m = 0; m <= 2; m++ ) {
469 imin = imax + 1 ;
470 imax = imax + nf[m] ;
471 for ( i = imin - 1; i < imax; i++ ) {
472 fw = fq[i] ;
473 f = fw * fx ;
474 cf = cos ( f ) ;
475 sf = sin ( f ) ;
476 v[1][0] += fw * ( sx[i] * cf - cx[i] * sf ) * wx ;
477 v[1][1] += fw * ( sy[i] * cf - cy[i] * sf ) * wx ;
478 v[1][2] += fw * ( sz[i] * cf - cz[i] * sf ) * wx ;
479 if ( m > 0 ) {
480 v[1][0] += m * wt * ( cx[i] * cf + sx[i] * sf ) * wy ;
481 v[1][1] += m * wt * ( cy[i] * cf + sy[i] * sf ) * wy ;
482 v[1][2] += m * wt * ( cz[i] * cf + sz[i] * sf ) * wy ;
483 }
484 }
485 wy = wx ;
486 wx *= x ;
487 }
488
489 /*
490 * Stock the results
491 */
492
493 gal_sxpv ( 1.0e-10, v, pvt ) ;
494
495 /*
496 * Convert to FK5 equitorial
497 */
498
499 gal_rxpv ( fk5[ref], pvt, pv ) ;
500
501 /*
502 * Finished.
503 */
504
505 }
506
507 /*
508 * gal - General Astrodynamics Library
509 * Copyright (C) 2008 Paul C. L. Willmott
510 *
511 * This program is free software; you can redistribute it and/or modify
512 * it under the terms of the GNU General Public License as published by
513 * the Free Software Foundation; either version 2 of the License, or
514 * (at your option) any later version.
515 *
516 * This program is distributed in the hope that it will be useful,
517 * but WITHOUT ANY WARRANTY; without even the implied warranty of
518 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
519 * GNU General Public License for more details.
520 *
521 * You should have received a copy of the GNU General Public License along
522 * with this program; if not, write to the Free Software Foundation, Inc.,
523 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
524 *
525 * Contact:
526 *
527 * Paul Willmott
528 * vp9mu@amsat.org
529 */
530