FACT++  1.0
palDmoon.c
Go to the documentation of this file.
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 
95 void 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 }
PAL__DS2R
static const double PAL__DS2R
Definition: palmac.h:93
PAL__DAS2R
static const double PAL__DAS2R
Definition: palmac.h:78
NL
#define NL
PAL__DD2R
static const double PAL__DD2R
Definition: palmac.h:72
palDmoon
void palDmoon(double date, double pv[6])
Definition: palDmoon.c:95
NB
#define NB
NP
#define NP