Home | History | Annotate | Line # | Download | only in pom
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