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