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