File Coverage

palsrc/palDmoon.c
Criterion Covered Total %
statement 128 128 100.0
branch 18 18 100.0
condition n/a
subroutine n/a
pod n/a
total 146 146 100.0


line stmt bran cond sub pod time code
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 190           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 190           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 190           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 190           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 190           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 190           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 190           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 190           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 190           ELP=PAL__DD2R*fmod(ELP0+(ELP1+(ELP2+ELP3*T)*T)*T,360.);
381 190           DELP=PAL__DD2R*(ELP1+(2.*ELP2+3*ELP3*T)*T);
382              
383             /* Sun's mean anomaly */
384 190           EM=PAL__DD2R*fmod(EM0+(EM1+(EM2+EM3*T)*T)*T,360.);
385 190           DEM=PAL__DD2R*(EM1+(2.*EM2+3*EM3*T)*T);
386              
387             /* Moon's mean anomaly */
388 190           EMP=PAL__DD2R*fmod(EMP0+(EMP1+(EMP2+EMP3*T)*T)*T,360.);
389 190           DEMP=PAL__DD2R*(EMP1+(2.*EMP2+3*EMP3*T)*T);
390              
391             /* Moon's mean elongation */
392 190           D=PAL__DD2R*fmod(D0+(D1+(D2+D3*T)*T)*T,360.);
393 190           DD=PAL__DD2R*(D1+(2.*D2+3.*D3*T)*T);
394              
395             /* Mean distance of the Moon from its ascending node */
396 190           F=PAL__DD2R*fmod(F0+(F1+(F2+F3*T)*T)*T,360.);
397 190           DF=PAL__DD2R*(F1+(2.*F2+3.*F3*T)*T);
398              
399             /* Longitude of the Moon's ascending node */
400 190           OM=PAL__DD2R*fmod(OM0+(OM1+(OM2+OM3*T)*T)*T,360.);
401 190           DOM=PAL__DD2R*(OM1+(2.*OM2+3.*OM3*T)*T);
402 190           SINOM=sin(OM);
403 190           COSOM=cos(OM);
404 190           DOMCOM=DOM*COSOM;
405              
406             /* Add the periodic variations */
407 190           THETA=PAL__DD2R*(PA0+PA1*T);
408 190           WA=sin(THETA);
409 190           DWA=PAL__DD2R*PA1*cos(THETA);
410 190           THETA=PAL__DD2R*(PE0+(PE1+PE2*T)*T);
411 190           WB=PEC*sin(THETA);
412 190           DWB=PAL__DD2R*PEC*(PE1+2.*PE2*T)*cos(THETA);
413 190           ELP=ELP+PAL__DD2R*(PAC*WA+WB+PFC*SINOM);
414 190           DELP=DELP+PAL__DD2R*(PAC*DWA+DWB+PFC*DOMCOM);
415 190           EM=EM+PAL__DD2R*PBC*WA;
416 190           DEM=DEM+PAL__DD2R*PBC*DWA;
417 190           EMP=EMP+PAL__DD2R*(PCC*WA+WB+PGC*SINOM);
418 190           DEMP=DEMP+PAL__DD2R*(PCC*DWA+DWB+PGC*DOMCOM);
419 190           D=D+PAL__DD2R*(PDC*WA+WB+PHC*SINOM);
420 190           DD=DD+PAL__DD2R*(PDC*DWA+DWB+PHC*DOMCOM);
421 190           WOM=OM+PAL__DD2R*(PJ0+PJ1*T);
422 190           DWOM=DOM+PAL__DD2R*PJ1;
423 190           SINWOM=sin(WOM);
424 190           COSWOM=cos(WOM);
425 190           F=F+PAL__DD2R*(WB+PIC*SINOM+PJC*SINWOM);
426 190           DF=DF+PAL__DD2R*(DWB+PIC*DOMCOM+PJC*DWOM*COSWOM);
427              
428             /* E-factor, and square */
429 190           E=1.+(E1+E2*T)*T;
430 190           DE=E1+2.*E2*T;
431 190           ESQ=E*E;
432 190           DESQ=2.*E*DE;
433              
434             /*
435             * Series expansions
436             */
437              
438             /* Longitude */
439             V=0.;
440             DV=0.;
441 9690 100         for (N=NL-1; N>=0; N--) { /* DO N=NL, 1, -1 */
442 9500           COEFF=TL[N];
443 9500           EMN=(double)(ITL[N][0]);
444 9500           EMPN=(double)(ITL[N][1]);
445 9500           DN=(double)(ITL[N][2]);
446 9500           FN=(double)(ITL[N][3]);
447 9500           I=ITL[N][4];
448 9500 100         if (I == 0) {
449             EN=1.;
450             DEN=0.;
451 4370 100         } else if (I == 1) {
452             EN=E;
453             DEN=DE;
454             } else {
455             EN=ESQ;
456             DEN=DESQ;
457             }
458 9500           THETA=EMN*EM+EMPN*EMP+DN*D+FN*F;
459 9500           DTHETA=EMN*DEM+EMPN*DEMP+DN*DD+FN*DF;
460 9500           FTHETA=sin(THETA);
461 9500           V=V+COEFF*FTHETA*EN;
462 9500           DV=DV+COEFF*(cos(THETA)*DTHETA*EN+FTHETA*DEN);
463             }
464 190           EL=ELP+PAL__DD2R*V;
465 190           DEL=(DELP+PAL__DD2R*DV)/CJ;
466              
467             /* Latitude */
468             V=0.;
469             DV=0.;
470 8740 100         for (N=NB-1; N>=0; N--) { /* DO N=NB,1,-1 */
471 8550           COEFF=TB[N];
472 8550           EMN=(double)(ITB[N][0]);
473 8550           EMPN=(double)(ITB[N][1]);
474 8550           DN=(double)(ITB[N][2]);
475 8550           FN=(double)(ITB[N][3]);
476 8550           I=ITB[N][4];
477 8550 100         if (I == 0 ) {
478             EN=1.;
479             DEN=0.;
480 3040 100         } else if (I == 1) {
481             EN=E;
482             DEN=DE;
483             } else {
484             EN=ESQ;
485             DEN=DESQ;
486             }
487 8550           THETA=EMN*EM+EMPN*EMP+DN*D+FN*F;
488 8550           DTHETA=EMN*DEM+EMPN*DEMP+DN*DD+FN*DF;
489 8550           FTHETA=sin(THETA);
490 8550           V=V+COEFF*FTHETA*EN;
491 8550           DV=DV+COEFF*(cos(THETA)*DTHETA*EN+FTHETA*DEN);
492             }
493 190           BF=1.-CW1*COSOM-CW2*COSWOM;
494 190           DBF=CW1*DOM*SINOM+CW2*DWOM*SINWOM;
495 190           B=PAL__DD2R*V*BF;
496 190           DB=PAL__DD2R*(DV*BF+V*DBF)/CJ;
497              
498             /* Parallax */
499             V=0.;
500             DV=0.;
501 6080 100         for (N=NP-1; N>=0; N--) { /* DO N=NP,1,-1 */
502 5890           COEFF=TP[N];
503 5890           EMN=(double)(ITP[N][0]);
504 5890           EMPN=(double)(ITP[N][1]);
505 5890           DN=(double)(ITP[N][2]);
506 5890           FN=(double)(ITP[N][3]);
507 5890           I=ITP[N][4];
508 5890 100         if (I == 0) {
509             EN=1.;
510             DEN=0.;
511 2470 100         } else if (I == 1) {
512             EN=E;
513             DEN=DE;
514             } else {
515             EN=ESQ;
516             DEN=DESQ;
517             }
518 5890           THETA=EMN*EM+EMPN*EMP+DN*D+FN*F;
519 5890           DTHETA=EMN*DEM+EMPN*DEMP+DN*DD+FN*DF;
520 5890           FTHETA=cos(THETA);
521 5890           V=V+COEFF*FTHETA*EN;
522 5890           DV=DV+COEFF*(-sin(THETA)*DTHETA*EN+FTHETA*DEN);
523             }
524 190           P=PAL__DD2R*V;
525 190           DP=PAL__DD2R*DV/CJ;
526              
527             /*
528             * Transformation into final form
529             */
530              
531             /* Parallax to distance (AU, AU/sec) */
532 190           SP=sin(P);
533 190           R=ERADAU/SP;
534 190           DR=-R*DP*cos(P)/SP;
535              
536             /* Longitude, latitude to x,y,z (AU) */
537 190           SEL=sin(EL);
538 190           CEL=cos(EL);
539 190           SB=sin(B);
540 190           CB=cos(B);
541 190           RCB=R*CB;
542 190           RBD=R*DB;
543 190           W=RBD*SB-CB*DR;
544 190           X=RCB*CEL;
545 190           Y=RCB*SEL;
546 190           Z=R*SB;
547 190           XD=-Y*DEL-W*CEL;
548 190           YD=X*DEL-W*SEL;
549 190           ZD=RBD*CB+SB*DR;
550              
551             /* Julian centuries since J2000 */
552 190           T=(date-51544.5)/36525.;
553              
554             /* Fricke equinox correction */
555 190           EPJ=2000.+T*100.;
556 190           EQCOR=PAL__DS2R*(0.035+0.00085*(EPJ-B1950));
557              
558             /* Mean obliquity (IAU 1976) */
559 190           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 190           SINEPS=sin(EPS);
563 190           COSEPS=cos(EPS);
564 190           ES=EQCOR*SINEPS;
565 190           EC=EQCOR*COSEPS;
566 190           pv[0]=X-EC*Y+ES*Z;
567 190           pv[1]=EQCOR*X+Y*COSEPS-Z*SINEPS;
568 190           pv[2]=Y*SINEPS+Z*COSEPS;
569 190           pv[3]=XD-EC*YD+ES*ZD;
570 190           pv[4]=EQCOR*XD+YD*COSEPS-ZD*SINEPS;
571 190           pv[5]=YD*SINEPS+ZD*COSEPS;
572              
573 190           }