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