1 /*
2
3 A Moon for OpenWindows
4
5 Release 3.0
6
7 Designed and implemented by John Walker in December 1987,
8 Revised and updated in February of 1988.
9 Revised and updated again in June of 1988 by Ron Hitchens.
10 Revised and updated yet again in July/August of 1989 by Ron Hitchens.
11 Converted to OpenWindows in December of 1991 by John Walker.
12
13 This program is an OpenWindows tool which displays, as the icon
14 for a closed window, the current phase of the Moon. A subtitle in
15 the icon gives the age of the Moon in days and hours. If called
16 with the "-t" switch, it rapidly increments forward through time
17 to display the cycle of phases.
18
19 If you open the window, additional information is displayed
20 regarding the Moon. The information is generally accurate to
21 within ten minutes.
22
23 The algorithms used in this program to calculate the positions Sun
24 and Moon as seen from the Earth are given in the book "Practical
25 Astronomy With Your Calculator" by Peter Duffett-Smith, Second
26 Edition, Cambridge University Press, 1981. Ignore the word
27 "Calculator" in the title; this is an essential reference if
28 you're interested in developing software which calculates
29 planetary positions, orbits, eclipses, and the like. If you're
30 interested in pursuing such programming, you should also obtain:
31
32 "Astronomical Formulae for Calculators" by Jean Meeus, Third
33 Edition, Willmann-Bell, 1985. A must-have.
34
35 "Planetary Programs and Tables from -4000 to +2800" by Pierre
36 Bretagnon and Jean-Louis Simon, Willmann-Bell, 1986. If you want
37 the utmost (outside of JPL) accuracy for the planets, it's here.
38
39 "Celestial BASIC" by Eric Burgess, Revised Edition, Sybex, 1985.
40 Very cookbook oriented, and many of the algorithms are hard to dig
41 out of the turgid BASIC code, but you'll probably want it anyway.
42
43 Many of these references can be obtained from Willmann-Bell, P.O.
44 Box 35025, Richmond, VA 23235, USA. Phone: (804) 320-7016. In
45 addition to their own publications, they stock most of the
46 standard references for mathematical and positional astronomy.
47
48 This program was written by:
49
50 John Walker
51 Autodesk Neuchbtel
52 Avenue des Champs-Montants 14b
53 CH-2074 MARIN
54 Switzerland
55 Usenet: kelvin@Autodesk.com
56 Fax: 038/33 88 15
57 Voice: 038/33 76 33
58
59 This program is in the public domain: "Do what thou wilt shall be
60 the whole of the law". I'd appreciate receiving any bug fixes
61 and/or enhancements, which I'll incorporate in future versions of
62 the program. Please leave the original attribution information
63 intact so that credit and blame may be properly apportioned.
64
65 History:
66 --------
67 June 1988 Version 2.0 posted to usenet by John Walker
68
69 June 1988 Modified by Ron Hitchens
70 ronbo@vixen.uucp
71 ...!uunet!cs.utah.edu!caeco!vixen!ronbo
72 hitchens@cs.utexas.edu
73 to produce version 2.1.
74 Modified icon generation to show surface texture
75 on visible moon face.
76 Added a menu to allow switching in and out of
77 test mode, for entertainment value mostly.
78 Modified layout of information in open window display
79 to reduce the amount of pixels modified in each
80 update.
81
82 July 1989 Modified further for color displays. On a color Sun,
83 four colors will be used for the canvas and the icon.
84 Rather than just show the illuminated portion of
85 the moon, a color icon will also show the darkened
86 portion in a dark blue shade. The text on the icon
87 will also be drawn in a nice "buff" color, since there
88 was one more color left to use.
89 Add two command line args, "-c" and "-m" to explicitly
90 specify color or monochrome mode.
91 Use getopt to parse the args.
92 Change the tool menu slightly to use only one item
93 for switching in and out of test mode.
94
95 July 1989 Modified a little bit more a few days later to use 8
96 colors and an accurate grey-scale moon face created
97 by Joe Hitchens on an Amiga.
98 Added The Apollo 11 Commemorative Red Dot, to show
99 where Neil and Buzz went on vacation a few years ago.
100 Updated man page.
101
102 August 1989 Received version 2.3 of John Walker's original code.
103 Rolled in bug fixes to astronomical algorithms:
104
105 2.1 6/16/88 Bug fix. Table of phases didn't update
106 at the moment of the new moon. Call on
107 phasehunt didn't convert civil Julian
108 date to astronomical Julian date.
109 Reported by Dag Bruck
110 (dag@control.lth.se).
111
112 2.2 N/A (Superseded by new moon icon.)
113
114 2.3 6/7/89 Bug fix. Table of phases skipped the
115 phases for July 1989. This occurred
116 due to sloppy maintenance of the
117 synodic month index in the interchange
118 of information between phasehunt() and
119 meanphase(). I simplified and
120 corrected the handling of the month
121 index as phasehunt() steps along and
122 removed unneeded code from meanphase().
123 Reported by Bill Randle of Tektronix.
124 (billr@saab.CNA.TEK.COM).
125
126 January 1990 Ported to OpenWindows by John Walker.
127 All of Ron Hitchens' additions which were not
128 Sun-specific were carried on into the
129 OpenWindows version.
130
131 September 1993 reported to Motif (as God intended) by Cary Sandvig.
132 Some window reformatting was done as I was unable to
133 view the existing setup.
134
135 */
136
137 #include <stdio.h>
138 #include <stdlib.h>
139 #include <unistd.h>
140 #include <time.h>
141 #include <X11/Intrinsic.h>
142 #include <X11/StringDefs.h>
143 #include <X11/Shell.h>
144 #include <Xm/MainW.h>
145 #include <Xm/DrawingA.h>
146 #include <Xm/RowColumn.h>
147 #include <Xm/ToggleB.h>
148 #include <Xm/Label.h>
149 #include <Xm/Frame.h>
150
151 #define PI 3.14159265358979323846 /* Assume not near black hole nor in
152 Tennessee */
153
154 #include <math.h>
155
156 #define MOON_TITLE "XMoontool by John Walker v3.0.3"
157
158 #define UNKNOWN -2
159
160 #include "moon_icon"
161 #include "color.pix"
162
163 static char moon_icon_work[sizeof moon_icon_bits];
164 static char moon_icon_last[sizeof moon_icon_bits] = "abcd";
165 static char icon_label_last[20] = "";
166
167 #define CLOSED_SECS 120 /* update interval when tool closed */
168 #define OPEN_SECS 1 /* update interval when open */
169
170 /* 0 1 2 3 4 5 6 7*/
171 #define COLOR_R {0x00, 0xe0e0, 0xc0c0, 0x9090, 0xffff, 0x1010, 0x1010, 0x1010}
172 #define COLOR_G {0x00, 0xf0f0, 0xd0d0, 0xa0a0, 0x0000, 0x1010, 0x1010, 0x1010}
173 #define COLOR_B {0x00, 0xffff, 0xe0e0, 0xb0b0, 0x0000, 0x9090, 0x7070, 0x5050}
174
175 static int c_r[] = COLOR_R;
176 static int c_g[] = COLOR_G;
177 static int c_b[] = COLOR_B;
178
179 /* Astronomical constants */
180 #define epoch 2444238.5 /* 1980 January 0.0 */
181
182 /* Constants defining the Sun's apparent orbit */
183 #define elonge 278.833540 /* Ecliptic longitude of the Sun
184 at epoch 1980.0 */
185 #define elongp 282.596403 /* Ecliptic longitude of the Sun at
186 perigee */
187 #define eccent 0.016718 /* Eccentricity of Earth's orbit */
188 #define sunsmax 1.495985e8 /* Semi-major axis of Earth's orbit, km */
189 #define sunangsiz 0.533128 /* Sun's angular size, degrees, at
190 semi-major axis distance */
191
192 /* Elements of the Moon's orbit, epoch 1980.0 */
193 #define mmlong 64.975464 /* Moon's mean lonigitude at the epoch */
194 #define mmlongp 349.383063 /* Mean longitude of the perigee at the
195 epoch */
196 #define mlnode 151.950429 /* Mean longitude of the node at the
197 epoch */
198 #define minc 5.145396 /* Inclination of the Moon's orbit */
199 #define mecc 0.054900 /* Eccentricity of the Moon's orbit */
200 #define mangsiz 0.5181 /* Moon's angular size at distance a
201 from Earth */
202 #define msmax 384401.0 /* Semi-major axis of Moon's orbit in km */
203 #define mparallax 0.9507 /* Parallax at distance a from Earth */
204 #define synmonth 29.53058868 /* Synodic month (new Moon to new Moon) */
205 #define lunatbase 2423436.0 /* Base date for E. W. Brown's numbered
206 series of lunations (1923 January 16) */
207
208 /* Properties of the Earth */
209 #define earthrad 6378.16 /* Radius of Earth in kilometres */
210
211 /* Distance units */
212 #define MILE 0.621371 /* Conversion factor km -> mi */
213 #define units(x) ((x) == MILE) ? "mi" : "km"
214
215 /* Handy mathematical functions */
216 #define sgn(x) (((x) < 0) ? -1 : ((x) > 0 ? 1 : 0)) /* Extract sign */
217 #define abs(x) ((x) < 0 ? (-(x)) : (x)) /* Absolute val */
218 #define fixangle(a) ((a) - 360.0 * (floor((a) / 360.0))) /* Fix angle */
219 #define torad(d) ((d) * (PI / 180.0)) /* Deg->Rad */
220 #define todeg(d) ((d) * (180.0 / PI)) /* Rad->Deg */
221 #define dsin(x) (sin(torad((x)))) /* Sin from deg */
222 #define dcos(x) (cos(torad((x)))) /* Cos from deg */
223
224 static double metric = 1.0; /* Length unit factor */
225 static int testmode = FALSE; /* Rapid warp through time for debugging */
226 static double faketime = 0.0; /* Time increment for test mode */
227 static int color_mode = UNKNOWN; /* indicates color/mono mode */
228 static Boolean iconic = False; /* is the window iconified? */
229 static double nptime = 0.0; /* Next new moon time */
230 static double pptime = 0.0; /* Previous new moon time */
231 static int col_vals[16]; /* pixel values for color mode */
232
233 static char *moname[] = {
234 "January", "February", "March", "April", "May",
235 "June", "July", "August", "September",
236 "October", "November", "December"
237 };
238
239 static char *labels[] = {
240 "Julian date: ",
241 "Universal time: ",
242 "Local time: ",
243 "",
244 "Age of moon: ",
245 "Moon phase: ",
246 "Moon's distance:",
247 "Moon subtends: ",
248 "",
249 "Sun's distance: ",
250 "Sun subtends: ",
251 "",
252 "Last new moon:",
253 "First quarter:",
254 "Full moon: ",
255 "Last quarter: ",
256 "Next new moon:"
257 };
258
259 #define Nlabels ((sizeof labels) / sizeof(char *))
260
261 static Widget wlabel[Nlabels]; /* Widgets for displaying values */
262 static Widget loonie[2]; /* Lunation display widgets */
263 static Widget toplevel; /* Top level application shell */
264 static Widget canvas; /* Canvas for drawing moon image */
265 static Pixmap p, p1; /* Pixmaps swapped into icon */
266 static XGCValues gcv; /* Dummy values for graphics context */
267 static GC gc; /* Graphics context */
268 static XtAppContext app; /* Application context */
269 static XtIntervalId tid = (int)NULL; /* Active timeout interval ID */
270
271 static char olabel[Nlabels][60]; /* Old label values */
272 static char luabel[2][60]; /* Old lunation values */
273
274 /* Forward functions */
275 double phase2();
276 static double jtime(), phase();
277 static void phasehunt(), fmt_phase_time();
278 static void ringgg(), jyear(), jhms();
279
280 /* EXPOSE -- Graphics area repaint procedure. */
281 static void
expose(w,cd,cld)282 expose(w, cd, cld)
283 Widget w;
284 XtPointer cd, cld;
285 {
286 moon_icon_last[0] = ~moon_icon_work[0]; /* Force screen refresh */
287 ringgg(toplevel, NULL);
288 }
289
290 /* RESIZE -- Graphics area resize procedure. */
291 static void
resize(w,cd,cld)292 resize(w, cd, cld)
293 Widget w;
294 XtPointer cd, cld;
295 {
296 if (XtIsRealized(w)) {
297 XClearArea(XtDisplay(w), XtWindow(w), 0, 0, 0, 0, TRUE);
298 }
299 }
300
301 /* EPROC -- Intercept expose and UnmapNotify events to discover when
302 the window is iconised and restored to full size. There
303 must be a better way to do this, but I'll be darned if I
304 can figger it out! */
305 static void
eproc(w,client_data,event,ctd)306 eproc(w, client_data, event, ctd)
307 Widget w;
308 XtPointer client_data;
309 XEvent *event;
310 Boolean *ctd;
311 {
312 if (event->type == MapNotify)
313 iconic = FALSE;
314 else if (event->type == UnmapNotify)
315 iconic = TRUE;
316 else
317 iconic = FALSE;
318 ringgg(toplevel, NULL);
319 }
320
321 /* TESTON -- Activate test mode. */
322 static void
teston(w,client_data,call_data)323 teston(w, client_data, call_data)
324 Widget w;
325 XtPointer client_data, call_data;
326 {
327 testmode = (testmode)?FALSE:TRUE;
328 if (testmode)
329 faketime = 0.0;
330 else {
331 nptime = 0.0; /* Force recomputation of phase table */
332 expose(NULL, NULL, NULL); /* Force complete redisplay */
333 }
334 }
335
336 /* TESTOFF -- Terminate test mode. */
337 static void
testoff(w,client_data,call_data)338 testoff(w, client_data, call_data)
339 Widget w;
340 XtPointer client_data, call_data;
341 {
342 testmode = FALSE;
343 nptime = 0.0; /* Force recomputation of phase table */
344 expose(NULL, NULL, NULL); /* Force complete redisplay */
345 }
346
347 /* MAIN -- Main program. */
348 int
main(argc,argv)349 main(argc, argv)
350 int argc;
351 char *argv[];
352 {
353 int c;
354 extern int opterr;
355 Widget master_rc, canvas_rc, tmp_w, txt1_rc, txt2_rc, txt3_rc, txt4_rc,
356 testbutton, col1_rc, col2_rc, txt1_f, txt2_f, txt3_f, txt4_f;
357
358 toplevel = XtVaAppInitialize(&app, "XMoontool",
359 (XrmOptionDescList) NULL,
360 0, &argc, argv, (String *) NULL,
361 (ArgList) NULL);
362
363 opterr = 0;
364 while ((c = getopt(argc, argv, "cmtU")) != EOF) {
365 switch (c) {
366 case 't': /* Jump into test mode */
367 testmode = TRUE;
368 break;
369 case 'c': /* Force to color mode */
370 color_mode = TRUE;
371 break;
372 case 'm': /* Force mono mode */
373 color_mode = FALSE;
374 break;
375 case 'U':
376 metric = MILE;
377 break;
378 }
379 }
380 if (opterr) {
381 return 2;
382 }
383
384 #ifdef DEFAULT_COLOR
385 if (color_mode == UNKNOWN) color_mode = TRUE;
386 #endif
387
388 if (color_mode == UNKNOWN) {
389 XVisualInfo tmpl, *ret;
390 int n;
391
392 tmpl.depth = 8;
393 tmpl.screen = DefaultScreen(XtDisplay(toplevel));
394 ret = XGetVisualInfo(XtDisplay(toplevel), VisualDepthMask, &tmpl, &n);
395 if (n == 0) {
396 tmpl.depth = 1;
397 ret = XGetVisualInfo(XtDisplay(toplevel), VisualDepthMask, &tmpl, &n);
398 if (n == 0) {
399 fprintf(stderr, "Fatal Error! No visuals that we can use\n");
400 exit (2);
401 } else
402 color_mode = FALSE;
403 } else
404 color_mode = TRUE;
405 XFree(ret);
406 }
407
408 gcv.function = GXcopy;
409 gc = XtGetGC(toplevel, GCForeground | GCBackground | GCFunction, &gcv);
410
411 if (color_mode) {
412 register int i, j;
413 Colormap cmap;
414 XColor ctmp;
415
416 XtVaGetValues(toplevel, XmNcolormap, &cmap, NULL);
417 for (i=0; i<4; i++) {
418 ctmp.red = color_moon_col_r[i];
419 ctmp.green = color_moon_col_g[i];
420 ctmp.blue = color_moon_col_b[i];
421 ctmp.flags = DoRed | DoGreen | DoBlue;
422 XAllocColor(XtDisplay(toplevel), cmap, &ctmp);
423 col_vals[i] = ctmp.pixel;
424 }
425 for (i=0; i<8; i++) {
426 ctmp.red = c_r[i];
427 ctmp.green = c_g[i];
428 ctmp.blue = c_b[i];
429 ctmp.flags = DoRed | DoGreen | DoBlue;
430 XAllocColor(XtDisplay(toplevel), cmap, &ctmp);
431
432 col_vals[i+4] = ctmp.pixel;
433 }
434 p = XCreatePixmap(XtDisplay(toplevel),
435 RootWindowOfScreen(XtScreen(toplevel)), 64, 64,
436 DefaultDepthOfScreen(XtScreen(toplevel)));
437 p1 = XCreatePixmap(XtDisplay(toplevel),
438 RootWindowOfScreen(XtScreen(toplevel)), 64, 64,
439 DefaultDepthOfScreen(XtScreen(toplevel)));
440 for (i=0; i<64; i++)
441 for (j=0; j<64; j++) {
442 XSetForeground(XtDisplay(toplevel), gc,
443 col_vals[color_moon_pixels[(64*j)+i]]);
444 XDrawPoint(XtDisplay(toplevel), p, gc, i, j);
445 XDrawPoint(XtDisplay(toplevel), p1, gc, i, j);
446 }
447 } else {
448 p = XCreatePixmapFromBitmapData(XtDisplay(toplevel),
449 RootWindowOfScreen(XtScreen(toplevel)),
450 moon_icon_bits, moon_icon_width,
451 moon_icon_height,
452 BlackPixelOfScreen(XtScreen(toplevel)),
453 WhitePixelOfScreen(XtScreen(toplevel)),
454 DefaultDepthOfScreen(XtScreen(toplevel)));
455 p1 = XCreatePixmapFromBitmapData(XtDisplay(toplevel),
456 RootWindowOfScreen(XtScreen(toplevel)),
457 moon_icon_bits, moon_icon_width,
458 moon_icon_height,
459 BlackPixelOfScreen(XtScreen(toplevel)),
460 WhitePixelOfScreen(XtScreen(toplevel)),
461 DefaultDepthOfScreen(XtScreen(toplevel)));
462 }
463
464 XtVaSetValues(toplevel, XmNtitle, MOON_TITLE, XmNiconPixmap, p, NULL);
465
466 XtAddEventHandler(toplevel, ExposureMask | StructureNotifyMask,
467 TRUE, eproc, NULL);
468
469 master_rc = XtVaCreateManagedWidget("master_rc", xmRowColumnWidgetClass,
470 toplevel, XmNorientation, XmHORIZONTAL,
471 NULL);
472 col1_rc = XtVaCreateManagedWidget("col1_rc", xmRowColumnWidgetClass,
473 master_rc, XmNorientation, XmVERTICAL,
474 NULL);
475 col2_rc = XtVaCreateManagedWidget("col2_rc", xmRowColumnWidgetClass,
476 master_rc, XmNorientation, XmVERTICAL,
477 NULL);
478
479 txt1_f = XtVaCreateManagedWidget("txt1_f", xmFrameWidgetClass, col1_rc,
480 XmNshadowType, XmSHADOW_ETCHED_IN, NULL);
481 txt1_rc = XtVaCreateManagedWidget("txt1_rc", xmRowColumnWidgetClass,
482 txt1_f, XmNorientation, XmVERTICAL,
483 NULL);
484 txt2_f = XtVaCreateManagedWidget("txt2_f", xmFrameWidgetClass, col1_rc,
485 XmNshadowType, XmSHADOW_ETCHED_IN, NULL);
486 txt2_rc = XtVaCreateManagedWidget("txt2_rc", xmRowColumnWidgetClass,
487 txt2_f, XmNorientation, XmVERTICAL,
488 NULL);
489 txt3_f = XtVaCreateManagedWidget("txt3_f", xmFrameWidgetClass, col1_rc,
490 XmNshadowType, XmSHADOW_ETCHED_IN, NULL);
491 txt3_rc = XtVaCreateManagedWidget("txt3_rc", xmRowColumnWidgetClass,
492 txt3_f, XmNorientation, XmVERTICAL,
493 NULL);
494 txt4_f = XtVaCreateManagedWidget("txt4_f", xmFrameWidgetClass, col2_rc,
495 XmNshadowType, XmSHADOW_ETCHED_IN, NULL);
496 txt4_rc = XtVaCreateManagedWidget("txt4_rc", xmRowColumnWidgetClass,
497 txt4_f, XmNorientation, XmVERTICAL,
498 NULL);
499 canvas_rc = XtVaCreateManagedWidget("slave_rc", xmRowColumnWidgetClass,
500 col2_rc, XmNorientation, XmHORIZONTAL,
501 NULL);
502
503 memset(olabel, 0, sizeof olabel);
504 for (c = 0; c < Nlabels; c++) {
505 Widget lpanel, lentry;
506
507 if ((c==3) || (c==8) || (c==11))
508 continue;
509
510 lentry = XtVaCreateManagedWidget("cell", xmRowColumnWidgetClass,
511 ((c<3)?txt1_rc:((c<8)?txt2_rc:
512 ((c<11)?txt3_rc:txt4_rc))),
513 XmNorientation, XmHORIZONTAL, NULL);
514 lpanel = XtVaCreateManagedWidget("llabel", xmLabelWidgetClass,
515 lentry, XmNlabelString,
516 XmStringCreateSimple(labels[c]),
517 XmNalignment, XmALIGNMENT_BEGINNING,
518 NULL);
519 wlabel[c] = XtVaCreateManagedWidget("lvalue", xmLabelWidgetClass,
520 lentry, XmNalignment,
521 XmALIGNMENT_END, XmNlabelString,
522 XmStringCreateSimple(""), NULL);
523 if (c > 11) {
524 lpanel = XtVaCreateManagedWidget("llabel", xmLabelWidgetClass,
525 lentry, XmNlabelString,
526 XmStringCreateSimple((c==12 || c==16)?
527 "Lunation:":""),
528 XmNalignment,
529 XmALIGNMENT_BEGINNING, NULL);
530 loonie[c>12?1:0] = XtVaCreateManagedWidget("lvalue",
531 xmLabelWidgetClass,
532 lentry, XmNalignment,
533 XmALIGNMENT_END,
534 XmNlabelString,
535 XmStringCreateSimple(""),
536 NULL);
537 }
538 }
539
540 tmp_w = XtVaCreateManagedWidget("frame", xmFrameWidgetClass, canvas_rc,
541 XmNshadowType, XmSHADOW_ETCHED_IN,
542 XmNmarginWidth, 3, XmNmarginHeight, 3,
543 NULL);
544 canvas = XtVaCreateManagedWidget("canvas", xmDrawingAreaWidgetClass, tmp_w,
545 XmNwidth, 64, XmNheight, 64,
546 XmNresizePolicy, XmNONE, NULL);
547
548 XtAddCallback(canvas, XmNexposeCallback, expose, NULL);
549 XtAddCallback(canvas, XmNresizeCallback, resize, NULL);
550
551 testbutton = XtVaCreateManagedWidget("test", xmToggleButtonWidgetClass,
552 col2_rc, XmNlabelString, XmStringCreateSimple("Test"),
553 NULL);
554
555 XtAddCallback(testbutton, XmNvalueChangedCallback, teston, NULL);
556
557 XtVaSetValues(testbutton, XmNset, testmode, NULL);
558
559 ringgg(toplevel, NULL);
560 XtRealizeWidget(toplevel);
561 XtAppMainLoop(app);
562
563 return 0;
564 }
565
566 /* ROP -- Perform raster op on icon. */
567
rop(sx,sy,lx)568 static void rop(sx, sy, lx)
569 int sx, sy, lx;
570 {
571 int rowoff = sy * (moon_icon_width / 8), i;
572
573 for (i = sx; i < (sx + lx); i++) {
574 moon_icon_work[rowoff + (i / 8)] =
575 (moon_icon_work[rowoff + (i / 8)] & ~(1 << (i & 7))) |
576 (moon_icon_bits[rowoff + (i / 8)] &
577 (1 << (i & 7)));
578 }
579 }
580
rop2(sx,sy,lx,p)581 static void rop2(sx, sy, lx, p)
582 int sx, sy, lx;
583 Pixmap p;
584 {
585 int i;
586
587 for (i = sx; i < (sx + lx); i++) {
588 XSetForeground(XtDisplay(toplevel), gc,
589 col_vals[color_moon_pixels[(64*sy)+i]+4]);
590 XDrawPoint(XtDisplay(toplevel), p, gc, i, sy);
591 }
592 }
593
594 /* DRAWMOON -- Construct icon for moon, given phase of moon. */
595
drawmoon(ph,aom_d,aom_h,aom_m)596 static void drawmoon(ph, aom_d, aom_h, aom_m)
597 double ph;
598 int aom_d, aom_h, aom_m;
599 {
600 register int i, lx, rx, j, iradius;
601 register double cp, xscale, radius;
602 static char tbuf[20];
603 Pixmap np;
604
605 #define RADIUS 27.0
606 #define IRADIUS 27
607 #define OFFSET 28
608 #define CENTER 32
609
610 if (!XtIsRealized(toplevel)) {
611 return;
612 }
613
614 if (aom_d == 0) {
615 sprintf(tbuf, "%dh %dm", aom_h, aom_m);
616 } else {
617 sprintf(tbuf, "%dd %dh", aom_d, aom_h);
618 }
619
620 if (iconic || (color_mode == FALSE)) {
621 memset(moon_icon_work, 0xFF, sizeof moon_icon_work);
622
623 xscale = cos(2 * PI * ph);
624 for (i = 0; i < IRADIUS; i++) {
625 cp = RADIUS * cos(asin((double) i / RADIUS));
626 if (ph < 0.5) {
627 rx = CENTER + cp;
628 lx = CENTER + xscale * cp;
629 } else {
630 lx = CENTER - cp;
631 rx = CENTER - xscale * cp;
632 }
633
634 /* We now know the left and right endpoints of the scan line
635 for this y coordinate. We raster-op the corresponding
636 scanlines from the source pixrect to the destination
637 pixrect, offsetting to properly place it in the pixrect and
638 reflecting vertically. */
639
640 rop(lx, OFFSET + i, (rx - lx) + 1);
641 if (i != 0) {
642 rop(lx, OFFSET - i, (rx - lx) + 1);
643 }
644 }
645
646 /* The following silly little dance where we double buffer the
647 icon in two pixmaps and swap them back and forth has nothing to
648 do with efficiency. Rather, it's a work-around for the fact
649 that X doesn't update the icon unless you pass it it different
650 address for the icon pixmap than the one it had before. Note
651 also that we don't update the icon unless it's actually
652 changed. This not only saves time, it also eliminates
653 unnecessary blinking of the icon on the display due to nugatory
654 refreshes. */
655
656 if ((memcmp(moon_icon_work, moon_icon_last, sizeof moon_icon_last) != 0) || iconic) {
657 memcpy(moon_icon_last, moon_icon_work, sizeof moon_icon_last);
658 np = XCreatePixmapFromBitmapData(XtDisplay(toplevel),
659 RootWindowOfScreen(XtScreen(toplevel)),
660 moon_icon_work, moon_icon_width,
661 moon_icon_height,
662 BlackPixelOfScreen(XtScreen(toplevel)),
663 WhitePixelOfScreen(XtScreen(toplevel)),
664 DefaultDepthOfScreen(XtScreen(toplevel)));
665 XCopyArea(XtDisplay(toplevel), np, p1, gc, 0, 0, moon_icon_width,
666 moon_icon_height, 0, 0);
667 if (XtIsRealized(canvas)) {
668 XCopyArea(XtDisplay(canvas), np, XtWindow(canvas), gc, 0, 0,
669 moon_icon_width, moon_icon_height, 0, 0);
670 }
671 XFreePixmap(XtDisplay(toplevel), np);
672 XtVaSetValues(toplevel, XmNiconPixmap, p1, NULL);
673
674 np = p;
675 p = p1;
676 p1 = np;
677 }
678 } else {
679 int lval;
680
681 memset(moon_icon_work, 0xFF, sizeof moon_icon_work);
682
683 /* This makes the moon appear completely dark when new. */
684 radius = (ph < 0.01 || ph > 0.99) ? RADIUS + 1.0 : RADIUS;
685 iradius = radius;
686
687 xscale = cos(2 * PI * ph);
688 for (i = 0; i < iradius; i++) {
689 cp = radius * cos(asin((double) i / radius));
690 if (ph < 0.5) {
691 rx = CENTER + cp;
692 lx = CENTER + xscale * cp;
693 } else {
694 lx = CENTER - cp;
695 rx = CENTER - xscale * cp;
696 }
697
698 /* We now know the left and right endpoints of the scan line
699 for this y coordinate. We raster-op the corresponding
700 scanlines from the source pixrect to the destination
701 pixrect, offsetting to properly place it in the pixrect and
702 reflecting vertically. */
703
704 rop(lx, OFFSET + i, (rx - lx) + 1);
705 if (i != 0) {
706 rop(lx, OFFSET - i, (rx - lx) + 1);
707 }
708 }
709 if (memcmp(moon_icon_work, moon_icon_last, sizeof moon_icon_last) != 0) {
710 memcpy(moon_icon_last, moon_icon_work, sizeof moon_icon_last);
711 np = XCreatePixmap(XtDisplay(toplevel),
712 RootWindowOfScreen(XtScreen(toplevel)), 64, 64,
713 DefaultDepthOfScreen(XtScreen(toplevel)));
714 for (i=0; i<64; i++)
715 for (j=0; j<64; j++) {
716 lval = color_moon_pixels[(64*j)+i];
717 XSetForeground(XtDisplay(toplevel), gc,
718 col_vals[(lval>0)?lval+8:4]);
719 XDrawPoint(XtDisplay(toplevel), np, gc, i, j);
720 }
721
722 for (i=0; i < iradius; i++) {
723 cp = radius * cos(asin((double)i / radius));
724 if (ph < 0.5) {
725 rx = CENTER + cp;
726 lx = CENTER + xscale * cp;
727 } else {
728 lx = CENTER - cp;
729 rx = CENTER - xscale * cp;
730 }
731 rop2(lx, OFFSET+i, (rx-lx)+1, np);
732 if (i != 0)
733 rop2(lx, OFFSET-i, (rx-lx)+1, np);
734 }
735 XSetForeground(XtDisplay(toplevel), gc, col_vals[8]);
736 XDrawPoint(XtDisplay(toplevel), np, gc, 41, 29);
737 XCopyArea(XtDisplay(toplevel), np, p1, gc, 0, 0, moon_icon_width,
738 moon_icon_height, 0, 0);
739 if (XtIsRealized(canvas))
740 XCopyArea(XtDisplay(canvas), np, XtWindow(canvas), gc, 0, 0,
741 moon_icon_width, moon_icon_height, 0, 0);
742 XFreePixmap(XtDisplay(toplevel), np);
743 XtVaSetValues(toplevel, XmNiconPixmap, p1, NULL);
744
745 np = p;
746 p = p1;
747 p1 = np;
748 }
749 }
750
751 /* Update the date of the moon in the icon label if it's changed. */
752
753 if (strcmp(tbuf, icon_label_last) != 0) {
754 strcpy(icon_label_last, tbuf);
755 XtVaSetValues(toplevel, XmNiconName, tbuf, NULL);
756 }
757 }
758
759 /* RINGGG -- Update status on interval timer ticks and redraw
760 window if needed. */
761
762 #define prt(y) if (strcmp(olabel[y - 1], tbuf) != 0) { \
763 XtVaSetValues(wlabel[y - 1], XmNlabelString, XmStringCreateSimple(tbuf), NULL); \
764 strcpy(olabel[y - 1], tbuf); }
765
766 #define prl(y) if (strcmp(luabel[y - 1], tbuf) != 0) { \
767 XtVaSetValues(loonie[y - 1], XmNlabelString, XmStringCreateSimple(tbuf), NULL); \
768 strcpy(luabel[y - 1], tbuf); }
769
770 #define EPL(x) (x), (x) == 1 ? "" : "s"
771 #define APOS(x) (x + 13)
772
ringgg(client_data,id)773 static void ringgg(client_data, id)
774 XtPointer client_data;
775 XtIntervalId *id;
776 {
777 int lunation;
778 int i, yy, mm, dd, hh, mmm, ss;
779 int aom_d, aom_h, aom_m;
780 time_t t;
781 double jd, p, aom, cphase, cdist, cangdia, csund, csuang;
782 double phasar[5];
783 char tbuf[80];
784 struct tm *gm;
785
786 time(&t);
787 jd = jtime((gm = gmtime(&t)));
788 if (testmode) {
789 if (faketime == 0.0) {
790 faketime = jd + 1;
791 } else {
792 faketime += 1.0 / 24;
793 }
794 jd = faketime;
795 }
796
797 p = phase2(jd, &cphase, &aom, &cdist, &cangdia, &csund, &csuang);
798
799 /* aom as returned by phase() is at least confusing. The aom_m update
800 rate can be observed to be fairly different from 1 min so that the
801 displayed information is not actually reliable.
802 Instead we let the current aom be the time elapsed since the recent
803 new moon.
804 */
805 if (jd > nptime || jd < pptime) {
806 phasehunt(jd, phasar);
807
808 /* Check the result and, if desired, try to correct it. */
809 if (phasar[0] > jd) {
810 phasehunt(jd - 1.5, phasar);
811 }
812 if (phasar[4] < jd) {
813 phasehunt(jd + 1.5, phasar);
814 }
815 pptime = phasar[0];
816 }
817
818 aom = jd - pptime;
819 aom_d = (int) aom;
820 aom_h = (int) (24 * (aom - floor(aom)));
821 aom_m = (int) (1440 * (aom - floor(aom))) % 60;
822
823 drawmoon(p, aom_d, aom_h, aom_m);
824
825 /* If we were called preemptively, cancel the pending unexpired
826 timeout. */
827
828 if (id == NULL && tid != (int)NULL) {
829 XtRemoveTimeOut(tid);
830 }
831
832 tid = XtAppAddTimeOut(app, testmode ? 150 :
833 (id == NULL ? 10 : ((iconic?CLOSED_SECS:OPEN_SECS) * 1000)),
834 ringgg, (Widget) client_data);
835
836 if (iconic)
837 return;
838
839 /* Update textual information for open window */
840
841 sprintf(tbuf, "%.5f", jd); /* Julian date */
842 prt(1);
843
844 if (testmode) { /* Universal time */
845 jyear(jd, &yy, &mm, &dd);
846 jhms(jd, &hh, &mmm, &ss);
847 sprintf(tbuf, "%02d:%02d:%02d %02d %s %d UTC",
848 hh, mmm, ss, dd, moname [mm - 1], yy);
849 } else {
850 strftime(tbuf, sizeof(tbuf)-1, "%H:%M:%S %d %B %Y %Z", gm);
851 }
852 prt(2);
853
854 if (testmode == FALSE) { /* Ignore local time in test mode */
855 gm = localtime(&t);
856
857 /* Local time */
858
859 strftime(tbuf, sizeof(tbuf)-1, "%H:%M:%S %d %B %Y %Z", gm);
860 prt(3);
861 }
862
863 /* Information about the Moon */
864
865 /* Age of moon */
866
867 sprintf(tbuf, "%d day%s, %d hour%s, %d minute%s.",
868 EPL(aom_d), EPL(aom_h), EPL(aom_m));
869 prt(5);
870
871 /* Moon phase */
872
873 sprintf(tbuf, "%.0f%% (0%% = New, 100%% = Full)", cphase * 100);
874 while (strlen(tbuf) < 40) {
875 strcat(tbuf, " ");
876 }
877 prt(6);
878
879 /* Moon distance */
880
881 sprintf(tbuf, "%.0f %s, %.1f Earth radii.",
882 cdist * metric, units(metric), cdist / earthrad);
883 prt(7);
884
885 /* Moon subtends */
886
887 sprintf(tbuf, "%.4f degrees.", cangdia);
888 prt(8);
889
890 /* Information about the Sun */
891
892 /* Sun's distance */
893
894 sprintf(tbuf, "%.0f %s, %.3f astronomical units.",
895 csund * metric, units(metric), csund / sunsmax);
896 prt(10);
897
898 /* Sun subtends */
899
900 sprintf(tbuf, "%.4f degrees.", csuang);
901 prt(11);
902
903 /* Calculate times of phases of this lunation. This is
904 sufficiently time-consuming that we only do it once a month. */
905
906 if (jd > nptime) {
907 phasehunt(jd, phasar);
908
909 /* Check the result and, if desired, try to correct it. */
910 if (phasar[0] > jd) {
911 phasehunt(jd - 1.5, phasar);
912 }
913 if (phasar[4] < jd) {
914 phasehunt(jd + 1.5, phasar);
915 }
916
917 lunation = floor(((phasar[0] + 7) - lunatbase) / synmonth) + 1;
918
919 for (i = 0; i < 5; i++) {
920 fmt_phase_time(phasar[i], tbuf, sizeof(tbuf));
921 while (strlen(tbuf) < 23) {
922 strcat(tbuf, " ");
923 }
924 prt(APOS(i));
925 }
926 nptime = phasar[4];
927
928 /* Edit lunation numbers into cells reserved for them. */
929
930 sprintf(tbuf, "%d", lunation);
931 prl(1);
932 sprintf(tbuf, "%d", lunation + 1);
933 prl(2);
934 }
935 return;
936 }
937 #undef APOS
938
939
940 /* FMT_PHASE_TIME -- Format the provided julian date into the
941 provided buffer in UTC format for screen
942 display */
943
fmt_phase_time(utime,buf,len)944 static void fmt_phase_time(utime, buf, len)
945 double utime;
946 char *buf;
947 size_t len;
948 {
949 int yy, mm, dd, hh, mmm, ss;
950 time_t pt;
951 struct tm ph, *lt;
952
953 jyear(utime, &yy, &mm, &dd);
954 jhms(utime, &hh, &mmm, &ss);
955
956 ph.tm_sec = ss;
957 ph.tm_min = mmm;
958 ph.tm_hour = hh;
959 ph.tm_year = yy - 1900;
960 ph.tm_mon = mm - 1;
961 ph.tm_mday = dd;
962
963 pt = timegm(&ph);
964 lt = localtime(&pt);
965 strftime(buf, len-1, "%d %b %Y %H:%M %Z", lt);
966 }
967
968
969 /* JDATE -- Convert internal GMT date and time to Julian day and
970 fraction. */
971
jdate(t)972 static long jdate(t)
973 struct tm *t;
974 {
975 long c, m, y;
976
977 y = t->tm_year + 1900;
978 m = t->tm_mon + 1;
979 if (m > 2) {
980 m = m - 3;
981 } else {
982 m = m + 9;
983 y--;
984 }
985 c = y / 100L; /* Compute century */
986 y -= 100L * c;
987 return (t->tm_mday + (c * 146097L) / 4 + (y * 1461L) / 4 +
988 (m * 153L + 2) / 5 + 1721119L);
989 }
990
991
992 /* JTIME -- Convert internal GMT date and time to astronomical
993 Julian time (i.e. Julian date plus day fraction,
994 expressed as a double). */
995
jtime(t)996 static double jtime(t)
997 struct tm *t;
998 {
999 return (jdate(t) - 0.5) +
1000 (t->tm_sec + 60 * (t->tm_min + 60 * t->tm_hour)) / 86400.0;
1001 }
1002
1003
1004 /* JYEAR -- Convert Julian date to year, month, day, which are
1005 returned via integer pointers to integers. */
1006
jyear(td,yy,mm,dd)1007 static void jyear(td, yy, mm, dd)
1008 double td;
1009 int *yy, *mm, *dd;
1010 {
1011 double j, d, y, m;
1012
1013 td += 0.5; /* Astronomical to civil */
1014 j = floor(td);
1015 j = j - 1721119.0;
1016 y = floor(((4 * j) - 1) / 146097.0);
1017 j = (j * 4.0) - (1.0 + (146097.0 * y));
1018 d = floor(j / 4.0);
1019 j = floor(((4.0 * d) + 3.0) / 1461.0);
1020 d = ((4.0 * d) + 3.0) - (1461.0 * j);
1021 d = floor((d + 4.0) / 4.0);
1022 m = floor(((5.0 * d) - 3) / 153.0);
1023 d = (5.0 * d) - (3.0 + (153.0 * m));
1024 d = floor((d + 5.0) / 5.0);
1025 y = (100.0 * y) + j;
1026 if (m < 10.0) {
1027 m = m + 3;
1028 } else {
1029 m = m - 9;
1030 y = y + 1;
1031 }
1032 *yy = y;
1033 *mm = m;
1034 *dd = d;
1035 }
1036
1037
1038 /* JHMS -- Convert Julian time to hour, minutes, and seconds. */
1039
jhms(j,h,m,s)1040 static void jhms(j, h, m, s)
1041 double j;
1042 int *h, *m, *s;
1043 {
1044 long ij;
1045
1046 j += 0.5; /* Astronomical to civil */
1047 ij = (j - floor(j)) * 86400.0;
1048 *h = ij / 3600L;
1049 *m = (ij / 60L) % 60L;
1050 *s = ij % 60L;
1051 }
1052
1053
1054 /* MEANPHASE -- Calculates time of the mean new Moon for a given
1055 base date. This argument K to this function is the
1056 precomputed synodic month index, given by:
1057
1058 K = (year - 1900) * 12.3685
1059
1060 where year is expressed as a year and fractional year. */
1061
meanphase(sdate,k)1062 static double meanphase(sdate, k)
1063 double sdate, k;
1064 {
1065 double t, t2, t3, nt1;
1066
1067 /* Time in Julian centuries from 1900 January 0.5 */
1068 t = (sdate - 2415020.0) / 36525;
1069 t2 = t * t; /* Square for frequent use */
1070 t3 = t2 * t; /* Cube for frequent use */
1071
1072 nt1 = 2415020.75933 + synmonth * k
1073 + 0.0001178 * t2
1074 - 0.000000155 * t3
1075 + 0.00033 * dsin(166.56 + 132.87 * t - 0.009173 * t2);
1076
1077 return nt1;
1078 }
1079
1080
1081 /* TRUEPHASE -- Given a K value used to determine the mean phase of
1082 the new moon, and a phase selector (0.0, 0.25, 0.5,
1083 0.75), obtain the true, corrected phase time. */
1084
truephase(k,phase)1085 static double truephase(k, phase)
1086 double k, phase;
1087 {
1088 double t, t2, t3, pt, m, mprime, f;
1089 int apcor = FALSE;
1090
1091 k += phase; /* Add phase to new moon time */
1092 t = k / 1236.85; /* Time in Julian centuries from
1093 1900 January 0.5 */
1094 t2 = t * t; /* Square for frequent use */
1095 t3 = t2 * t; /* Cube for frequent use */
1096 pt = 2415020.75933 /* Mean time of phase */
1097 + synmonth * k
1098 + 0.0001178 * t2
1099 - 0.000000155 * t3
1100 + 0.00033 * dsin(166.56 + 132.87 * t - 0.009173 * t2);
1101
1102 m = 359.2242 /* Sun's mean anomaly */
1103 + 29.10535608 * k
1104 - 0.0000333 * t2
1105 - 0.00000347 * t3;
1106 mprime = 306.0253 /* Moon's mean anomaly */
1107 + 385.81691806 * k
1108 + 0.0107306 * t2
1109 + 0.00001236 * t3;
1110 f = 21.2964 /* Moon's argument of latitude */
1111 + 390.67050646 * k
1112 - 0.0016528 * t2
1113 - 0.00000239 * t3;
1114 if ((phase < 0.01) || (abs(phase - 0.5) < 0.01)) {
1115
1116 /* Corrections for New and Full Moon */
1117
1118 pt += (0.1734 - 0.000393 * t) * dsin(m)
1119 + 0.0021 * dsin(2 * m)
1120 - 0.4068 * dsin(mprime)
1121 + 0.0161 * dsin(2 * mprime)
1122 - 0.0004 * dsin(3 * mprime)
1123 + 0.0104 * dsin(2 * f)
1124 - 0.0051 * dsin(m + mprime)
1125 - 0.0074 * dsin(m - mprime)
1126 + 0.0004 * dsin(2 * f + m)
1127 - 0.0004 * dsin(2 * f - m)
1128 - 0.0006 * dsin(2 * f + mprime)
1129 + 0.0010 * dsin(2 * f - mprime)
1130 + 0.0005 * dsin(m + 2 * mprime);
1131 apcor = TRUE;
1132 } else if ((abs(phase - 0.25) < 0.01 || (abs(phase - 0.75) < 0.01))) {
1133 pt += (0.1721 - 0.0004 * t) * dsin(m)
1134 + 0.0021 * dsin(2 * m)
1135 - 0.6280 * dsin(mprime)
1136 + 0.0089 * dsin(2 * mprime)
1137 - 0.0004 * dsin(3 * mprime)
1138 + 0.0079 * dsin(2 * f)
1139 - 0.0119 * dsin(m + mprime)
1140 - 0.0047 * dsin(m - mprime)
1141 + 0.0003 * dsin(2 * f + m)
1142 - 0.0004 * dsin(2 * f - m)
1143 - 0.0006 * dsin(2 * f + mprime)
1144 + 0.0021 * dsin(2 * f - mprime)
1145 + 0.0003 * dsin(m + 2 * mprime)
1146 + 0.0004 * dsin(m - 2 * mprime)
1147 - 0.0003 * dsin(2 * m + mprime);
1148 if (phase < 0.5)
1149 /* First quarter correction */
1150 pt += 0.0028 - 0.0004 * dcos(m) + 0.0003 * dcos(mprime);
1151 else
1152 /* Last quarter correction */
1153 pt += -0.0028 + 0.0004 * dcos(m) - 0.0003 * dcos(mprime);
1154 apcor = TRUE;
1155 }
1156 if (!apcor) {
1157 fprintf(stderr,
1158 "TRUEPHASE called with invalid phase selector.\n");
1159 abort();
1160 }
1161 return pt;
1162 }
1163
1164
1165 /* PHASEHUNT -- Find time of phases of the moon which surround the
1166 current date. Five phases are found, starting and
1167 ending with the new moons which bound the current
1168 lunation. */
1169
phasehunt(sdate,phases)1170 static void phasehunt(sdate, phases)
1171 double sdate;
1172 double phases[5];
1173 {
1174 double adate, k1, k2, nt1, nt2;
1175 int yy, mm, dd;
1176
1177 adate = sdate - 45;
1178
1179 jyear(adate, &yy, &mm, &dd);
1180 k1 = floor((yy + ((mm - 1) * (1.0 / 12.0)) - 1900) * 12.3685);
1181
1182 adate = nt1 = meanphase(adate, k1);
1183 while (TRUE) {
1184 adate += synmonth;
1185 k2 = k1 + 1;
1186 nt2 = meanphase(adate, k2);
1187 if (nt1 <= sdate && nt2 > sdate)
1188 break;
1189 nt1 = nt2;
1190 k1 = k2;
1191 }
1192 phases[0] = truephase(k1, 0.0);
1193 phases[1] = truephase(k1, 0.25);
1194 phases[2] = truephase(k1, 0.5);
1195 phases[3] = truephase(k1, 0.75);
1196 phases[4] = truephase(k2, 0.0);
1197 }
1198
1199
1200 /* KEPLER -- Solve the equation of Kepler. */
1201
kepler(m,ecc)1202 static double kepler(m, ecc)
1203 double m, ecc;
1204 {
1205 double e, delta;
1206 #define EPSILON 1E-6
1207
1208 e = m = torad(m);
1209 do {
1210 delta = e - ecc * sin(e) - m;
1211 e -= delta / (1 - ecc * cos(e));
1212 } while (abs(delta) > EPSILON);
1213 return e;
1214 }
1215
1216 /* PHASE -- Calculate phase of moon as a fraction:
1217
1218 The argument is the time for which the phase is requested,
1219 expressed as a Julian date and fraction. Returns the terminator
1220 phase angle as a percentage of a full circle (i.e., 0 to 1), and
1221 stores into pointer arguments the illuminated fraction of the
1222 Moon's disc, the Moon's age in days and fraction, the distance of
1223 the Moon from the centre of the Earth, and the angular diameter
1224 subtended by the Moon as seen by an observer at the centre of the
1225 Earth.
1226 */
1227
phase(pdate,pphase,mage,dist,angdia,sudist,suangdia)1228 static double phase(pdate, pphase, mage, dist, angdia, sudist, suangdia)
1229 double pdate;
1230 double *pphase; /* Illuminated fraction */
1231 double *mage; /* Age of moon in days */
1232 double *dist; /* Distance in kilometres */
1233 double *angdia; /* Angular diameter in degrees */
1234 double *sudist; /* Distance to Sun */
1235 double *suangdia; /* Sun's angular diameter */
1236 {
1237
1238 double Day, N, M, Ec, Lambdasun, ml, MM, MN, Ev, Ae, A3, MmP,
1239 mEc, A4, lP, V, lPP, NP, y, x, Lambdamoon, BetaM,
1240 MoonAge, MoonPhase,
1241 MoonDist, MoonDFrac, MoonAng, MoonPar,
1242 F, SunDist, SunAng;
1243
1244 /* Calculation of the Sun's position */
1245
1246 Day = pdate - epoch; /* Date within epoch */
1247 N = fixangle((360 / 365.2422) * Day); /* Mean anomaly of the Sun */
1248 M = fixangle(N + elonge - elongp); /* Convert from perigee
1249 co-ordinates to epoch 1980.0 */
1250 Ec = kepler(M, eccent); /* Solve equation of Kepler */
1251 Ec = sqrt((1 + eccent) / (1 - eccent)) * tan(Ec / 2);
1252 Ec = 2 * todeg(atan(Ec)); /* True anomaly */
1253 Lambdasun = fixangle(Ec + elongp); /* Sun's geocentric ecliptic
1254 longitude */
1255 /* Orbital distance factor */
1256 F = ((1 + eccent * cos(torad(Ec))) / (1 - eccent * eccent));
1257 SunDist = sunsmax / F; /* Distance to Sun in km */
1258 SunAng = F * sunangsiz; /* Sun's angular size in degrees */
1259
1260
1261 /* Calculation of the Moon's position */
1262
1263 /* Moon's mean longitude */
1264 ml = fixangle(13.1763966 * Day + mmlong);
1265
1266 /* Moon's mean anomaly */
1267 MM = fixangle(ml - 0.1114041 * Day - mmlongp);
1268
1269 /* Moon's ascending node mean longitude */
1270 MN = fixangle(mlnode - 0.0529539 * Day);
1271
1272 /* Evection */
1273 Ev = 1.2739 * sin(torad(2 * (ml - Lambdasun) - MM));
1274
1275 /* Annual equation */
1276 Ae = 0.1858 * sin(torad(M));
1277
1278 /* Correction term */
1279 A3 = 0.37 * sin(torad(M));
1280
1281 /* Corrected anomaly */
1282 MmP = MM + Ev - Ae - A3;
1283
1284 /* Correction for the equation of the centre */
1285 mEc = 6.2886 * sin(torad(MmP));
1286
1287 /* Another correction term */
1288 A4 = 0.214 * sin(torad(2 * MmP));
1289
1290 /* Corrected longitude */
1291 lP = ml + Ev + mEc - Ae + A4;
1292
1293 /* Variation */
1294 V = 0.6583 * sin(torad(2 * (lP - Lambdasun)));
1295
1296 /* True longitude */
1297 lPP = lP + V;
1298
1299 /* Corrected longitude of the node */
1300 NP = MN - 0.16 * sin(torad(M));
1301
1302 /* Y inclination coordinate */
1303 y = sin(torad(lPP - NP)) * cos(torad(minc));
1304
1305 /* X inclination coordinate */
1306 x = cos(torad(lPP - NP));
1307
1308 /* Ecliptic longitude */
1309 Lambdamoon = todeg(atan2(y, x));
1310 Lambdamoon += NP;
1311
1312 /* Ecliptic latitude */
1313 BetaM = todeg(asin(sin(torad(lPP - NP)) * sin(torad(minc))));
1314
1315 /* Calculation of the phase of the Moon */
1316
1317 /* Age of the Moon in degrees */
1318 MoonAge = lPP - Lambdasun;
1319
1320 /* Phase of the Moon */
1321 MoonPhase = (1 - cos(torad(MoonAge))) / 2;
1322
1323 /* Calculate distance of moon from the centre of the Earth */
1324
1325 MoonDist = (msmax * (1 - mecc * mecc)) /
1326 (1 + mecc * cos(torad(MmP + mEc)));
1327
1328 /* Calculate Moon's angular diameter */
1329
1330 MoonDFrac = MoonDist / msmax;
1331 MoonAng = mangsiz / MoonDFrac;
1332
1333 /* Calculate Moon's parallax */
1334
1335 MoonPar = mparallax / MoonDFrac;
1336
1337 *pphase = MoonPhase;
1338 *mage = synmonth * (fixangle(MoonAge) / 360.0);
1339 *dist = MoonDist;
1340 *angdia = MoonAng;
1341 *sudist = SunDist;
1342 *suangdia = SunAng;
1343 return fixangle(MoonAge) / 360.0;
1344 }
1345