1 #include <cmath>
2 #include <cstdio>
3 #include <cstdlib>
4 using namespace std;
5
6 #include "elp82b.h"
7
8 #include "xpUtil.h"
9
10 void
moon(const double jd,double & X,double & Y,double & Z)11 moon(const double jd, double &X, double &Y, double &Z)
12 {
13 const double sec_to_rad = deg_to_rad / 3600;
14
15 double R[3] = { 0, 0, 0 };
16
17 // w[0] is the mean mean longitude of the moon
18 // w[1] is the mean longitude of the lunar perigee
19 // w[2] is the mean longitude of the lunar ascending node
20 // Units are seconds of arc
21 double w[3][5];
22 w[0][0] = (3600 * 218. + 60 * 18. + 59.95571);
23 w[1][0] = (3600 * 83. + 60 * 21. + 11.67475);
24 w[2][0] = (3600 * 125. + 60 * 2. + 40.39816);
25
26 w[0][1] = 1732559343.73604;
27 w[1][1] = 14643420.2632;
28 w[2][1] = -6967919.3622;
29
30 w[0][2] = -5.8883;
31 w[1][2] = -38.2776;
32 w[2][2] = 6.3622;
33
34 w[0][3] = 0.6604e-2;
35 w[1][3] = -0.45047e-1;
36 w[2][3] = 0.7625e-2;
37
38 w[0][4] = -0.3169e-4;
39 w[1][4] = 0.21301e-3;
40 w[2][4] = -0.3586e-4;
41
42 // T and omega are angles of the inertial mean ecliptic of J2000
43 // referred to the inertial mean equinox of J2000
44 double T[5];
45 T[0] = (3600 * 100. + 60 * 27. + 59.22059);
46 T[1] = 129597742.2758;
47 T[2] = -0.0202;
48 T[3] = 0.9e-5;
49 T[4] = 0.15e-6;
50
51 double omega[5];
52 omega[0] = (3600 * 102. + 60 * 56. + 14.42753);
53 omega[1] = 1161.2283;
54 omega[2] = 0.5327;
55 omega[3] = -0.138e-3;
56 omega[4] = 0;
57
58 // precession between J2000 and the date
59 double preces[5];
60 preces[0] = 0;
61 preces[1] = 5029.0966;
62 preces[2] = 1.1120;
63 preces[3] = .77e-4;
64 preces[4] = -.2353e-4;
65
66 // planetary longitudes in J2000
67 double plon[8][2];
68 plon[0][0] = (3600 * 252. + 60 * 15. + 3.25986);
69 plon[1][0] = (3600 * 181. + 60 * 58. + 47.28305);
70 plon[2][0] = T[0];
71 plon[3][0] = (3600 * 355. + 60 * 25. + 59.78866);
72 plon[4][0] = (3600 * 34. + 60 * 21. + 5.34212);
73 plon[5][0] = (3600 * 50. + 60 * 4. + 38.89694);
74 plon[6][0] = (3600 * 314. + 60 * 3. + 18.01841);
75 plon[7][0] = (3600 * 304. + 60 * 20. + 55.19575);
76
77 plon[0][1] = 538101628.68898;
78 plon[1][1] = 210664136.43355;
79 plon[2][1] = T[1];
80 plon[3][1] = 68905077.59284;
81 plon[4][1] = 10925660.42861;
82 plon[5][1] = 4399609.65932;
83 plon[6][1] = 1542481.19393;
84 plon[7][1] = 786550.32074;
85
86 // Convert to radians
87 for (int i = 0; i < 8; i++)
88 {
89 plon[i][0] *= sec_to_rad;
90 plon[i][1] *= sec_to_rad;
91 }
92
93 // Corrections of the constants (fit to DE200/LE200)
94 double delnu = 0.55604 / w[0][1];
95 double dele = 0.01789 * sec_to_rad;
96 double delg = -0.08066 * sec_to_rad;
97 double delnp = -0.06424 / w[0][1];
98 double delep = -0.12879 * sec_to_rad;
99
100 // Delaunay's arguments
101 double del[4][5];
102 for (int i = 0; i < 5; i++)
103 {
104 del[0][i] = w[0][i] - T[i];
105 del[1][i] = T[i] - omega[i];
106 del[2][i] = w[0][i] - w[1][i];
107 del[3][i] = w[0][i] - w[2][i];
108 for (int j = 0; j < 4; j++)
109 del[j][i] *= sec_to_rad;
110 }
111 del[0][0] += M_PI;
112
113 double zeta[2];
114 zeta[0] = w[0][0] * sec_to_rad;
115 zeta[1] = (w[0][1] + preces[1]) * sec_to_rad;
116
117 // precession matrix
118 double p[5], q[5];
119 p[0] = 0.10180391e-4;
120 p[1] = 0.47020439e-6;
121 p[2] = -0.5417367e-9;
122 p[3] = -0.2507948e-11;
123 p[4] = 0.463486e-14;
124 q[0] = -0.113469002e-3;
125 q[1] = 0.12372674e-6;
126 q[2] = 0.1265417e-8;
127 q[3] = -0.1371808e-11;
128 q[4] = -0.320334e-14;
129
130 double t[5];
131 t[0] = 1;
132 t[1] = (jd - 2451545.0) / 36525;
133 t[2] = t[1] * t[1];
134 t[3] = t[2] * t[1];
135 t[4] = t[2] * t[2];
136
137 double ath = 384747.9806743165;
138 double a0 = 384747.9806448954;
139 double am = 0.074801329518;
140 double alfa = 0.002571881335;
141 double dtasm = 2 * alfa / (3 * am);
142
143 // Main problem
144 for (int i = 1; i < 4; i++)
145 {
146 int iv = ((i - 1) % 3);
147
148 const int *ilu = ILU[i-1];
149 const double *p_coef = COEF[i-1];
150
151 for (int ii = 0; ii < NUM[i-1]; ii++)
152 {
153 if (ii > 0)
154 {
155 ilu += 4;
156 p_coef += 7;
157 }
158
159 double coef[7];
160 for (int j = 0; j < 7; j++)
161 coef[j] = p_coef[j];
162
163 double x = coef[0];
164
165 double tgv = coef[1] + dtasm * coef[5];
166
167 if (i == 3)
168 coef[0] -= 2 * coef[0] * delnu / 3;
169
170 x = (coef[0] + tgv * (delnp - am * delnu)
171 + coef[2] * delg + coef[3] * dele
172 + coef[4] * delep);
173 double y = 0;
174
175 for (int k = 0; k < 5; k++)
176 {
177 for (int j = 0; j < 4; j++)
178 {
179 y += ilu[j] * del[j][k] * t[k];
180 }
181 }
182
183 if (i == 3) y += M_PI_2;
184 y = fmod(y, 2 * M_PI);
185
186 R[iv] += x * sin(y);
187
188 }
189 }
190
191 // Figures - Tides - Relativity - Solar eccentricity
192 for (int i = 4; i < 37; i++)
193 {
194 if (i > 9 && i < 22) continue;
195 int iv = ((i - 1) % 3);
196
197 const int *ilu = ILU[i-1];
198 const int *p_iz = IZ[i-1];
199 const double *p_pha = PHA[i-1];
200 const double *p_x = XX[i-1];
201
202 for (int ii = 0; ii < NUM[i-1]; ii++)
203 {
204 if (ii > 0)
205 {
206 ilu += 4;
207 p_iz++;
208 p_pha++;
209 p_x++;
210 }
211
212 int iz = *p_iz;
213 double pha = *p_pha;
214 double x = *p_x;
215
216 if ((i > 6 && i < 10) || (i > 24 && i < 28)) x *= t[1];
217 if (i > 33 && i < 37) x *= t[2];
218
219 double y = pha * deg_to_rad;
220
221 for (int k = 0; k < 2; k++)
222 {
223 y += iz * zeta[k] * t[k];
224 for (int l = 0; l < 4; l++)
225 {
226 y += ilu[l] * del[l][k] * t[k];
227 }
228 }
229 y = fmod(y, 2*M_PI);
230 R[iv] += x * sin(y);
231 }
232 }
233
234 // Planetary perturbations
235 for (int i = 10; i < 22; i++)
236 {
237 int iv = ((i - 1) % 3);
238
239 const int *ipla = IPLA[i-1];
240 const double *p_pha = PHA[i-1];
241 const double *p_x = XX[i-1];
242
243 for (int ii = 0; ii < NUM[i-1]; ii++)
244 {
245 if (ii > 0)
246 {
247 ipla += 11;
248 p_pha++;
249 p_x++;
250 }
251
252 double pha = *p_pha;
253 double x = *p_x;
254
255 if ((i > 12 && i < 16) || (i > 18 && i < 22)) x *= t[1];
256
257 double y = pha * deg_to_rad;
258
259 if (i < 16)
260 {
261 for (int k = 0; k < 2; k++)
262 {
263 y += (ipla[8] * del[0][k] + ipla[9] * del[2][k]
264 + ipla[10] * del[3][k]) * t[k];
265 for (int l = 0; l < 8; l++)
266 {
267 y += ipla[l] * plon[l][k] * t[k];
268 }
269 }
270 }
271 else
272 {
273 for (int k = 0; k < 2; k++)
274 {
275 for (int l = 0; l < 4; l++)
276 y += ipla[l+7] * del[l][k] * t[k];
277 for (int l = 0; l < 7; l++)
278 {
279 y += ipla[l] * plon[l][k] * t[k];
280 }
281 }
282 }
283
284 y = fmod(y, 2*M_PI);
285 R[iv] += x * sin(y);
286 }
287 }
288
289 // Change of coordinates
290 for (int i = 0; i < 5; i++)
291 R[0] += w[0][i] * t[i];
292 R[0] *= sec_to_rad;
293 R[1] *= sec_to_rad;
294 R[2] *= a0/ath;
295
296 double x1 = R[2] * cos(R[1]);
297 double x2 = x1 * sin(R[0]);
298 x1 *= cos(R[0]);
299 double x3 = R[2] * sin(R[1]);
300
301 double pw = 0;
302 double qw = 0;
303 for (int i = 0; i < 5; i++)
304 {
305 pw += p[i] * t[i];
306 qw += q[i] * t[i];
307 }
308
309 pw *= t[1];
310 qw *= t[1];
311
312 double ra = 2 * sqrt(1 - pw * pw - qw * qw);
313 double pwqw = 2 * pw * qw;
314 double pw2 = 1 - 2 * pw * pw;
315 double qw2 = 1 - 2 * qw * qw;
316 pw *= ra;
317 qw *= ra;
318
319 R[0] = pw2*x1+pwqw*x2+pw*x3;
320 R[1] = pwqw*x1+qw2*x2-qw*x3;
321 R[2] = -pw*x1+qw*x2+(pw2+qw2-1)*x3;
322
323 X = R[0] / AU_to_km;
324 Y = R[1] / AU_to_km;
325 Z = R[2] / AU_to_km;
326
327 // rotate to earth equator J2000
328 const double eps = 23.4392911 * deg_to_rad;
329 rotateX(X, Y, Z, -eps);
330 }
331