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
|
|
|
|
|
|
} |