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