source: trunk/FACT++/pal/palDmoon.c@ 19633

Last change on this file since 19633 was 18347, checked in by tbretz, 9 years ago
File size: 15.3 KB
Line 
1/*
2*+
3* Name:
4* palDmoon
5
6* Purpose:
7* Approximate geocentric position and velocity of the Moon
8
9* Language:
10* Starlink ANSI C
11
12* Type of Module:
13* Library routine
14
15* Invocation:
16* void palDmoon( double date, double pv[6] );
17
18* Arguments:
19* date = double (Given)
20* TDB as a Modified Julian Date (JD-2400000.5)
21* pv = double [6] (Returned)
22* Moon x,y,z,xdot,ydot,zdot, mean equator and
23* equinox of date (AU, AU/s)
24
25* Description:
26* Calculate the approximate geocentric position of the Moon
27* using a full implementation of the algorithm published by
28* Meeus (l'Astronomie, June 1984, p348).
29
30* Authors:
31* TIMJ: Tim Jenness (JAC, Hawaii)
32* PTW: Patrick T. Wallace
33* {enter_new_authors_here}
34
35* Notes:
36* - Meeus quotes accuracies of 10 arcsec in longitude, 3 arcsec in
37* latitude and 0.2 arcsec in HP (equivalent to about 20 km in
38* distance). Comparison with JPL DE200 over the interval
39* 1960-2025 gives RMS errors of 3.7 arcsec and 83 mas/hour in
40* longitude, 2.3 arcsec and 48 mas/hour in latitude, 11 km
41* and 81 mm/s in distance. The maximum errors over the same
42* interval are 18 arcsec and 0.50 arcsec/hour in longitude,
43* 11 arcsec and 0.24 arcsec/hour in latitude, 40 km and 0.29 m/s
44* in distance.
45* - The original algorithm is expressed in terms of the obsolete
46* timescale Ephemeris Time. Either TDB or TT can be used, but
47* not UT without incurring significant errors (30 arcsec at
48* the present time) due to the Moon's 0.5 arcsec/sec movement.
49* - The algorithm is based on pre IAU 1976 standards. However,
50* the result has been moved onto the new (FK5) equinox, an
51* adjustment which is in any case much smaller than the
52* intrinsic accuracy of the procedure.
53* - Velocity is obtained by a complete analytical differentiation
54* of the Meeus model.
55
56* History:
57* 2012-03-07 (TIMJ):
58* Initial version based on a direct port of the SLA/F code.
59* Adapted with permission from the Fortran SLALIB library.
60* {enter_further_changes_here}
61
62* Copyright:
63* Copyright (C) 1998 Rutherford Appleton Laboratory
64* Copyright (C) 2012 Science and Technology Facilities Council.
65* All Rights Reserved.
66
67* Licence:
68* This program is free software; you can redistribute it and/or
69* modify it under the terms of the GNU General Public License as
70* published by the Free Software Foundation; either version 3 of
71* the License, or (at your option) any later version.
72*
73* This program is distributed in the hope that it will be
74* useful, but WITHOUT ANY WARRANTY; without even the implied
75* warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
76* PURPOSE. See the GNU General Public License for more details.
77*
78* You should have received a copy of the GNU General Public License
79* along with this program; if not, write to the Free Software
80* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
81* MA 02110-1301, USA.
82
83* Bugs:
84* {note_any_bugs_here}
85*-
86*/
87
88#include "pal.h"
89#include "pal1sofa.h"
90#include "palmac.h"
91
92/* Autoconf can give us -DPIC */
93#undef PIC
94
95void palDmoon( double date, double pv[6] ) {
96
97 /* Seconds per Julian century (86400*36525) */
98 const double CJ = 3155760000.0;
99
100 /* Julian epoch of B1950 */
101 const double B1950 = 1949.9997904423;
102
103 /* Earth equatorial radius in AU ( = 6378.137 / 149597870 ) */
104 const double ERADAU=4.2635212653763e-5;
105
106 double T,THETA,SINOM,COSOM,DOMCOM,WA,DWA,WB,DWB,WOM,
107 DWOM,SINWOM,COSWOM,V,DV,COEFF,EMN,EMPN,DN,FN,EN,
108 DEN,DTHETA,FTHETA,EL,DEL,B,DB,BF,DBF,P,DP,SP,R,
109 DR,X,Y,Z,XD,YD,ZD,SEL,CEL,SB,CB,RCB,RBD,W,EPJ,
110 EQCOR,EPS,SINEPS,COSEPS,ES,EC;
111 double ELP, DELP;
112 double EM, DEM, EMP, DEMP, D, DD, F, DF, OM, DOM, E, DESQ, ESQ, DE;
113 int N,I;
114
115 /*
116 * Coefficients for fundamental arguments
117 *
118 * at J1900: T**0, T**1, T**2, T**3
119 * at epoch: T**0, T**1
120 *
121 * Units are degrees for position and Julian centuries for time
122 *
123 */
124
125 /* Moon's mean longitude */
126 const double ELP0=270.434164;
127 const double ELP1=481267.8831;
128 const double ELP2=-0.001133;
129 const double ELP3=0.0000019;
130
131 /* Sun's mean anomaly */
132 const double EM0=358.475833;
133 const double EM1=35999.0498;
134 const double EM2=-0.000150;
135 const double EM3=-0.0000033;
136
137 /* Moon's mean anomaly */
138 const double EMP0=296.104608;
139 const double EMP1=477198.8491;
140 const double EMP2=0.009192;
141 const double EMP3=0.0000144;
142
143 /* Moon's mean elongation */
144 const double D0=350.737486;
145 const double D1=445267.1142;
146 const double D2=-0.001436;
147 const double D3=0.0000019;
148
149 /* Mean distance of the Moon from its ascending node */
150 const double F0=11.250889;
151 const double F1=483202.0251;
152 const double F2=-0.003211;
153 const double F3=-0.0000003;
154
155 /* Longitude of the Moon's ascending node */
156 const double OM0=259.183275;
157 const double OM1=-1934.1420;
158 const double OM2=0.002078;
159 const double OM3=0.0000022;
160
161 /* Coefficients for (dimensionless) E factor */
162 const double E1=-0.002495;
163 const double E2=-0.00000752;
164
165 /* Coefficients for periodic variations etc */
166 const double PAC=0.000233;
167 const double PA0=51.2;
168 const double PA1=20.2;
169 const double PBC=-0.001778;
170 const double PCC=0.000817;
171 const double PDC=0.002011;
172 const double PEC=0.003964;
173 const double PE0=346.560;
174 const double PE1=132.870;
175 const double PE2=-0.0091731;
176 const double PFC=0.001964;
177 const double PGC=0.002541;
178 const double PHC=0.001964;
179 const double PIC=-0.024691;
180 const double PJC=-0.004328;
181 const double PJ0=275.05;
182 const double PJ1=-2.30;
183 const double CW1=0.0004664;
184 const double CW2=0.0000754;
185
186 /*
187 * Coefficients for Moon position
188 *
189 * Tx(N) = coefficient of L, B or P term (deg)
190 * ITx(N,1-5) = coefficients of M, M', D, F, E**n in argument
191 */
192#define NL 50
193#define NB 45
194#define NP 31
195
196 /*
197 * Longitude
198 */
199 const double TL[NL] = {
200 6.28875,1.274018,.658309,.213616,-.185596,
201 -.114336,.058793,.057212,.05332,.045874,.041024,-.034718,-.030465,
202 .015326,-.012528,-.01098,.010674,.010034,.008548,-.00791,-.006783,
203 .005162,.005,.004049,.003996,.003862,.003665,.002695,.002602,
204 .002396,-.002349,.002249,-.002125,-.002079,.002059,-.001773,
205 -.001595,.00122,-.00111,8.92e-4,-8.11e-4,7.61e-4,7.17e-4,7.04e-4,
206 6.93e-4,5.98e-4,5.5e-4,5.38e-4,5.21e-4,4.86e-4
207 };
208
209 const int ITL[NL][5] = {
210 /* M M' D F n */
211 { +0, +1, +0, +0, 0 },
212 { +0, -1, +2, +0, 0 },
213 { +0, +0, +2, +0, 0 },
214 { +0, +2, +0, +0, 0 },
215 { +1, +0, +0, +0, 1 },
216 { +0, +0, +0, +2, 0 },
217 { +0, -2, +2, +0, 0 },
218 { -1, -1, +2, +0, 1 },
219 { +0, +1, +2, +0, 0 },
220 { -1, +0, +2, +0, 1 },
221 { -1, +1, +0, +0, 1 },
222 { +0, +0, +1, +0, 0 },
223 { +1, +1, +0, +0, 1 },
224 { +0, +0, +2, -2, 0 },
225 { +0, +1, +0, +2, 0 },
226 { +0, -1, +0, +2, 0 },
227 { +0, -1, +4, +0, 0 },
228 { +0, +3, +0, +0, 0 },
229 { +0, -2, +4, +0, 0 },
230 { +1, -1, +2, +0, 1 },
231 { +1, +0, +2, +0, 1 },
232 { +0, +1, -1, +0, 0 },
233 { +1, +0, +1, +0, 1 },
234 { -1, +1, +2, +0, 1 },
235 { +0, +2, +2, +0, 0 },
236 { +0, +0, +4, +0, 0 },
237 { +0, -3, +2, +0, 0 },
238 { -1, +2, +0, +0, 1 },
239 { +0, +1, -2, -2, 0 },
240 { -1, -2, +2, +0, 1 },
241 { +0, +1, +1, +0, 0 },
242 { -2, +0, +2, +0, 2 },
243 { +1, +2, +0, +0, 1 },
244 { +2, +0, +0, +0, 2 },
245 { -2, -1, +2, +0, 2 },
246 { +0, +1, +2, -2, 0 },
247 { +0, +0, +2, +2, 0 },
248 { -1, -1, +4, +0, 1 },
249 { +0, +2, +0, +2, 0 },
250 { +0, +1, -3, +0, 0 },
251 { +1, +1, +2, +0, 1 },
252 { -1, -2, +4, +0, 1 },
253 { -2, +1, +0, +0, 2 },
254 { -2, +1, -2, +0, 2 },
255 { +1, -2, +2, +0, 1 },
256 { -1, +0, +2, -2, 1 },
257 { +0, +1, +4, +0, 0 },
258 { +0, +4, +0, +0, 0 },
259 { -1, +0, +4, +0, 1 },
260 { +0, +2, -1, +0, 0 }
261 };
262
263 /*
264 * Latitude
265 */
266 const double TB[NB] = {
267 5.128189,.280606,.277693,.173238,.055413,
268 .046272,.032573,.017198,.009267,.008823,.008247,.004323,.0042,
269 .003372,.002472,.002222,.002072,.001877,.001828,-.001803,-.00175,
270 .00157,-.001487,-.001481,.001417,.00135,.00133,.001106,.00102,
271 8.33e-4,7.81e-4,6.7e-4,6.06e-4,5.97e-4,4.92e-4,4.5e-4,4.39e-4,
272 4.23e-4,4.22e-4,-3.67e-4,-3.53e-4,3.31e-4,3.17e-4,3.06e-4,
273 -2.83e-4
274 };
275
276 const int ITB[NB][5] = {
277 /* M M' D F n */
278 { +0, +0, +0, +1, 0 },
279 { +0, +1, +0, +1, 0 },
280 { +0, +1, +0, -1, 0 },
281 { +0, +0, +2, -1, 0 },
282 { +0, -1, +2, +1, 0 },
283 { +0, -1, +2, -1, 0 },
284 { +0, +0, +2, +1, 0 },
285 { +0, +2, +0, +1, 0 },
286 { +0, +1, +2, -1, 0 },
287 { +0, +2, +0, -1, 0 },
288 { -1, +0, +2, -1, 1 },
289 { +0, -2, +2, -1, 0 },
290 { +0, +1, +2, +1, 0 },
291 { -1, +0, -2, +1, 1 },
292 { -1, -1, +2, +1, 1 },
293 { -1, +0, +2, +1, 1 },
294 { -1, -1, +2, -1, 1 },
295 { -1, +1, +0, +1, 1 },
296 { +0, -1, +4, -1, 0 },
297 { +1, +0, +0, +1, 1 },
298 { +0, +0, +0, +3, 0 },
299 { -1, +1, +0, -1, 1 },
300 { +0, +0, +1, +1, 0 },
301 { +1, +1, +0, +1, 1 },
302 { -1, -1, +0, +1, 1 },
303 { -1, +0, +0, +1, 1 },
304 { +0, +0, -1, +1, 0 },
305 { +0, +3, +0, +1, 0 },
306 { +0, +0, +4, -1, 0 },
307 { +0, -1, +4, +1, 0 },
308 { +0, +1, +0, -3, 0 },
309 { +0, -2, +4, +1, 0 },
310 { +0, +0, +2, -3, 0 },
311 { +0, +2, +2, -1, 0 },
312 { -1, +1, +2, -1, 1 },
313 { +0, +2, -2, -1, 0 },
314 { +0, +3, +0, -1, 0 },
315 { +0, +2, +2, +1, 0 },
316 { +0, -3, +2, -1, 0 },
317 { +1, -1, +2, +1, 1 },
318 { +1, +0, +2, +1, 1 },
319 { +0, +0, +4, +1, 0 },
320 { -1, +1, +2, +1, 1 },
321 { -2, +0, +2, -1, 2 },
322 { +0, +1, +0, +3, 0 }
323 };
324
325 /*
326 * Parallax
327 */
328 const double TP[NP] = {
329 .950724,.051818,.009531,.007843,.002824,
330 8.57e-4,5.33e-4,4.01e-4,3.2e-4,-2.71e-4,-2.64e-4,-1.98e-4,1.73e-4,
331 1.67e-4,-1.11e-4,1.03e-4,-8.4e-5,-8.3e-5,7.9e-5,7.2e-5,6.4e-5,
332 -6.3e-5,4.1e-5,3.5e-5,-3.3e-5,-3e-5,-2.9e-5,-2.9e-5,2.6e-5,
333 -2.3e-5,1.9e-5
334 };
335
336 const int ITP[NP][5] = {
337 /* M M' D F n */
338 { +0, +0, +0, +0, 0 },
339 { +0, +1, +0, +0, 0 },
340 { +0, -1, +2, +0, 0 },
341 { +0, +0, +2, +0, 0 },
342 { +0, +2, +0, +0, 0 },
343 { +0, +1, +2, +0, 0 },
344 { -1, +0, +2, +0, 1 },
345 { -1, -1, +2, +0, 1 },
346 { -1, +1, +0, +0, 1 },
347 { +0, +0, +1, +0, 0 },
348 { +1, +1, +0, +0, 1 },
349 { +0, -1, +0, +2, 0 },
350 { +0, +3, +0, +0, 0 },
351 { +0, -1, +4, +0, 0 },
352 { +1, +0, +0, +0, 1 },
353 { +0, -2, +4, +0, 0 },
354 { +0, +2, -2, +0, 0 },
355 { +1, +0, +2, +0, 1 },
356 { +0, +2, +2, +0, 0 },
357 { +0, +0, +4, +0, 0 },
358 { -1, +1, +2, +0, 1 },
359 { +1, -1, +2, +0, 1 },
360 { +1, +0, +1, +0, 1 },
361 { -1, +2, +0, +0, 1 },
362 { +0, +3, -2, +0, 0 },
363 { +0, +1, +1, +0, 0 },
364 { +0, +0, -2, +2, 0 },
365 { +1, +2, +0, +0, 1 },
366 { -2, +0, +2, +0, 2 },
367 { +0, +1, -2, +2, 0 },
368 { -1, -1, +4, +0, 1 }
369 };
370
371 /* Centuries since J1900 */
372 T=(date-15019.5)/36525.;
373
374 /*
375 * Fundamental arguments (radians) and derivatives (radians per
376 * Julian century) for the current epoch
377 */
378
379 /* Moon's mean longitude */
380 ELP=PAL__DD2R*fmod(ELP0+(ELP1+(ELP2+ELP3*T)*T)*T,360.);
381 DELP=PAL__DD2R*(ELP1+(2.*ELP2+3*ELP3*T)*T);
382
383 /* Sun's mean anomaly */
384 EM=PAL__DD2R*fmod(EM0+(EM1+(EM2+EM3*T)*T)*T,360.);
385 DEM=PAL__DD2R*(EM1+(2.*EM2+3*EM3*T)*T);
386
387 /* Moon's mean anomaly */
388 EMP=PAL__DD2R*fmod(EMP0+(EMP1+(EMP2+EMP3*T)*T)*T,360.);
389 DEMP=PAL__DD2R*(EMP1+(2.*EMP2+3*EMP3*T)*T);
390
391 /* Moon's mean elongation */
392 D=PAL__DD2R*fmod(D0+(D1+(D2+D3*T)*T)*T,360.);
393 DD=PAL__DD2R*(D1+(2.*D2+3.*D3*T)*T);
394
395 /* Mean distance of the Moon from its ascending node */
396 F=PAL__DD2R*fmod(F0+(F1+(F2+F3*T)*T)*T,360.);
397 DF=PAL__DD2R*(F1+(2.*F2+3.*F3*T)*T);
398
399 /* Longitude of the Moon's ascending node */
400 OM=PAL__DD2R*fmod(OM0+(OM1+(OM2+OM3*T)*T)*T,360.);
401 DOM=PAL__DD2R*(OM1+(2.*OM2+3.*OM3*T)*T);
402 SINOM=sin(OM);
403 COSOM=cos(OM);
404 DOMCOM=DOM*COSOM;
405
406 /* Add the periodic variations */
407 THETA=PAL__DD2R*(PA0+PA1*T);
408 WA=sin(THETA);
409 DWA=PAL__DD2R*PA1*cos(THETA);
410 THETA=PAL__DD2R*(PE0+(PE1+PE2*T)*T);
411 WB=PEC*sin(THETA);
412 DWB=PAL__DD2R*PEC*(PE1+2.*PE2*T)*cos(THETA);
413 ELP=ELP+PAL__DD2R*(PAC*WA+WB+PFC*SINOM);
414 DELP=DELP+PAL__DD2R*(PAC*DWA+DWB+PFC*DOMCOM);
415 EM=EM+PAL__DD2R*PBC*WA;
416 DEM=DEM+PAL__DD2R*PBC*DWA;
417 EMP=EMP+PAL__DD2R*(PCC*WA+WB+PGC*SINOM);
418 DEMP=DEMP+PAL__DD2R*(PCC*DWA+DWB+PGC*DOMCOM);
419 D=D+PAL__DD2R*(PDC*WA+WB+PHC*SINOM);
420 DD=DD+PAL__DD2R*(PDC*DWA+DWB+PHC*DOMCOM);
421 WOM=OM+PAL__DD2R*(PJ0+PJ1*T);
422 DWOM=DOM+PAL__DD2R*PJ1;
423 SINWOM=sin(WOM);
424 COSWOM=cos(WOM);
425 F=F+PAL__DD2R*(WB+PIC*SINOM+PJC*SINWOM);
426 DF=DF+PAL__DD2R*(DWB+PIC*DOMCOM+PJC*DWOM*COSWOM);
427
428 /* E-factor, and square */
429 E=1.+(E1+E2*T)*T;
430 DE=E1+2.*E2*T;
431 ESQ=E*E;
432 DESQ=2.*E*DE;
433
434 /*
435 * Series expansions
436 */
437
438 /* Longitude */
439 V=0.;
440 DV=0.;
441 for (N=NL-1; N>=0; N--) { /* DO N=NL, 1, -1 */
442 COEFF=TL[N];
443 EMN=(double)(ITL[N][0]);
444 EMPN=(double)(ITL[N][1]);
445 DN=(double)(ITL[N][2]);
446 FN=(double)(ITL[N][3]);
447 I=ITL[N][4];
448 if (I == 0) {
449 EN=1.;
450 DEN=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=sin(THETA);
461 V=V+COEFF*FTHETA*EN;
462 DV=DV+COEFF*(cos(THETA)*DTHETA*EN+FTHETA*DEN);
463 }
464 EL=ELP+PAL__DD2R*V;
465 DEL=(DELP+PAL__DD2R*DV)/CJ;
466
467 /* Latitude */
468 V=0.;
469 DV=0.;
470 for (N=NB-1; N>=0; N--) { /* DO N=NB,1,-1 */
471 COEFF=TB[N];
472 EMN=(double)(ITB[N][0]);
473 EMPN=(double)(ITB[N][1]);
474 DN=(double)(ITB[N][2]);
475 FN=(double)(ITB[N][3]);
476 I=ITB[N][4];
477 if (I == 0 ) {
478 EN=1.;
479 DEN=0.;
480 } else if (I == 1) {
481 EN=E;
482 DEN=DE;
483 } else {
484 EN=ESQ;
485 DEN=DESQ;
486 }
487 THETA=EMN*EM+EMPN*EMP+DN*D+FN*F;
488 DTHETA=EMN*DEM+EMPN*DEMP+DN*DD+FN*DF;
489 FTHETA=sin(THETA);
490 V=V+COEFF*FTHETA*EN;
491 DV=DV+COEFF*(cos(THETA)*DTHETA*EN+FTHETA*DEN);
492 }
493 BF=1.-CW1*COSOM-CW2*COSWOM;
494 DBF=CW1*DOM*SINOM+CW2*DWOM*SINWOM;
495 B=PAL__DD2R*V*BF;
496 DB=PAL__DD2R*(DV*BF+V*DBF)/CJ;
497
498 /* Parallax */
499 V=0.;
500 DV=0.;
501 for (N=NP-1; N>=0; N--) { /* DO N=NP,1,-1 */
502 COEFF=TP[N];
503 EMN=(double)(ITP[N][0]);
504 EMPN=(double)(ITP[N][1]);
505 DN=(double)(ITP[N][2]);
506 FN=(double)(ITP[N][3]);
507 I=ITP[N][4];
508 if (I == 0) {
509 EN=1.;
510 DEN=0.;
511 } else if (I == 1) {
512 EN=E;
513 DEN=DE;
514 } else {
515 EN=ESQ;
516 DEN=DESQ;
517 }
518 THETA=EMN*EM+EMPN*EMP+DN*D+FN*F;
519 DTHETA=EMN*DEM+EMPN*DEMP+DN*DD+FN*DF;
520 FTHETA=cos(THETA);
521 V=V+COEFF*FTHETA*EN;
522 DV=DV+COEFF*(-sin(THETA)*DTHETA*EN+FTHETA*DEN);
523 }
524 P=PAL__DD2R*V;
525 DP=PAL__DD2R*DV/CJ;
526
527 /*
528 * Transformation into final form
529 */
530
531 /* Parallax to distance (AU, AU/sec) */
532 SP=sin(P);
533 R=ERADAU/SP;
534 DR=-R*DP*cos(P)/SP;
535
536 /* Longitude, latitude to x,y,z (AU) */
537 SEL=sin(EL);
538 CEL=cos(EL);
539 SB=sin(B);
540 CB=cos(B);
541 RCB=R*CB;
542 RBD=R*DB;
543 W=RBD*SB-CB*DR;
544 X=RCB*CEL;
545 Y=RCB*SEL;
546 Z=R*SB;
547 XD=-Y*DEL-W*CEL;
548 YD=X*DEL-W*SEL;
549 ZD=RBD*CB+SB*DR;
550
551 /* Julian centuries since J2000 */
552 T=(date-51544.5)/36525.;
553
554 /* Fricke equinox correction */
555 EPJ=2000.+T*100.;
556 EQCOR=PAL__DS2R*(0.035+0.00085*(EPJ-B1950));
557
558 /* Mean obliquity (IAU 1976) */
559 EPS=PAL__DAS2R*(84381.448+(-46.8150+(-0.00059+0.001813*T)*T)*T);
560
561 /* To the equatorial system, mean of date, FK5 system */
562 SINEPS=sin(EPS);
563 COSEPS=cos(EPS);
564 ES=EQCOR*SINEPS;
565 EC=EQCOR*COSEPS;
566 pv[0]=X-EC*Y+ES*Z;
567 pv[1]=EQCOR*X+Y*COSEPS-Z*SINEPS;
568 pv[2]=Y*SINEPS+Z*COSEPS;
569 pv[3]=XD-EC*YD+ES*ZD;
570 pv[4]=EQCOR*XD+YD*COSEPS-ZD*SINEPS;
571 pv[5]=YD*SINEPS+ZD*COSEPS;
572
573}
Note: See TracBrowser for help on using the repository browser.