source: trunk/MagicSoft/slalib/dmoon.c

Last change on this file was 731, checked in by tbretz, 24 years ago
*** empty log message ***
  • Property svn:executable set to *
File size: 20.6 KB
Line 
1#include "slalib.h"
2#include "slamac.h"
3void slaDmoon ( double date, double pv[6] )
4/*
5** - - - - - - - - -
6** s l a D m o o n
7** - - - - - - - - -
8**
9** Approximate geocentric position and velocity of the Moon
10** (double precision).
11**
12** Given:
13** date double TDB (loosely ET) as a Modified Julian Date
14** (JD-2400000.5)
15**
16** Returned:
17** pv double[6] Moon x,y,z,xdot,ydot,zdot, mean equator
18** and equinox of date (AU, AU/s)
19**
20** Notes:
21**
22** 1 This routine is a full implementation of the algorithm
23** published by Meeus (see reference).
24**
25** 2 Meeus quotes accuracies of 10 arcsec in longitude, 3 arcsec in
26** latitude and 0.2 arcsec in HP (equivalent to about 20 km in
27** distance). Comparison with JPL DE200 over the interval
28** 1960-2025 gives RMS errors of 3.7 arcsec and 83 mas/hour in
29** longitude, 2.3 arcsec and 48 mas/hour in latitude, 11 km
30** and 81 mm/s in distance. The maximum errors over the same
31** interval are 18 arcsec and 0.50 arcsec/hour in longitude,
32** 11 arcsec and 0.24 arcsec/hour in latitude, 40 km and 0.29 m/s
33** in distance.
34**
35** 3 The original algorithm is expressed in terms of the obsolete
36** timescale Ephemeris Time. Either TDB or TT can be used, but
37** not UT without incurring significant errors (30 arcsec at
38** the present time) due to the Moon's 0.5 arcsec/sec movement.
39**
40** 4 The algorithm is based on pre IAU 1976 standards. However,
41** the result has been moved onto the new (FK5) equinox, an
42** adjustment which is in any case much smaller than the
43** intrinsic accuracy of the procedure.
44**
45** 5 Velocity is obtained by a complete analytical differentiation
46** of the Meeus model.
47**
48** Reference:
49** Meeus, l'Astronomie, June 1984, p348.
50**
51** Defined in slamac.h: DD2R, DAS2R, DS2R, dmod
52**
53** Last revision: 22 January 1998
54**
55** Copyright P.T.Wallace. All rights reserved.
56*/
57
58#define CJ 3155760000.0 /* Seconds per Julian century */
59 /* ( = 86400 * 36525 ) */
60
61#define ERADAU 4.2635212653763e-5 /* Earth equatorial radius in AU */
62 /* ( = 6378.137 / 149597870 ) */
63
64#define B1950 1949.9997904423 /* Julian epoch of B1950 */
65
66{
67 double t, theta, sinom, cosom, domcom, wa, dwa, wb, dwb, wom,
68 dwom, sinwom, coswom, v, dv, coeff, emn, empn, dn, fn, en,
69 den, dtheta, ftheta, el, del, b, db, bf, dbf, p, dp, sp, r,
70 dr, x, y, z, xd, yd, zd, sel, cel, sb, cb, rcb, rbd, w, epj,
71 eqcor, eps, sineps, coseps, es, ec;
72 int n, i;
73
74/*
75** Coefficients for fundamental arguments
76**
77** at J1900: T^0, T^1, T^2, T^3
78** at epoch: T^0, T^1
79**
80** Units are degrees for position and Julian centuries for time.
81*/
82
83/* Moon's mean longitude */
84 static double elp0 = 270.434164,
85 elp1 = 481267.8831,
86 elp2 = -0.001133,
87 elp3 = 0.0000019;
88 double elp, delp;
89
90/* Sun's mean anomaly */
91 static double em0 = 358.475833,
92 em1 = 35999.0498,
93 em2 = -0.000150,
94 em3 = -0.0000033;
95 double em, dem;
96
97/* Moon's mean anomaly */
98 static double emp0 = 296.104608,
99 emp1 = 477198.8491,
100 emp2 = 0.009192,
101 emp3 = 0.0000144;
102 double emp, demp;
103
104/* Moon's mean elongation */
105 static double d0 = 350.737486,
106 d1 = 445267.1142,
107 d2 = -0.001436,
108 d3 = 0.0000019;
109 double d, dd;
110
111/* Mean distance of the Moon from its ascending node */
112 static double f0 = 11.250889,
113 f1 = 483202.0251,
114 f2 = -0.003211,
115 f3 = -0.0000003;
116 double f, df;
117
118/* Longitude of the Moon's ascending node */
119 static double om0 = 259.183275,
120 om1 = -1934.1420,
121 om2 = 0.002078,
122 om3 = 0.0000022;
123 double om, dom;
124
125/* Coefficients for (dimensionless) E factor */
126 static double e1 = -0.002495,
127 e2 = -0.00000752;
128 double e, de, esq, desq;
129
130/* Coefficients for periodic variations etc */
131 static double pac = 0.000233, pa0 = 51.2,
132 pa1 = 20.2;
133 static double pbc = -0.001778;
134 static double pcc = 0.000817;
135 static double pdc = 0.002011;
136 static double pec = 0.003964, pe0 = 346.560,
137 pe1 = 132.870,
138 pe2 = -0.0091731;
139 static double pfc = 0.001964;
140 static double pgc = 0.002541;
141 static double phc = 0.001964;
142 static double pic = -0.024691;
143 static double pjc = -0.004328, pj0 = 275.05,
144 pj1 = -2.30;
145 static double cw1 = 0.0004664;
146 static double cw2 = 0.0000754;
147
148/*
149** Coefficients for Moon longitude, latitude, parallax series
150*/
151 struct term {
152 double coef; /* coefficient of L, B or P term (deg) */
153 int nem; /* multiple of M in argument */
154 int nemp; /* " " M' " " */
155 int nd; /* " " D " " */
156 int nf; /* " " F " " */
157 int ne; /* power of e to multiply term by */
158 };
159
160/*
161** Longitude coeff M M' D F n
162*/
163 static struct term tl[] = { { 6.288750, 0, 1, 0, 0, 0 },
164 { 1.274018, 0, -1, 2, 0, 0 },
165 { 0.658309, 0, 0, 2, 0, 0 },
166 { 0.213616, 0, 2, 0, 0, 0 },
167 { -0.185596, 1, 0, 0, 0, 1 },
168 { -0.114336, 0, 0, 0, 2, 0 },
169 { 0.058793, 0, -2, 2, 0, 0 },
170 { 0.057212, -1, -1, 2, 0, 1 },
171 { 0.053320, 0, 1, 2, 0, 0 },
172 { 0.045874, -1, 0, 2, 0, 1 },
173 { 0.041024, -1, 1, 0, 0, 1 },
174 { -0.034718, 0, 0, 1, 0, 0 },
175 { -0.030465, 1, 1, 0, 0, 1 },
176 { 0.015326, 0, 0, 2, -2, 0 },
177 { -0.012528, 0, 1, 0, 2, 0 },
178 { -0.010980, 0, -1, 0, 2, 0 },
179 { 0.010674, 0, -1, 4, 0, 0 },
180 { 0.010034, 0, 3, 0, 0, 0 },
181 { 0.008548, 0, -2, 4, 0, 0 },
182 { -0.007910, 1, -1, 2, 0, 1 },
183 { -0.006783, 1, 0, 2, 0, 1 },
184 { 0.005162, 0, 1, -1, 0, 0 },
185 { 0.005000, 1, 0, 1, 0, 1 },
186 { 0.004049, -1, 1, 2, 0, 1 },
187 { 0.003996, 0, 2, 2, 0, 0 },
188 { 0.003862, 0, 0, 4, 0, 0 },
189 { 0.003665, 0, -3, 2, 0, 0 },
190 { 0.002695, -1, 2, 0, 0, 1 },
191 { 0.002602, 0, 1, -2, -2, 0 },
192 { 0.002396, -1, -2, 2, 0, 1 },
193 { -0.002349, 0, 1, 1, 0, 0 },
194 { 0.002249, -2, 0, 2, 0, 2 },
195 { -0.002125, 1, 2, 0, 0, 1 },
196 { -0.002079, 2, 0, 0, 0, 2 },
197 { 0.002059, -2, -1, 2, 0, 2 },
198 { -0.001773, 0, 1, 2, -2, 0 },
199 { -0.001595, 0, 0, 2, 2, 0 },
200 { 0.001220, -1, -1, 4, 0, 1 },
201 { -0.001110, 0, 2, 0, 2, 0 },
202 { 0.000892, 0, 1, -3, 0, 0 },
203 { -0.000811, 1, 1, 2, 0, 1 },
204 { 0.000761, -1, -2, 4, 0, 1 },
205 { 0.000717, -2, 1, 0, 0, 2 },
206 { 0.000704, -2, 1, -2, 0, 2 },
207 { 0.000693, 1, -2, 2, 0, 1 },
208 { 0.000598, -1, 0, 2, -2, 1 },
209 { 0.000550, 0, 1, 4, 0, 0 },
210 { 0.000538, 0, 4, 0, 0, 0 },
211 { 0.000521, -1, 0, 4, 0, 1 },
212 { 0.000486, 0, 2, -1, 0, 0 } };
213 static int NL = ( sizeof tl / sizeof ( struct term ) );
214
215/*
216** Latitude coeff M M' D F n
217*/
218 static struct term tb[] = { { 5.128189, 0, 0, 0, 1, 0 },
219 { 0.280606, 0, 1, 0, 1, 0 },
220 { 0.277693, 0, 1, 0, -1, 0 },
221 { 0.173238, 0, 0, 2, -1, 0 },
222 { 0.055413, 0, -1, 2, 1, 0 },
223 { 0.046272, 0, -1, 2, -1, 0 },
224 { 0.032573, 0, 0, 2, 1, 0 },
225 { 0.017198, 0, 2, 0, 1, 0 },
226 { 0.009267, 0, 1, 2, -1, 0 },
227 { 0.008823, 0, 2, 0, -1, 0 },
228 { 0.008247, -1, 0, 2, -1, 1 },
229 { 0.004323, 0, -2, 2, -1, 0 },
230 { 0.004200, 0, 1, 2, 1, 0 },
231 { 0.003372, -1, 0, -2, 1, 1 },
232 { 0.002472, -1, -1, 2, 1, 1 },
233 { 0.002222, -1, 0, 2, 1, 1 },
234 { 0.002072, -1, -1, 2, -1, 1 },
235 { 0.001877, -1, 1, 0, 1, 1 },
236 { 0.001828, 0, -1, 4, -1, 0 },
237 { -0.001803, 1, 0, 0, 1, 1 },
238 { -0.001750, 0, 0, 0, 3, 0 },
239 { 0.001570, -1, 1, 0, -1, 1 },
240 { -0.001487, 0, 0, 1, 1, 0 },
241 { -0.001481, 1, 1, 0, 1, 1 },
242 { 0.001417, -1, -1, 0, 1, 1 },
243 { 0.001350, -1, 0, 0, 1, 1 },
244 { 0.001330, 0, 0, -1, 1, 0 },
245 { 0.001106, 0, 3, 0, 1, 0 },
246 { 0.001020, 0, 0, 4, -1, 0 },
247 { 0.000833, 0, -1, 4, 1, 0 },
248 { 0.000781, 0, 1, 0, -3, 0 },
249 { 0.000670, 0, -2, 4, 1, 0 },
250 { 0.000606, 0, 0, 2, -3, 0 },
251 { 0.000597, 0, 2, 2, -1, 0 },
252 { 0.000492, -1, 1, 2, -1, 1 },
253 { 0.000450, 0, 2, -2, -1, 0 },
254 { 0.000439, 0, 3, 0, -1, 0 },
255 { 0.000423, 0, 2, 2, 1, 0 },
256 { 0.000422, 0, -3, 2, -1, 0 },
257 { -0.000367, 1, -1, 2, 1, 1 },
258 { -0.000353, 1, 0, 2, 1, 1 },
259 { 0.000331, 0, 0, 4, 1, 0 },
260 { 0.000317, -1, 1, 2, 1, 1 },
261 { 0.000306, -2, 0, 2, -1, 2 },
262 { -0.000283, 0, 1, 0, 3, 0 } };
263 static int NB = ( sizeof tb / sizeof ( struct term ) );
264
265/*
266** Parallax coeff M M' D F n
267*/
268 static struct term tp[] = { { 0.950724, 0, 0, 0, 0, 0 },
269 { 0.051818, 0, 1, 0, 0, 0 },
270 { 0.009531, 0, -1, 2, 0, 0 },
271 { 0.007843, 0, 0, 2, 0, 0 },
272 { 0.002824, 0, 2, 0, 0, 0 },
273 { 0.000857, 0, 1, 2, 0, 0 },
274 { 0.000533, -1, 0, 2, 0, 1 },
275 { 0.000401, -1, -1, 2, 0, 1 },
276 { 0.000320, -1, 1, 0, 0, 1 },
277 { -0.000271, 0, 0, 1, 0, 0 },
278 { -0.000264, 1, 1, 0, 0, 1 },
279 { -0.000198, 0, -1, 0, 2, 0 },
280 { 0.000173, 0, 3, 0, 0, 0 },
281 { 0.000167, 0, -1, 4, 0, 0 },
282 { -0.000111, 1, 0, 0, 0, 1 },
283 { 0.000103, 0, -2, 4, 0, 0 },
284 { -0.000084, 0, 2, -2, 0, 0 },
285 { -0.000083, 1, 0, 2, 0, 1 },
286 { 0.000079, 0, 2, 2, 0, 0 },
287 { 0.000072, 0, 0, 4, 0, 0 },
288 { 0.000064, -1, 1, 2, 0, 1 },
289 { -0.000063, 1, -1, 2, 0, 1 },
290 { 0.000041, 1, 0, 1, 0, 1 },
291 { 0.000035, -1, 2, 0, 0, 1 },
292 { -0.000033, 0, 3, -2, 0, 0 },
293 { -0.000030, 0, 1, 1, 0, 0 },
294 { -0.000029, 0, 0, -2, 2, 0 },
295 { -0.000029, 1, 2, 0, 0, 1 },
296 { 0.000026, -2, 0, 2, 0, 2 },
297 { -0.000023, 0, 1, -2, 2, 0 },
298 { 0.000019, -1, -1, 4, 0, 1 } };
299 static int NP = ( sizeof tp / sizeof ( struct term ) );
300
301
302
303/* --------------------- */
304/* Execution starts here */
305/* --------------------- */
306
307/* Centuries since J1900 */
308 t = ( date - 15019.5 ) / 36525.0;
309
310/* --------------------- */
311/* Fundamental arguments */
312/* --------------------- */
313
314/* Arguments (radians) and derivatives (radians per Julian century)
315 for the current epoch */
316
317/* Moon's mean longitude */
318 elp = DD2R * dmod ( elp0 + ( elp1 + ( elp2 + elp3 * t ) * t ) * t,
319 360.0 );
320 delp = DD2R * ( elp1 + ( 2.0 * elp2 + 3.0 *elp3 * t ) * t );
321
322/* Sun's mean anomaly */
323 em = DD2R * dmod ( em0 + ( em1 + ( em2 + em3 * t ) * t ) * t, 360.0 );
324 dem = DD2R * ( em1 + ( 2.0 * em2 + 3.0 * em3 * t ) * t );
325
326/* Moon's mean anomaly */
327 emp = DD2R * dmod ( emp0 + ( emp1 + ( emp2 + emp3 * t ) * t ) * t,
328 360.0 );
329 demp = DD2R * ( emp1 + ( 2.0 * emp2 + 3.0 * emp3 * t ) * t );
330
331/* Moon's mean elongation */
332 d = DD2R * dmod ( d0 + ( d1 + ( d2 + d3 * t ) * t ) * t, 360.0 );
333 dd = DD2R * ( d1 + ( 2.0 * d2 + 3.0 * d3 * t ) * t );
334
335/* Mean distance of the Moon from its ascending node */
336 f = DD2R * dmod ( f0 + ( f1 + ( f2 + f3 * t ) * t ) * t, 360.0 );
337 df = DD2R * ( f1 + ( 2.0 * f2 + 3.0 * f3 * t ) * t );
338
339/* Longitude of the Moon's ascending node */
340 om = DD2R * dmod ( om0 + ( om1 + ( om2 + om3 * t ) * t ) * t, 360.0 );
341 dom = DD2R * ( om1 + ( 2.0 * om2 + 3.0 * om3 * t ) * t );
342 sinom = sin ( om );
343 cosom = cos ( om );
344 domcom = dom * cosom;
345
346/* Add the periodic variations */
347 theta = DD2R * ( pa0 + pa1 * t );
348 wa = sin ( theta );
349 dwa = DD2R * pa1 * cos ( theta );
350 theta = DD2R * ( pe0 + ( pe1 + pe2 * t ) * t );
351 wb = pec * sin ( theta );
352 dwb = DD2R * pec*( pe1 + 2.0 * pe2 * t ) * cos ( theta );
353 elp += DD2R * ( pac * wa + wb + pfc * sinom );
354 delp += DD2R * ( pac * dwa + dwb + pfc * domcom );
355 em += DD2R * pbc * wa;
356 dem += DD2R * pbc * dwa;
357 emp += DD2R * ( pcc * wa + wb + pgc * sinom );
358 demp += DD2R * ( pcc * dwa + dwb + pgc * domcom );
359 d += DD2R * ( pdc * wa + wb + phc * sinom );
360 dd += DD2R * ( pdc * dwa + dwb + phc * domcom );
361 wom = om + DD2R * ( pj0 + pj1 * t );
362 dwom = dom + DD2R * pj1;
363 sinwom = sin ( wom );
364 coswom = cos ( wom );
365 f += DD2R * ( wb + pic * sinom + pjc * sinwom );
366 df += DD2R * ( dwb + pic * domcom + pjc * dwom * coswom );
367
368/* E-factor, and square */
369 e = 1.0 + ( e1 + e2 * t ) * t;
370 de = e1 + 2.0 * e2 * t;
371 esq = e * e;
372 desq = 2.0 * e * de;
373
374/* ----------------- */
375/* Series expansions */
376/* ----------------- */
377
378/* Longitude */
379 v = 0.0;
380 dv = 0.0;
381 for ( n = NL -1; n >= 0; n-- ) {
382 coeff = tl[n].coef;
383 emn = (double) tl[n].nem;
384 empn = (double) tl[n].nemp;
385 dn = (double) tl[n].nd;
386 fn = (double) tl[n].nf;
387 i = tl[n].ne;
388 if ( i == 0 ) {
389 en = 1.0;
390 den = 0.0;
391 } else if ( i == 1 ) {
392 en = e;
393 den = de;
394 } else {
395 en = esq;
396 den = desq;
397 }
398 theta = emn * em + empn * emp + dn * d + fn * f;
399 dtheta = emn * dem + empn * demp + dn * dd + fn * df;
400 ftheta = sin ( theta );
401 v += coeff * ftheta * en;
402 dv += coeff * ( cos ( theta ) * dtheta * en + ftheta * den );
403 }
404 el = elp + DD2R * v;
405 del = ( delp + DD2R * dv ) / CJ;
406
407/* Latitude */
408 v = 0.0;
409 dv = 0.0;
410 for ( n = NB - 1; n >= 0; n-- ) {
411 coeff = tb[n].coef;
412 emn = (double) tb[n].nem;
413 empn = (double) tb[n].nemp;
414 dn = (double) tb[n].nd;
415 fn = (double) tb[n].nf;
416 i = tb[n].ne;
417 if ( i == 0 ) {
418 en = 1.0;
419 den = 0.0;
420 } else if ( i == 1 ) {
421 en = e;
422 den = de;
423 } else {
424 en = esq;
425 den = desq;
426 }
427 theta = emn * em + empn * emp + dn * d + fn * f;
428 dtheta = emn * dem + empn * demp + dn * dd + fn * df;
429 ftheta = sin ( theta );
430 v += coeff * ftheta * en;
431 dv += coeff * ( cos ( theta ) * dtheta * en + ftheta * den );
432 }
433 bf = 1.0 - cw1 * cosom - cw2 * coswom;
434 dbf = cw1 * dom * sinom + cw2 * dwom * sinwom;
435 b = DD2R * v * bf;
436 db = DD2R * ( dv * bf + v * dbf ) / CJ;
437
438/* Parallax */
439 v = 0.0;
440 dv = 0.0;
441 for ( n = NP - 1; n >= 0; n-- ) {
442 coeff = tp[n].coef;
443 emn = (double) tp[n].nem;
444 empn = (double) tp[n].nemp;
445 dn = (double) tp[n].nd;
446 fn = (double) tp[n].nf;
447 i = tp[n].ne;
448 if ( i == 0 ) {
449 en = 1.0;
450 den = 0.0;
451 } else if ( i == 1 ) {
452 en = e;
453 den = de;
454 } else {
455 en = esq;
456 den = desq;
457 }
458 theta = emn * em + empn * emp + dn * d + fn * f;
459 dtheta = emn * dem + empn * demp + dn * dd + fn * df;
460 ftheta = cos ( theta );
461 v += coeff * ftheta * en;
462 dv += coeff* ( - sin ( theta ) * dtheta * en + ftheta * den );
463 }
464 p = DD2R * v;
465 dp = DD2R * dv / CJ;
466
467/* ------------------------------ */
468/* Transformation into final form */
469/* ------------------------------ */
470
471/* Parallax to distance (AU, AU/sec) */
472 sp = sin ( p );
473 r = ERADAU / sp;
474 dr = - r * dp * cos ( p ) / sp;
475
476/* Longitude, latitude to x, y, z (AU) */
477 sel = sin ( el );
478 cel = cos ( el );
479 sb = sin ( b );
480 cb = cos ( b );
481 rcb = r * cb;
482 rbd = r * db;
483 w = rbd * sb - cb * dr;
484 x = rcb * cel;
485 y = rcb * sel;
486 z = r * sb;
487 xd = - y * del - w * cel;
488 yd = x * del - w * sel;
489 zd = rbd * cb + sb * dr;
490
491/* Julian centuries since J2000 */
492 t = ( date - 51544.5 ) / 36525.0;
493
494/* Fricke equinox correction */
495 epj = 2000.0 + t * 100.0;
496 eqcor = DS2R * ( 0.035 + 0.00085 * ( epj - B1950 ) );
497
498/* Mean obliquity (IAU 1976) */
499 eps = DAS2R *
500 ( 84381.448 + ( - 46.8150 + ( - 0.00059 + 0.001813 * t ) * t ) * t );
501
502/* To the equatorial system, mean of date, FK5 system */
503 sineps = sin ( eps );
504 coseps = cos ( eps );
505 es = eqcor * sineps;
506 ec = eqcor * coseps;
507 pv[0] = x - ec * y + es * z;
508 pv[1] = eqcor * x + y * coseps - z * sineps;
509 pv[2] = y * sineps + z * coseps;
510 pv[3] = xd - ec * yd + es * zd;
511 pv[4] = eqcor * xd + yd * coseps - zd * sineps;
512 pv[5] = yd * sineps + zd * coseps;
513}
Note: See TracBrowser for help on using the repository browser.