1 /*
2 * Copyright (c) 1989, 1993
3 * The Regents of the University of California. All rights reserved.
4 *
5 * This code is derived from software posted to USENET.
6 *
7 * %sccs.include.redist.c%
8 */
9
10 #ifndef lint
11 static char copyright[] =
12 "@(#) Copyright (c) 1989, 1993\n\
13 The Regents of the University of California. All rights reserved.\n";
14 #endif /* not lint */
15
16 #ifndef lint
17 static char sccsid[] = "@(#)pom.c 8.1 (Berkeley) 05/31/93";
18 #endif /* not lint */
19
20 /*
21 * Phase of the Moon. Calculates the current phase of the moon.
22 * Based on routines from `Practical Astronomy with Your Calculator',
23 * by Duffett-Smith. Comments give the section from the book that
24 * particular piece of code was adapted from.
25 *
26 * -- Keith E. Brandt VIII 1984
27 *
28 */
29
30 #include <sys/time.h>
31 #include <stdio.h>
32 #include <tzfile.h>
33 #include <math.h>
34
35 #define PI 3.141592654
36 #define EPOCH 85
37 #define EPSILONg 279.611371 /* solar ecliptic long at EPOCH */
38 #define RHOg 282.680403 /* solar ecliptic long of perigee at EPOCH */
39 #define ECCEN 0.01671542 /* solar orbit eccentricity */
40 #define lzero 18.251907 /* lunar mean long at EPOCH */
41 #define Pzero 192.917585 /* lunar mean long of perigee at EPOCH */
42 #define Nzero 55.204723 /* lunar mean long of node at EPOCH */
43
44 double dtor(), potm(), adj360();
45
main()46 main()
47 {
48 extern int errno;
49 struct timeval tp;
50 struct timezone tzp;
51 struct tm *GMT, *gmtime();
52 double days, today, tomorrow;
53 int cnt;
54 char *strerror();
55
56 if (gettimeofday(&tp,&tzp)) {
57 (void)fprintf(stderr, "pom: %s\n", strerror(errno));
58 exit(1);
59 }
60 GMT = gmtime(&tp.tv_sec);
61 days = (GMT->tm_yday + 1) + ((GMT->tm_hour +
62 (GMT->tm_min / 60.0) + (GMT->tm_sec / 3600.0)) / 24.0);
63 for (cnt = EPOCH; cnt < GMT->tm_year; ++cnt)
64 days += isleap(cnt) ? 366 : 365;
65 today = potm(days) + .5;
66 (void)printf("The Moon is ");
67 if ((int)today == 100)
68 (void)printf("Full\n");
69 else if (!(int)today)
70 (void)printf("New\n");
71 else {
72 tomorrow = potm(days + 1);
73 if ((int)today == 50)
74 (void)printf("%s\n", tomorrow > today ?
75 "at the First Quarter" : "at the Last Quarter");
76 else {
77 (void)printf("%s ", tomorrow > today ?
78 "Waxing" : "Waning");
79 if (today > 50)
80 (void)printf("Gibbous (%1.0f%% of Full)\n",
81 today);
82 else if (today < 50)
83 (void)printf("Crescent (%1.0f%% of Full)\n",
84 today);
85 }
86 }
87 }
88
89 /*
90 * potm --
91 * return phase of the moon
92 */
93 double
potm(days)94 potm(days)
95 double days;
96 {
97 double N, Msol, Ec, LambdaSol, l, Mm, Ev, Ac, A3, Mmprime;
98 double A4, lprime, V, ldprime, D, Nm;
99
100 N = 360 * days / 365.2422; /* sec 42 #3 */
101 adj360(&N);
102 Msol = N + EPSILONg - RHOg; /* sec 42 #4 */
103 adj360(&Msol);
104 Ec = 360 / PI * ECCEN * sin(dtor(Msol)); /* sec 42 #5 */
105 LambdaSol = N + Ec + EPSILONg; /* sec 42 #6 */
106 adj360(&LambdaSol);
107 l = 13.1763966 * days + lzero; /* sec 61 #4 */
108 adj360(&l);
109 Mm = l - (0.1114041 * days) - Pzero; /* sec 61 #5 */
110 adj360(&Mm);
111 Nm = Nzero - (0.0529539 * days); /* sec 61 #6 */
112 adj360(&Nm);
113 Ev = 1.2739 * sin(dtor(2*(l - LambdaSol) - Mm)); /* sec 61 #7 */
114 Ac = 0.1858 * sin(dtor(Msol)); /* sec 61 #8 */
115 A3 = 0.37 * sin(dtor(Msol));
116 Mmprime = Mm + Ev - Ac - A3; /* sec 61 #9 */
117 Ec = 6.2886 * sin(dtor(Mmprime)); /* sec 61 #10 */
118 A4 = 0.214 * sin(dtor(2 * Mmprime)); /* sec 61 #11 */
119 lprime = l + Ev + Ec - Ac + A4; /* sec 61 #12 */
120 V = 0.6583 * sin(dtor(2 * (lprime - LambdaSol))); /* sec 61 #13 */
121 ldprime = lprime + V; /* sec 61 #14 */
122 D = ldprime - LambdaSol; /* sec 63 #2 */
123 return(50 * (1 - cos(dtor(D)))); /* sec 63 #3 */
124 }
125
126 /*
127 * dtor --
128 * convert degrees to radians
129 */
130 double
dtor(deg)131 dtor(deg)
132 double deg;
133 {
134 return(deg * PI / 180);
135 }
136
137 /*
138 * adj360 --
139 * adjust value so 0 <= deg <= 360
140 */
141 double
adj360(deg)142 adj360(deg)
143 double *deg;
144 {
145 for (;;)
146 if (*deg < 0)
147 *deg += 360;
148 else if (*deg > 360)
149 *deg -= 360;
150 else
151 break;
152 }
153