pom.c revision 1.10
1/* $OpenBSD: pom.c,v 1.10 2003/07/10 00:03:01 david 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#ifndef lint 36static char copyright[] = 37"@(#) Copyright (c) 1989, 1993\n\ 38 The Regents of the University of California. All rights reserved.\n"; 39#endif /* not lint */ 40 41#ifndef lint 42#if 0 43static char sccsid[] = "@(#)pom.c 8.1 (Berkeley) 5/31/93"; 44#else 45static char rcsid[] = "$OpenBSD: pom.c,v 1.10 2003/07/10 00:03:01 david Exp $"; 46#endif 47#endif /* not lint */ 48 49/* 50 * Phase of the Moon. Calculates the current phase of the moon. 51 * Based on routines from `Practical Astronomy with Your Calculator', 52 * by Duffett-Smith. Comments give the section from the book that 53 * particular piece of code was adapted from. 54 * 55 * -- Keith E. Brandt VIII 1984 56 * 57 * Updated to the Third Edition of Duffett-Smith's book, IX 1998 58 * 59 */ 60 61#include <sys/time.h> 62#include <sys/types.h> 63#include <ctype.h> 64#include <stdio.h> 65#include <stdlib.h> 66#include <string.h> 67#include <math.h> 68#include <err.h> 69#include <tzfile.h> 70#include <unistd.h> 71 72#ifndef M_PI 73#define M_PI 3.14159265358979323846 74#endif 75 76#define EPOCH 90 77#define EPSILONg 279.403303 /* solar ecliptic long at EPOCH */ 78#define RHOg 282.768422 /* solar ecliptic long of perigee at EPOCH */ 79#define ECCEN 0.016713 /* solar orbit eccentricity */ 80#define lzero 318.351648 /* lunar mean long at EPOCH */ 81#define Pzero 36.340410 /* lunar mean long of perigee at EPOCH */ 82#define Nzero 318.510107 /* lunar mean long of node at EPOCH */ 83 84void adj360(double *); 85double dtor(double); 86double potm(double); 87time_t parsetime(char *); 88void badformat(void); 89 90int 91main(argc, argv) 92 int argc; 93 char *argv[]; 94{ 95 struct timeval tp; 96 struct timezone tzp; 97 struct tm *GMT; 98 time_t tmpt; 99 double days, today, tomorrow; 100 int cnt; 101 char buf[1024]; 102 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 if (gettimeofday(&tp,&tzp)) 110 err(1, "gettimeofday"); 111 tmpt = tp.tv_sec; 112 } 113 GMT = gmtime(&tmpt); 114 days = (GMT->tm_yday + 1) + ((GMT->tm_hour + 115 (GMT->tm_min / 60.0) + (GMT->tm_sec / 3600.0)) / 24.0); 116 for (cnt = EPOCH; cnt < GMT->tm_year; ++cnt) 117 days += isleap(cnt + TM_YEAR_BASE) ? 366 : 365; 118 /* Selected time could be before EPOCH */ 119 for (cnt = GMT->tm_year; cnt < EPOCH; ++cnt) 120 days -= isleap(cnt + TM_YEAR_BASE) ? 366 : 365; 121 today = potm(days) + 0.5; 122 (void)printf("The Moon is "); 123 if ((int)today == 100) 124 (void)printf("Full\n"); 125 else if (!(int)today) 126 (void)printf("New\n"); 127 else { 128 tomorrow = potm(days + 1); 129 if ((int)today == 50) 130 (void)printf("%s\n", tomorrow > today ? 131 "at the First Quarter" : "at the Last Quarter"); 132 /* today is 0.5 too big, but it doesn't matter here 133 * since the phase is changing fast enough 134 */ 135 else { 136 today -= 0.5; /* Now it might matter */ 137 (void)printf("%s ", tomorrow > today ? 138 "Waxing" : "Waning"); 139 if (today > 50) 140 (void)printf("Gibbous (%1.0f%% of Full)\n", 141 today); 142 else if (today < 50) 143 (void)printf("Crescent (%1.0f%% of Full)\n", 144 today); 145 } 146 } 147 exit(0); 148} 149 150/* 151 * potm -- 152 * return phase of the moon 153 */ 154double 155potm(days) 156 double days; 157{ 158 double N, Msol, Ec, LambdaSol, l, Mm, Ev, Ac, A3, Mmprime; 159 double A4, lprime, V, ldprime, D, Nm; 160 161 N = 360.0 * days / 365.242191; /* sec 46 #3 */ 162 adj360(&N); 163 Msol = N + EPSILONg - RHOg; /* sec 46 #4 */ 164 adj360(&Msol); 165 Ec = 360 / M_PI * ECCEN * sin(dtor(Msol)); /* sec 46 #5 */ 166 LambdaSol = N + Ec + EPSILONg; /* sec 46 #6 */ 167 adj360(&LambdaSol); 168 l = 13.1763966 * days + lzero; /* sec 65 #4 */ 169 adj360(&l); 170 Mm = l - (0.1114041 * days) - Pzero; /* sec 65 #5 */ 171 adj360(&Mm); 172 Nm = Nzero - (0.0529539 * days); /* sec 65 #6 */ 173 adj360(&Nm); 174 Ev = 1.2739 * sin(dtor(2*(l - LambdaSol) - Mm)); /* sec 65 #7 */ 175 Ac = 0.1858 * sin(dtor(Msol)); /* sec 65 #8 */ 176 A3 = 0.37 * sin(dtor(Msol)); 177 Mmprime = Mm + Ev - Ac - A3; /* sec 65 #9 */ 178 Ec = 6.2886 * sin(dtor(Mmprime)); /* sec 65 #10 */ 179 A4 = 0.214 * sin(dtor(2 * Mmprime)); /* sec 65 #11 */ 180 lprime = l + Ev + Ec - Ac + A4; /* sec 65 #12 */ 181 V = 0.6583 * sin(dtor(2 * (lprime - LambdaSol))); /* sec 65 #13 */ 182 ldprime = lprime + V; /* sec 65 #14 */ 183 D = ldprime - LambdaSol; /* sec 67 #2 */ 184 return(50.0 * (1 - cos(dtor(D)))); /* sec 67 #3 */ 185} 186 187/* 188 * dtor -- 189 * convert degrees to radians 190 */ 191double 192dtor(deg) 193 double deg; 194{ 195 return(deg * M_PI / 180); 196} 197 198/* 199 * adj360 -- 200 * adjust value so 0 <= deg <= 360 201 */ 202void 203adj360(deg) 204 double *deg; 205{ 206 for (;;) 207 if (*deg < 0.0) 208 *deg += 360.0; 209 else if (*deg > 360.0) 210 *deg -= 360.0; 211 else 212 break; 213} 214 215#define ATOI2(ar) ((ar)[0] - '0') * 10 + ((ar)[1] - '0'); (ar) += 2; 216time_t 217parsetime(p) 218 char *p; 219{ 220 struct tm *lt; 221 int bigyear; 222 int yearset = 0; 223 time_t tval; 224 char *t; 225 226 for (t = p; *t; ++t) { 227 if (isdigit(*t)) 228 continue; 229 badformat(); 230 } 231 232 tval = time(NULL); 233 lt = localtime(&tval); 234 lt->tm_sec = 0; 235 lt->tm_min = 0; 236 237 switch (strlen(p)) { 238 case 10: /* yyyy */ 239 bigyear = ATOI2(p); 240 lt->tm_year = bigyear * 100 - TM_YEAR_BASE; 241 yearset = 1; 242 /* FALLTHROUGH */ 243 case 8: /* yy */ 244 if (yearset) { 245 lt->tm_year += ATOI2(p); 246 } else { 247 lt->tm_year = ATOI2(p) + 1900 - TM_YEAR_BASE; 248 if (lt->tm_year < 69) /* hack for 2000 */ 249 lt->tm_year += 100; 250 } 251 /* FALLTHROUGH */ 252 case 6: /* mm */ 253 lt->tm_mon = ATOI2(p); 254 if ((lt->tm_mon > 12) || !lt->tm_mon) 255 badformat(); 256 --lt->tm_mon; /* time struct is 0 - 11 */ 257 /* FALLTHROUGH */ 258 case 4: /* dd */ 259 lt->tm_mday = ATOI2(p); 260 if ((lt->tm_mday > 31) || !lt->tm_mday) 261 badformat(); 262 /* FALLTHROUGH */ 263 case 2: /* HH */ 264 lt->tm_hour = ATOI2(p); 265 if (lt->tm_hour > 23) 266 badformat(); 267 break; 268 default: 269 badformat(); 270 } 271 /* The calling code needs a valid tm_ydays and this is the easiest 272 * way to get one */ 273 if ((tval = mktime(lt)) == -1) 274 errx(1, "specified date is outside allowed range"); 275 return (tval); 276} 277 278void 279badformat() 280{ 281 warnx("illegal time format"); 282 (void)fprintf(stderr, "usage: pom [[[[[cc]yy]mm]dd]HH]\n"); 283 exit(1); 284} 285