1 #include "astro.h"
2
3 double k1, k2, k3, k4;
4 double mnom, msun, noded, dmoon;
5
6 void
moon(void)7 moon(void)
8 {
9 Moontab *mp;
10 double dlong, lsun, psun;
11 double eccm, eccs, chp, cpe;
12 double v0, t0, m0, j0;
13 double arg1, arg2, arg3, arg4, arg5, arg6, arg7;
14 double arg8, arg9, arg10;
15 double dgamma, k5, k6;
16 double lterms, sterms, cterms, nterms, pterms, spterms;
17 double gamma1, gamma2, gamma3, arglat;
18 double xmp, ymp, zmp;
19 double obl2;
20
21 /*
22 * the fundamental elements - all referred to the epoch of
23 * Jan 0.5, 1900 and to the mean equinox of date.
24 */
25
26 dlong = 270.434164 + 13.1763965268*eday - .001133*capt2
27 + 2.e-6*capt3;
28 argp = 334.329556 + .1114040803*eday - .010325*capt2
29 - 12.e-6*capt3;
30 node = 259.183275 - .0529539222*eday + .002078*capt2
31 + 2.e-6*capt3;
32 lsun = 279.696678 + .9856473354*eday + .000303*capt2;
33 psun = 281.220833 + .0000470684*eday + .000453*capt2
34 + 3.e-6*capt3;
35
36 dlong = fmod(dlong, 360.);
37 argp = fmod(argp, 360.);
38 node = fmod(node, 360.);
39 lsun = fmod(lsun, 360.);
40 psun = fmod(psun, 360.);
41
42 eccm = 22639.550;
43 eccs = .01675104 - .00004180*capt;
44 incl = 18461.400;
45 cpe = 124.986;
46 chp = 3422.451;
47
48 /*
49 * some subsidiary elements - they are all longitudes
50 * and they are referred to the epoch 1/0.5 1900 and
51 * to the fixed mean equinox of 1850.0.
52 */
53
54 v0 = 342.069128 + 1.6021304820*eday;
55 t0 = 98.998753 + 0.9856091138*eday;
56 m0 = 293.049675 + 0.5240329445*eday;
57 j0 = 237.352319 + 0.0830912295*eday;
58
59 /*
60 * the following are periodic corrections to the
61 * fundamental elements and constants.
62 * arg3 is the "Great Venus Inequality".
63 */
64
65 arg1 = 41.1 + 20.2*(capt+.5);
66 arg2 = dlong - argp + 33. + 3.*t0 - 10.*v0 - 2.6*(capt+.5);
67 arg3 = dlong - argp + 151.1 + 16.*t0 - 18.*v0 - (capt+.5);
68 arg4 = node;
69 arg5 = node + 276.2 - 2.3*(capt+.5);
70 arg6 = 313.9 + 13.*t0 - 8.*v0;
71 arg7 = dlong - argp + 112.0 + 29.*t0 - 26.*v0;
72 arg8 = dlong + argp - 2.*lsun + 273. + 21.*t0 - 20.*v0;
73 arg9 = node + 290.1 - 0.9*(capt+.5);
74 arg10 = 115. + 38.5*(capt+.5);
75 arg1 *= radian;
76 arg2 *= radian;
77 arg3 *= radian;
78 arg4 *= radian;
79 arg5 *= radian;
80 arg6 *= radian;
81 arg7 *= radian;
82 arg8 *= radian;
83 arg9 *= radian;
84 arg10 *= radian;
85
86 dlong +=
87 (0.84 *sin(arg1)
88 + 0.31 *sin(arg2)
89 + 14.27 *sin(arg3)
90 + 7.261*sin(arg4)
91 + 0.282*sin(arg5)
92 + 0.237*sin(arg6)
93 + 0.108*sin(arg7)
94 + 0.126*sin(arg8))/3600.;
95
96 argp +=
97 (- 2.10 *sin(arg1)
98 - 0.118*sin(arg3)
99 - 2.076*sin(arg4)
100 - 0.840*sin(arg5)
101 - 0.593*sin(arg6))/3600.;
102
103 node +=
104 (0.63*sin(arg1)
105 + 0.17*sin(arg3)
106 + 95.96*sin(arg4)
107 + 15.58*sin(arg5)
108 + 1.86*sin(arg9))/3600.;
109
110 t0 +=
111 (- 6.40*sin(arg1)
112 - 1.89*sin(arg6))/3600.;
113
114 psun +=
115 (6.40*sin(arg1)
116 + 1.89*sin(arg6))/3600.;
117
118 dgamma = - 4.318*cos(arg4)
119 - 0.698*cos(arg5)
120 - 0.083*cos(arg9);
121
122 j0 +=
123 0.33*sin(arg10);
124
125 /*
126 * the following factors account for the fact that the
127 * eccentricity, solar eccentricity, inclination and
128 * parallax used by Brown to make up his coefficients
129 * are both wrong and out of date. Brown did the same
130 * thing in a different way.
131 */
132
133 k1 = eccm/22639.500;
134 k2 = eccs/.01675104;
135 k3 = 1. + 2.708e-6 + .000108008*dgamma;
136 k4 = cpe/125.154;
137 k5 = chp/3422.700;
138
139 /*
140 * the principal arguments that are used to compute
141 * perturbations are the following differences of the
142 * fundamental elements.
143 */
144
145 mnom = dlong - argp;
146 msun = lsun - psun;
147 noded = dlong - node;
148 dmoon = dlong - lsun;
149
150 /*
151 * solar terms in longitude
152 */
153
154 lterms = 0.0;
155 mp = moontab;
156 for(;;) {
157 if(mp->f == 0.0)
158 break;
159 lterms += sinx(mp->f,
160 mp->c[0], mp->c[1],
161 mp->c[2], mp->c[3], 0.0);
162 mp++;
163 }
164 mp++;
165
166 /*
167 * planetary terms in longitude
168 */
169
170 lterms += sinx(0.822, 0,0,0,0, t0-v0);
171 lterms += sinx(0.307, 0,0,0,0, 2.*t0-2.*v0+179.8);
172 lterms += sinx(0.348, 0,0,0,0, 3.*t0-2.*v0+272.9);
173 lterms += sinx(0.176, 0,0,0,0, 4.*t0-3.*v0+271.7);
174 lterms += sinx(0.092, 0,0,0,0, 5.*t0-3.*v0+199.);
175 lterms += sinx(0.129, 1,0,0,0, -t0+v0+180.);
176 lterms += sinx(0.152, 1,0,0,0, t0-v0);
177 lterms += sinx(0.127, 1,0,0,0, 3.*t0-3.*v0+180.);
178 lterms += sinx(0.099, 0,0,0,2, t0-v0);
179 lterms += sinx(0.136, 0,0,0,2, 2.*t0-2.*v0+179.5);
180 lterms += sinx(0.083, -1,0,0,2, -4.*t0+4.*v0+180.);
181 lterms += sinx(0.662, -1,0,0,2, -3.*t0+3.*v0+180.0);
182 lterms += sinx(0.137, -1,0,0,2, -2.*t0+2.*v0);
183 lterms += sinx(0.133, -1,0,0,2, t0-v0);
184 lterms += sinx(0.157, -1,0,0,2, 2.*t0-2.*v0+179.6);
185 lterms += sinx(0.079, -1,0,0,2, -8.*t0+6.*v0+162.6);
186 lterms += sinx(0.073, 2,0,0,-2, 3.*t0-3.*v0+180.);
187 lterms += sinx(0.643, 0,0,0,0, -t0+j0+178.8);
188 lterms += sinx(0.187, 0,0,0,0, -2.*t0+2.*j0+359.6);
189 lterms += sinx(0.087, 0,0,0,0, j0+289.9);
190 lterms += sinx(0.165, 0,0,0,0, -t0+2.*j0+241.5);
191 lterms += sinx(0.144, 1,0,0,0, t0-j0+1.0);
192 lterms += sinx(0.158, 1,0,0,0, -t0+j0+179.0);
193 lterms += sinx(0.190, 1,0,0,0, -2.*t0+2.*j0+180.0);
194 lterms += sinx(0.096, 1,0,0,0, -2.*t0+3.*j0+352.5);
195 lterms += sinx(0.070, 0,0,0,2, 2.*t0-2.*j0+180.);
196 lterms += sinx(0.167, 0,0,0,2, -t0+j0+178.5);
197 lterms += sinx(0.085, 0,0,0,2, -2.*t0+2.*j0+359.2);
198 lterms += sinx(1.137, -1,0,0,2, 2.*t0-2.*j0+180.3);
199 lterms += sinx(0.211, -1,0,0,2, -t0+j0+178.4);
200 lterms += sinx(0.089, -1,0,0,2, -2.*t0+2.*j0+359.2);
201 lterms += sinx(0.436, -1,0,0,2, 2.*t0-3.*j0+7.5);
202 lterms += sinx(0.240, 2,0,0,-2, -2.*t0+2.*j0+179.9);
203 lterms += sinx(0.284, 2,0,0,-2, -2.*t0+3.*j0+172.5);
204 lterms += sinx(0.195, 0,0,0,0, -2.*t0+2.*m0+180.2);
205 lterms += sinx(0.327, 0,0,0,0, -t0+2.*m0+224.4);
206 lterms += sinx(0.093, 0,0,0,0, -2.*t0+4.*m0+244.8);
207 lterms += sinx(0.073, 1,0,0,0, -t0+2.*m0+223.3);
208 lterms += sinx(0.074, 1,0,0,0, t0-2.*m0+306.3);
209 lterms += sinx(0.189, 0,0,0,0, node+180.);
210
211 /*
212 * solar terms in latitude
213 */
214
215 sterms = 0;
216 for(;;) {
217 if(mp->f == 0)
218 break;
219 sterms += sinx(mp->f,
220 mp->c[0], mp->c[1],
221 mp->c[2], mp->c[3], 0);
222 mp++;
223 }
224 mp++;
225
226 cterms = 0;
227 for(;;) {
228 if(mp->f == 0)
229 break;
230 cterms += cosx(mp->f,
231 mp->c[0], mp->c[1],
232 mp->c[2], mp->c[3], 0);
233 mp++;
234 }
235 mp++;
236
237 nterms = 0;
238 for(;;) {
239 if(mp->f == 0)
240 break;
241 nterms += sinx(mp->f,
242 mp->c[0], mp->c[1],
243 mp->c[2], mp->c[3], 0);
244 mp++;
245 }
246 mp++;
247
248 /*
249 * planetary terms in latitude
250 */
251
252 pterms =
253 sinx(0.215, 0,0,0,0, dlong);
254
255 /*
256 * solar terms in parallax
257 */
258
259 spterms = 3422.700;
260 for(;;) {
261 if(mp->f == 0)
262 break;
263 spterms += cosx(mp->f,
264 mp->c[0], mp->c[1],
265 mp->c[2], mp->c[3], 0);
266 mp++;
267 }
268
269 /*
270 * planetary terms in parallax
271 */
272
273 //spterms = spterms;
274
275 /*
276 * computation of longitude
277 */
278
279 lambda = (dlong + lterms/3600.)*radian;
280
281 /*
282 * computation of latitude
283 */
284
285 arglat = (noded + sterms/3600.)*radian;
286 gamma1 = 18519.700 * k3;
287 gamma2 = -6.241 * k3*k3*k3;
288 gamma3 = 0.004 * k3*k3*k3*k3*k3;
289
290 k6 = (gamma1 + cterms) / gamma1;
291
292 beta = k6 * (gamma1*sin(arglat) + gamma2*sin(3.*arglat)
293 + gamma3*sin(5.*arglat) + nterms)
294 + pterms;
295 if(flags['o'])
296 beta -= 0.6;
297 beta *= radsec;
298
299 /*
300 * computation of parallax
301 */
302
303 spterms = k5 * spterms *radsec;
304 hp = spterms + (spterms*spterms*spterms)/6.;
305
306 rad = hp/radsec;
307 rp = 1.;
308 semi = .0799 + .272453*(hp/radsec);
309 if(dmoon < 0.)
310 dmoon += 360.;
311 mag = dmoon/360.;
312
313 /*
314 * change to equatorial coordinates
315 */
316
317 lambda += phi;
318 obl2 = obliq + eps;
319 xmp = rp*cos(lambda)*cos(beta);
320 ymp = rp*(sin(lambda)*cos(beta)*cos(obl2) - sin(obl2)*sin(beta));
321 zmp = rp*(sin(lambda)*cos(beta)*sin(obl2) + cos(obl2)*sin(beta));
322
323 alpha = atan2(ymp, xmp);
324 delta = atan2(zmp, sqrt(xmp*xmp+ymp*ymp));
325 meday = eday;
326 mhp = hp;
327
328 geo();
329 }
330
331 double
sinx(double coef,int i,int j,int k,int m,double angle)332 sinx(double coef, int i, int j, int k, int m, double angle)
333 {
334 double x;
335
336 x = i*mnom + j*msun + k*noded + m*dmoon + angle;
337 x = coef*sin(x*radian);
338 if(i < 0)
339 i = -i;
340 for(; i>0; i--)
341 x *= k1;
342 if(j < 0)
343 j = -j;
344 for(; j>0; j--)
345 x *= k2;
346 if(k < 0)
347 k = -k;
348 for(; k>0; k--)
349 x *= k3;
350 if(m & 1)
351 x *= k4;
352
353 return x;
354 }
355
356 double
cosx(double coef,int i,int j,int k,int m,double angle)357 cosx(double coef, int i, int j, int k, int m, double angle)
358 {
359 double x;
360
361 x = i*mnom + j*msun + k*noded + m*dmoon + angle;
362 x = coef*cos(x*radian);
363 if(i < 0)
364 i = -i;
365 for(; i>0; i--)
366 x *= k1;
367 if(j < 0)
368 j = -j;
369 for(; j>0; j--)
370 x *= k2;
371 if(k < 0)
372 k = -k;
373 for(; k>0; k--)
374 x *= k3;
375 if(m & 1)
376 x *= k4;
377
378 return x;
379 }
380