1 /* $NetBSD: pom.c,v 1.21 2021/05/02 12:50:46 rillig Exp $ */
2
3 /*
4 * Copyright (c) 1989, 1993
5 * The Regents of the University of California. All rights reserved.
6 *
7 * This code is derived from software posted to USENET.
8 *
9 * Redistribution and use in source and binary forms, with or without
10 * modification, are permitted provided that the following conditions
11 * are met:
12 * 1. Redistributions of source code must retain the above copyright
13 * notice, this list of conditions and the following disclaimer.
14 * 2. Redistributions in binary form must reproduce the above copyright
15 * notice, this list of conditions and the following disclaimer in the
16 * documentation and/or other materials provided with the distribution.
17 * 3. Neither the name of the University nor the names of its contributors
18 * may be used to endorse or promote products derived from this software
19 * without specific prior written permission.
20 *
21 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
22 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
23 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
24 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
25 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
26 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
27 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
28 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
29 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
30 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
31 * SUCH DAMAGE.
32 */
33
34 #include <sys/cdefs.h>
35 #ifndef lint
36 __COPYRIGHT("@(#) Copyright (c) 1989, 1993\
37 The Regents of the University of California. All rights reserved.");
38 #endif /* not lint */
39
40 #ifndef lint
41 #if 0
42 static char sccsid[] = "@(#)pom.c 8.1 (Berkeley) 5/31/93";
43 #else
44 __RCSID("$NetBSD: pom.c,v 1.21 2021/05/02 12:50:46 rillig Exp $");
45 #endif
46 #endif /* not lint */
47
48 /*
49 * Phase of the Moon. Calculates the current phase of the moon.
50 * Based on routines from `Practical Astronomy with Your Calculator',
51 * by Duffett-Smith. Comments give the section from the book that
52 * particular piece of code was adapted from.
53 *
54 * -- Keith E. Brandt VIII 1984
55 *
56 * Updated to the Third Edition of Duffett-Smith's book, Paul Janzen, IX 1998
57 *
58 */
59
60 #include <ctype.h>
61 #include <err.h>
62 #include <math.h>
63 #include <stdio.h>
64 #include <string.h>
65 #include <stdlib.h>
66 #include <time.h>
67 #include <unistd.h>
68
69 #ifndef PI
70 #define PI 3.14159265358979323846
71 #endif
72
73 /*
74 * The EPOCH in the third edition of the book is 1990 Jan 0.0 TDT.
75 * In this program, we do not bother to correct for the differences
76 * between UTC (as shown by the UNIX clock) and TDT. (TDT = TAI + 32.184s;
77 * TAI-UTC = 32s in Jan 1999.)
78 */
79 #define EPOCH_MINUS_1970 (20 * 365 + 5 - 1) /* 20 years, 5 leaps, back 1 day to Jan 0 */
80 #define EPSILONg 279.403303 /* solar ecliptic long at EPOCH */
81 #define RHOg 282.768422 /* solar ecliptic long of perigee at EPOCH */
82 #define ECCEN 0.016713 /* solar orbit eccentricity */
83 #define lzero 318.351648 /* lunar mean long at EPOCH */
84 #define Pzero 36.340410 /* lunar mean long of perigee at EPOCH */
85 #define Nzero 318.510107 /* lunar mean long of node at EPOCH */
86
87 int main(int, char *[]);
88 static void adj360(double *);
89 static double dtor(double);
90 static double potm(double);
91 static time_t parsetime(char *);
92 static void badformat(void) __dead;
93
94 int
main(int argc,char * argv[])95 main(int argc, char *argv[])
96 {
97 time_t tmpt, now;
98 double days, today, tomorrow;
99 char buf[1024];
100
101 if (time(&now) == (time_t)-1)
102 err(1, "time");
103 if (argc > 1) {
104 tmpt = parsetime(argv[1]);
105 strftime(buf, sizeof(buf), "%a %Y %b %e %H:%M:%S (%Z)",
106 localtime(&tmpt));
107 printf("%s: ", buf);
108 } else {
109 tmpt = now;
110 }
111 days = (tmpt - EPOCH_MINUS_1970 * 86400) / 86400.0;
112 today = potm(days) + .5;
113 if (tmpt < now)
114 (void)printf("The Moon was ");
115 else if (tmpt == now)
116 (void)printf("The Moon is ");
117 else
118 (void)printf("The Moon will be ");
119 if ((int)today == 100)
120 (void)printf("Full\n");
121 else if (!(int)today)
122 (void)printf("New\n");
123 else {
124 tomorrow = potm(days + 1);
125 if ((int)today == 50)
126 (void)printf("%s\n", tomorrow > today ?
127 "at the First Quarter" : "at the Last Quarter");
128 /* today is 0.5 too big, but it doesn't matter here
129 * since the phase is changing fast enough
130 */
131 else {
132 today -= 0.5; /* Now it might matter */
133 (void)printf("%s ", tomorrow > today ?
134 "Waxing" : "Waning");
135 if (today > 50)
136 (void)printf("Gibbous (%1.0f%% of Full)\n",
137 today);
138 else if (today < 50)
139 (void)printf("Crescent (%1.0f%% of Full)\n",
140 today);
141 }
142 }
143
144 return EXIT_SUCCESS;
145 }
146
147 /*
148 * potm --
149 * return phase of the moon
150 */
151 static double
potm(double days)152 potm(double days)
153 {
154 double N, Msol, Ec, LambdaSol, l, Mm, Ev, Ac, A3, Mmprime;
155 double A4, lprime, V, ldprime, D, Nm;
156
157 N = 360 * days / 365.242191; /* sec 46 #3 */
158 adj360(&N);
159 Msol = N + EPSILONg - RHOg; /* sec 46 #4 */
160 adj360(&Msol);
161 Ec = 360 / PI * ECCEN * sin(dtor(Msol)); /* sec 46 #5 */
162 LambdaSol = N + Ec + EPSILONg; /* sec 46 #6 */
163 adj360(&LambdaSol);
164 l = 13.1763966 * days + lzero; /* sec 65 #4 */
165 adj360(&l);
166 Mm = l - (0.1114041 * days) - Pzero; /* sec 65 #5 */
167 adj360(&Mm);
168 Nm = Nzero - (0.0529539 * days); /* sec 65 #6 */
169 adj360(&Nm);
170 Ev = 1.2739 * sin(dtor(2*(l - LambdaSol) - Mm)); /* sec 65 #7 */
171 Ac = 0.1858 * sin(dtor(Msol)); /* sec 65 #8 */
172 A3 = 0.37 * sin(dtor(Msol));
173 Mmprime = Mm + Ev - Ac - A3; /* sec 65 #9 */
174 Ec = 6.2886 * sin(dtor(Mmprime)); /* sec 65 #10 */
175 A4 = 0.214 * sin(dtor(2 * Mmprime)); /* sec 65 #11 */
176 lprime = l + Ev + Ec - Ac + A4; /* sec 65 #12 */
177 V = 0.6583 * sin(dtor(2 * (lprime - LambdaSol))); /* sec 65 #13 */
178 ldprime = lprime + V; /* sec 65 #14 */
179 D = ldprime - LambdaSol; /* sec 67 #2 */
180 return(50.0 * (1 - cos(dtor(D)))); /* sec 67 #3 */
181 }
182
183 /*
184 * dtor --
185 * convert degrees to radians
186 */
187 static double
dtor(double deg)188 dtor(double deg)
189 {
190 return(deg * PI / 180);
191 }
192
193 /*
194 * adj360 --
195 * adjust value so 0 <= deg <= 360
196 */
197 static void
adj360(double * deg)198 adj360(double *deg)
199 {
200 for (;;)
201 if (*deg < 0)
202 *deg += 360;
203 else if (*deg > 360)
204 *deg -= 360;
205 else
206 break;
207 }
208
209 #define ATOI2(ar) ((ar)[0] - '0') * 10 + ((ar)[1] - '0'); (ar) += 2;
210 static time_t
parsetime(char * p)211 parsetime(char *p)
212 {
213 struct tm *lt;
214 int bigyear;
215 int yearset = 0;
216 time_t tval;
217 char *t;
218
219 for (t = p; *t; ++t) {
220 if (isdigit((unsigned char) *t))
221 continue;
222 badformat();
223 }
224
225 tval = time(NULL);
226 lt = localtime(&tval);
227 lt->tm_sec = 0;
228 lt->tm_min = 0;
229
230 switch (strlen(p)) {
231 case 10: /* yyyy */
232 bigyear = ATOI2(p);
233 lt->tm_year = bigyear * 100 - 1900;
234 yearset = 1;
235 /* FALLTHROUGH */
236 case 8: /* yy */
237 if (yearset) {
238 lt->tm_year += ATOI2(p);
239 } else {
240 lt->tm_year = ATOI2(p);
241 if (lt->tm_year < 69) /* hack for 2000 */
242 lt->tm_year += 100;
243 }
244 /* FALLTHROUGH */
245 case 6: /* mm */
246 lt->tm_mon = ATOI2(p);
247 if ((lt->tm_mon > 12) || !lt->tm_mon)
248 badformat();
249 --lt->tm_mon; /* time struct is 0 - 11 */
250 /* FALLTHROUGH */
251 case 4: /* dd */
252 lt->tm_mday = ATOI2(p);
253 if ((lt->tm_mday > 31) || !lt->tm_mday)
254 badformat();
255 /* FALLTHROUGH */
256 case 2: /* HH */
257 lt->tm_hour = ATOI2(p);
258 if (lt->tm_hour > 23)
259 badformat();
260 break;
261 default:
262 badformat();
263 }
264 /* The calling code needs a valid tm_ydays and this is the easiest
265 * way to get one */
266 if ((tval = mktime(lt)) == -1)
267 errx(1, "specified date is outside allowed range");
268 return (tval);
269 }
270
271 static void
badformat(void)272 badformat(void)
273 {
274 warnx("illegal time format");
275 (void)fprintf(stderr, "usage: %s [[[[[cc]yy]mm]dd]HH]\n",
276 getprogname());
277 exit(EXIT_FAILURE);
278 }
279