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 42static 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 87int main(int, char *[]); 88static void adj360(double *); 89static double dtor(double); 90static double potm(double); 91static time_t parsetime(char *); 92static void badformat(void) __dead; 93 94int 95main(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 */ 151static double 152potm(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 */ 187static double 188dtor(double deg) 189{ 190 return(deg * PI / 180); 191} 192 193/* 194 * adj360 -- 195 * adjust value so 0 <= deg <= 360 196 */ 197static void 198adj360(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; 210static time_t 211parsetime(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 271static void 272badformat(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