![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to moon.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [Research Unix] / researchv10no / cmd / sky|
Return to moon.c CVS log | Up to [Research Unix] / researchv10no / cmd / sky |
1.1 root 1: #include "sky.h"
2:
3: /*
4: * References:
5: * Brown,
6: * Improved Lunar Ephemeris
7: * Supplement to A.E. 1968
8: * Transformation corrections.
9: */
10:
11: double k1, k2, k3, k4;
12: double mnom, msun, noded, dmoon;
13: double sinx();
14: double cosx();
15:
16: extern struct moontab
17: {
18: float f;
19: char c[4];
20: } moontab[];
21:
22: extern struct moont1
23: {
24: float f1[2];
25: char c1[7];
26: } moont1[];
27:
28: moon()
29: {
30: register struct moontab *mp;
31: register struct moont1 *mp1;
32: double dlong, lsun, psun;
33: double eccm, eccs, chp, cpe;
34: double q0, v0, t0, m0, j0, sn0, l0;
35: double arg1, arg2, arg3, arg4, arg5, arg6, arg7;
36: double arg8, arg9, arg10, arg11, arg12, arg13;
37: double arg14, arg15, arg16, arg17;
38: double arg18, arg19, arg20, arg21, arg22;
39: double dgamma, k5, k6;
40: double lterms, sterms, cterms, nterms, pterms, spterms;
41: double gamma1, gamma2, gamma3, arglat;
42: double xmp, ymp, zmp;
43: double temp1, temp2;
44: double shsd;
45: double obl2;
46: double planp[7];
47:
48: object = "moon ";
49:
50: /*
51: * the fundamental elements - all referred to the epoch of
52: * Jan 0.5, 1900 and to the mean equinox of date.
53: */
54:
55: dlong = 270.434164 + 13.1763965268*eday - .001133*capt2
56: + 1.9e-6*capt3;
57: dlong -= .000086; /* empirical*/
58: argp = 334.329556 + .1114040803*eday - .010325*capt2
59: - 12.5e-6*capt3;
60: node = 259.183275 - .0529539222*eday + .002078*capt2
61: + 2.2e-6*capt3;
62: node += .000100; /* empirical */
63: lsun = 279.696678 + .9856473354*eday + .000303*capt2;
64: psun = 281.220844 + .0000470684*eday + .000453*capt2
65: + 3.3e-6*capt3;
66:
67: dlong = fmod(dlong, 360.);
68: argp = fmod(argp, 360.);
69: node = fmod(node, 360.);
70: lsun = fmod(lsun, 360.);
71: psun = fmod(psun, 360.);
72:
73: eccm = 22639.550;
74: eccs = .01675104 - .00004180*capt;
75: incl = 18461.400;
76: cpe = 124.986;
77: chp = 3422.451;
78:
79: /*
80: * some subsidiary elements - they are all longitudes
81: * and they are referred to the epoch 1/0.5 1900 and
82: * to the fixed mean equinox of 1850.0.
83: */
84:
85: q0 = 177.481153 + 4.0923388020*eday;
86: v0 = 342.069128 + 1.6021304820*eday;
87: t0 = 98.998753 + 0.9856091138*eday;
88: m0 = 293.049675 + 0.5240329445*eday;
89: j0 = 237.352319 + 0.0830912295*eday;
90: sn0 = 265.869508 + 0.03345974279*eday;
91: l0 = 269.736239 +13.1763583100*eday;
92:
93: /*
94: * the following are periodic corrections to the
95: * fundamental elements and constants.
96: * arg4 is the "Great Venus Inequality".
97: */
98:
99: arg1 = 41.1 + 20.2*(capt+.5);
100: arg2 = dlong - argp + 33. + 3.*t0 - 10.*v0 - 2.6*(capt+.5);
101: arg3 = dlong + 3.*argp - 4.*lsun + 67. - 23.*t0 + 25.*m0;
102: arg4 = dlong - argp + 151.1 + 16.*t0 - 18.*v0 - (capt+.5);
103: arg5 = node;
104: arg6 = node + 276.2 - 2.3*(capt+.5);
105: arg7 = 152. + 119.*(capt+0.5);
106: arg8 = dlong + argp - 2.*lsun + 105. + t0 - 3.*q0;
107: arg9 = 313.9 + 13.*t0 - 8.*v0;
108: arg10 = dlong - argp + 112.0 + 29.*t0 - 26.*v0;
109: arg11 = dlong - argp + 21.*t0 - 21.*v0;
110: arg12 = dlong + argp - 2.*lsun + 273. + 21.*t0 - 20.*v0;
111: arg13 = dlong + argp - 2.*lsun + 303. + 8.*t0 - 12.*v0;
112: arg14 = 2.*lsun - 2.*node + 270. + 6.*t0 - 5.*v0;
113: arg15 = dlong + 2.*lsun - 3.*argp + 24.*t0 - 24.*v0;
114: arg16 = dlong - lsun + 262. + 12.*t0 - 15.*v0;
115: arg17 = dlong - lsun + 190. + 25.*t0 - 23.*v0;
116: arg18 = 43. - 8.*t0 + 15.*m0;
117: arg19 = node - lsun + 165. + 2.*m0;
118: arg20 = node + 290.1 - 0.9*(capt+.5);
119: arg21 = 115. + 38.5*(capt+.5);
120: arg22 = 216.0 - 8.*t0 + 15.*m0;
121: arg1 *= radian;
122: arg2 *= radian;
123: arg3 *= radian;
124: arg4 *= radian;
125: arg5 *= radian;
126: arg6 *= radian;
127: arg7 *= radian;
128: arg8 *= radian;
129: arg9 *= radian;
130: arg10 *= radian;
131: arg11 *= radian;
132: arg12 *= radian;
133: arg13 *= radian;
134: arg14 *= radian;
135: arg15 *= radian;
136: arg16 *= radian;
137: arg17 *= radian;
138: arg18 *= radian;
139: arg19 *= radian;
140: arg20 *= radian;
141: arg21 *= radian;
142: arg22 *= radian;
143:
144: dlong += (
145: 0.84 *sin(arg1)
146: + 0.31 *sin(arg2)
147: + 0.04 *sin(arg3)
148: + 14.27 *sin(arg4)
149: + 7.261*sin(arg5)
150: + 0.282*sin(arg6)
151: + 0.04 *sin(arg7)
152: + 0.075*sin(arg8)
153: + 0.237*sin(arg9)
154: + 0.108*sin(arg10)
155: + 0.030*sin(arg11)
156: + 0.126*sin(arg12)
157: + 0.033*sin(arg13)
158: + 0.054*sin(arg14)
159: + 0.010*sin(arg15)
160: + 0.013*sin(arg16)
161: + 0.013*sin(arg17)
162: + 0.026*sin(arg18)
163: + 0.017*sin(arg19)
164: )/3600.;
165:
166: argp += (
167: - 2.10 *sin(arg1)
168: - 0.118*sin(arg4)
169: - 2.076*sin(arg5)
170: - 0.840*sin(arg6)
171: - 0.100*sin(arg7)
172: - 0.593*sin(arg9)
173: - 0.065*sin(arg18)
174: )/3600.;
175:
176: node += (
177: 0.63*sin(arg1)
178: + 0.17*sin(arg4)
179: + 95.96*sin(arg5)
180: + 15.58*sin(arg6)
181: + 1.86*sin(arg20)
182: )/3600.;
183:
184: t0 += (
185: -6.40*sin(arg1)
186: -0.27*sin(arg7)
187: -1.89*sin(arg9)
188: +0.20*sin(arg22)
189: )/3600.;
190:
191: lsun += (
192: -6.40*sin(arg1)
193: -0.27*sin(arg7)
194: -1.89*sin(arg9)
195: +0.20*sin(arg22)
196: )/3600.;
197:
198: dgamma = - 4.318*cos(arg5)
199: - 0.698*cos(arg6)
200: - 0.083*cos(arg20);
201:
202: j0 +=
203: 0.33*sin(arg21);
204:
205: sn0 +=
206: - 0.83*sin(arg21);
207:
208: /*
209: * the following factors account for the fact that the
210: * eccentricity, solar eccentricity, inclination and
211: * parallax used by Brown to make up his coefficients
212: * are both wrong and out of date. Brown did the same
213: * thing in a different way.
214: */
215:
216: k1 = eccm/22639.500;
217: k2 = eccs/.01675104;
218: k3 = 1. + 2.708e-6 + .000108008*dgamma;
219: k4 = cpe/125.154;
220: k5 = chp/3422.700;
221:
222: /*
223: * the principal arguments that are used to compute
224: * perturbations are the following differences of the
225: * fundamental elements.
226: */
227:
228: mnom = dlong - argp;
229: msun = lsun - psun;
230: noded = dlong - node;
231: dmoon = dlong - lsun;
232:
233: /*
234: * solar terms in longitude
235: */
236:
237: lterms = 0.0;
238: mp = moontab;
239: for(;;) {
240: if(mp->f == 0.0){
241: mp++;
242: break;
243: }
244: lterms += sinx(mp->f,
245: mp->c[0], mp->c[1],
246: mp->c[2], mp->c[3], 0.0);
247: mp++;
248: }
249:
250: lterms +=
251: (294.e-9*eday - 9193./1296000.*dgamma + .0013)*sinx(1.0,0,0,0,2,0.)
252: +(-220.e-9*eday + 5282./1296000.*dgamma)*sinx(1.0,1,0,0,-2,0.);
253:
254: planp[1] = q0;
255: planp[2] = v0;
256: planp[3] = t0;
257: planp[4] = m0;
258: planp[5] = j0;
259: planp[6] = sn0;
260:
261: /*
262: * planetary terms in longitude
263: */
264:
265: mp1 = moont1;
266: for(;;){
267: if(mp1->f1[0] == 0.){
268: mp1++;
269: break;
270: }
271: lterms += sinx(mp1->f1[0], mp1->c1[0], mp1->c1[1],
272: mp1->c1[2], mp1->c1[3],
273: mp1->c1[4]*t0+mp1->c1[5]*planp[mp1->c1[6]]+mp1->f1[1]);
274: mp1++;
275: }
276:
277: lterms += sinx(-.189, 0,0,0,0, node) /*IAU*/
278: + sinx(-.013, -1,0,0,0, node) /*IAU*/
279: + sinx(-.013, 1,0,0,0, node); /*IAU*/
280: lterms += sinx(0.019, 0,0,0,0, 5.*t0-3.*v0+node+216.0);
281: lterms += sinx(0.016, 0,0,0,0, -5.*t0+3.*v0+node+075.0);
282: lterms += sinx(-.038, 0,0,0,0, 2.*node);
283:
284: /*
285: * solar terms in latitude
286: */
287:
288: sterms = 0.0;
289: for(;;) {
290: if(mp->f == 0.0){
291: mp++;
292: break;
293: }
294: sterms += sinx(mp->f,
295: mp->c[0], mp->c[1],
296: mp->c[2], mp->c[3], 0.0);
297: mp++;
298: }
299:
300: cterms = 0.0;
301: for(;;) {
302: if(mp->f == 0.0){
303: mp++;
304: break;
305: }
306: cterms += cosx(mp->f,
307: mp->c[0], mp->c[1],
308: mp->c[2], mp->c[3], 0.0);
309: mp++;
310: }
311:
312: nterms = 0.0;
313: for(;;) {
314: if(mp->f == 0.0){
315: mp++;
316: break;
317: }
318: nterms += sinx(mp->f,
319: mp->c[0], mp->c[1],
320: mp->c[2], mp->c[3], 0.0);
321: mp++;
322: }
323:
324: /*
325: * planetary terms in latitude
326: */
327:
328: pterms = 0.;
329: for(;;){
330: if(mp1->f1[0] == 0.){
331: mp1++;
332: break;
333: }
334: pterms += sinx(mp1->f1[0], mp1->c1[0], mp1->c1[1],
335: mp1->c1[2], mp1->c1[3],
336: mp1->c1[4]*t0+mp1->c1[5]*planp[mp1->c1[6]]+mp1->f1[1]);
337: mp1++;
338: }
339:
340: pterms +=
341: sinx(0.014, 0,0,0,0, -2.*t0+2.*v0+l0+285.0)
342: + sinx(0.027, 0,0,0,0, -1.*t0+1.*v0+l0+285.0)
343: + sinx(0.015, 0,0,0,0, t0-v0+l0+105.0)
344: + sinx(0.077, 0,0,0,0, 5.*t0-3.*v0+l0+215.6)
345: + sinx(0.025, 0,0,0,0, -6.*t0+4.*v0+l0+255.0)
346: + sinx(0.074, 0,0,0,0, -5.*t0+3.*v0+l0+051.6)
347: + sinx(0.018, 0,0,0,0, -4.*t0+2.*v0+l0+075.0)
348: + sinx(0.010, 0,0,0,0, -3.*t0+v0+l0+075.0)
349: + sinx(0.030, 0,0,0,0, 8.*t0-5.*v0+l0+125.0);
350: pterms +=
351: sinx(0.010, 0,0,0,0, -t0+2.*m0+l0+345.0)
352: + sinx(0.035, 0,0,0,0, 2.*j0+l0+168.0)
353: + sinx(0.018, 0,0,0,0, -2.*j0+l0+024.0)
354: + sinx(-.017, 0,0,0,0, l0)
355: + sinx(0.083, 0,0,1,0, 2.*node)
356: + sinx(0.215, 0,0,0,0, dlong) /*IAU*/
357: + sinx(-.012, -1,0,0,0, dlong) /*IAU*/
358: + sinx(0.011, 1,0,0,0, dlong); /*IAU*/
359:
360: /*
361: * solar terms in parallax
362: */
363:
364: spterms = 3422.700;
365: for(;;) {
366: if(mp->f == 0.0){
367: mp++;
368: break;
369: }
370: spterms += cosx(mp->f,
371: mp->c[0], mp->c[1],
372: mp->c[2], mp->c[3], 0.0);
373: mp++;
374: }
375:
376: /*
377: * planetary terms in parallax
378: */
379:
380: for(;;){
381: if(mp1->f1[0] == 0.){
382: mp1++;
383: break;
384: }
385: spterms += cosx(mp1->f1[0], mp1->c1[0], mp1->c1[1],
386: mp1->c1[2], mp1->c1[3],
387: mp1->c1[4]*t0+mp1->c1[5]*planp[mp1->c1[6]]+mp1->f1[1]);
388: mp1++;
389: }
390:
391: /*
392: * computation of longitude
393: */
394:
395: lambda = (dlong + lterms/3600.)*radian;
396:
397: /*
398: * computation of latitude
399: */
400:
401: arglat = (noded + sterms/3600.)*radian;
402: gamma1 = 18519.700 * k3;
403: gamma2 = -6.241 * k3*k3*k3;
404: gamma3 = 0.004 * k3*k3*k3*k3*k3;
405:
406: k6 = (gamma1 + cterms) / gamma1;
407:
408: beta = k6 * (gamma1*sin(arglat) + gamma2*sin(3.*arglat)
409: + gamma3*sin(5.*arglat) + nterms)
410: + pterms;
411: if(flags & OCCULT)
412: beta -= 0.6;
413: beta *= radsec;
414:
415: /*
416: * computation of parallax
417: */
418:
419: spterms = k5 * spterms *radsec;
420: hp = spterms + (spterms*spterms*spterms)/6.;
421:
422: rad = hp/radsec;
423: georad = 1.;
424: semi = .0799 + .272453*(hp/radsec);
425: if(dmoon < 0.)
426: dmoon += 360.;
427: mag = dmoon/360.;
428:
429: /*
430: * change to equatorial coordinates
431: */
432:
433: lambda += psi;
434: obl2 = obliq + eps;
435: xmp = georad*cos(lambda)*cos(beta);
436: ymp = georad*(sin(lambda)*cos(beta)*cos(obl2) - sin(obl2)*sin(beta));
437: zmp = georad*(sin(lambda)*cos(beta)*sin(obl2) + cos(obl2)*sin(beta));
438:
439: alpha = atan2(ymp, xmp);
440: delta = atan2(zmp, sqrt(xmp*xmp+ymp*ymp));
441:
442: /*
443: *c lunar eclipse computation
444: */
445:
446: shsd = 1.0183*hp/radsec - 969.85/rps;
447: temp1 = sin(shdecl)*sin(delta) + cos(shdecl)*cos(delta)
448: *cos(shra - alpha);
449: temp2 = atan2(sqrt(1.-temp1*temp1),temp1)/radsec;
450: temp2 = semi + shsd - temp2;
451: temp2 = temp2/(2.*semi);
452: if(temp2 >= 0.){
453: if(temp2 < 1.)
454: printf("Partial Lunar Eclipse: Mag. = %.4f\n", temp2);
455: else
456: printf("Total Lunar Eclipse: Mag. = %.4f\n", temp2);
457: }
458:
459:
460: geolam = lambda;
461: geobet = beta;
462:
463: geo();
464:
465: /*
466: * solar eclipse computation
467: */
468:
469: if(!((flags&GEO)||(flags&HELIO))){
470: temp1 = sin(sundec)*sin(decl2) + cos(sundec)*cos(decl2)
471: *cos(sunra-ra);
472: temp1 = atan2(sqrt(1.-temp1*temp1),temp1)/radsec;
473: if(temp1 <= (semi2+sunsd)){
474: temp2 = (semi2+sunsd-temp1)/(2.*sunsd);
475: if(temp1 >= fabs(sunsd-semi2))
476: printf("Partial Solar Eclipse: Mag. = %.4f\n", temp2);
477: else if(temp2 > 1.)
478: printf("Total Solar Eclipse: Mag. = %.4f\n", temp2);
479: else
480: printf("Annular Solar Eclipse: Mag. = %.4f\n", temp2);
481: }
482: }
483:
484: /*
485: * constants for occultation computations
486: */
487:
488: moonra = ra;
489: moonde = decl2;
490: moonsd = semi2;
491:
492: }
493:
494: double
495: sinx(coef,i,j,k,m,angle)
496: double coef, angle;
497: {
498: double x;
499:
500: x = i*mnom + j*msun + k*noded + m*dmoon + angle;
501: x = coef*sin(x*radian);
502: if(i < 0)
503: i = -i;
504: for(; i>0; i--)
505: x *= k1;
506: if(j < 0)
507: j = -j;
508: for(; j>0; j--)
509: x *= k2;
510: if(k < 0)
511: k = -k;
512: for(; k>0; k--)
513: x *= k3;
514: if(m & 1)
515: x *= k4;
516:
517: return(x);
518: }
519:
520: double
521: cosx(coef,i,j,k,m,angle)
522: double coef, angle;
523: {
524: double x;
525:
526: x = i*mnom + j*msun + k*noded + m*dmoon + angle;
527: x = coef*cos(x*radian);
528: if(i < 0)
529: i = -i;
530: for(; i>0; i--)
531: x *= k1;
532: if(j < 0)
533: j = -j;
534: for(; j>0; j--)
535: x *= k2;
536: if(k < 0)
537: k = -k;
538: for(; k>0; k--)
539: x *= k3;
540: if(m & 1)
541: x *= k4;
542:
543: return(x);
544: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.