pom.c revision 203932
1214501Srpaulo/* 2214501Srpaulo * Copyright (c) 1989, 1993 3214501Srpaulo * The Regents of the University of California. All rights reserved. 4214501Srpaulo * 5252726Srpaulo * This code is derived from software posted to USENET. 6252726Srpaulo * 7214501Srpaulo * Redistribution and use in source and binary forms, with or without 8214501Srpaulo * modification, are permitted provided that the following conditions 9214501Srpaulo * are met: 10214501Srpaulo * 1. Redistributions of source code must retain the above copyright 11214501Srpaulo * notice, this list of conditions and the following disclaimer. 12214501Srpaulo * 2. Redistributions in binary form must reproduce the above copyright 13214501Srpaulo * notice, this list of conditions and the following disclaimer in the 14214501Srpaulo * documentation and/or other materials provided with the distribution. 15214501Srpaulo * 3. Neither the name of the University nor the names of its contributors 16214501Srpaulo * may be used to endorse or promote products derived from this software 17214501Srpaulo * without specific prior written permission. 18214501Srpaulo * 19214501Srpaulo * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND 20214501Srpaulo * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 21214501Srpaulo * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 22214501Srpaulo * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE 23214501Srpaulo * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 24214501Srpaulo * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 25252726Srpaulo * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 26252726Srpaulo * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 27252726Srpaulo * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 28214501Srpaulo * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 29214501Srpaulo * SUCH DAMAGE. 30214501Srpaulo */ 31214501Srpaulo 32214501Srpaulo#if 0 33214501Srpaulo#ifndef lint 34214501Srpaulostatic const char copyright[] = 35214501Srpaulo"@(#) Copyright (c) 1989, 1993\n\ 36214501Srpaulo The Regents of the University of California. All rights reserved.\n"; 37214501Srpaulo#endif /* not lint */ 38214501Srpaulo 39214501Srpaulo#ifndef lint 40214501Srpaulostatic const char sccsid[] = "@(#)pom.c 8.1 (Berkeley) 5/31/93"; 41214501Srpaulo#endif /* not lint */ 42214501Srpaulo#endif 43214501Srpaulo#include <sys/cdefs.h> 44214501Srpaulo__FBSDID("$FreeBSD: head/games/pom/pom.c 203932 2010-02-15 18:46:02Z imp $"); 45214501Srpaulo 46214501Srpaulo/* 47214501Srpaulo * Phase of the Moon. Calculates the current phase of the moon. 48214501Srpaulo * Based on routines from `Practical Astronomy with Your Calculator', 49214501Srpaulo * by Duffett-Smith. Comments give the section from the book that 50214501Srpaulo * particular piece of code was adapted from. 51214501Srpaulo * 52214501Srpaulo * -- Keith E. Brandt VIII 1984 53214501Srpaulo * 54214501Srpaulo */ 55214501Srpaulo 56214501Srpaulo#include <stdio.h> 57214501Srpaulo#include <stdlib.h> 58214501Srpaulo#include <math.h> 59214501Srpaulo#include <string.h> 60214501Srpaulo#include <sysexits.h> 61214501Srpaulo#include <time.h> 62214501Srpaulo#include <unistd.h> 63214501Srpaulo 64214501Srpaulo#ifndef PI 65214501Srpaulo#define PI 3.14159265358979323846 66214501Srpaulo#endif 67214501Srpaulo#define EPOCH 85 68214501Srpaulo#define EPSILONg 279.611371 /* solar ecliptic long at EPOCH */ 69214501Srpaulo#define RHOg 282.680403 /* solar ecliptic long of perigee at EPOCH */ 70#define ECCEN 0.01671542 /* solar orbit eccentricity */ 71#define lzero 18.251907 /* lunar mean long at EPOCH */ 72#define Pzero 192.917585 /* lunar mean long of perigee at EPOCH */ 73#define Nzero 55.204723 /* lunar mean long of node at EPOCH */ 74#define isleap(y) ((((y) % 4) == 0 && ((y) % 100) != 0) || ((y) % 400) == 0) 75 76static void adj360(double *); 77static double dtor(double); 78static double potm(double); 79static void usage(char *progname); 80 81int 82main(int argc, char **argv) 83{ 84 time_t tt; 85 struct tm GMT, tmd; 86 double days, today, tomorrow; 87 int ch, cnt; 88 char *odate = NULL, *otime = NULL; 89 90 while ((ch = getopt(argc, argv, "d:t:")) != -1) 91 switch (ch) { 92 case 'd': 93 odate = optarg; 94 break; 95 case 't': 96 otime = optarg; 97 break; 98 default: 99 usage(argv[0]); 100 } 101 102 argc -= optind; 103 argv += optind; 104 105 if (argc) 106 usage(argv[0]); 107 108 /* Adjust based on users preferences */ 109 time(&tt); 110 if (otime != NULL || odate != NULL) { 111 /* Save today in case -d isn't specified */ 112 localtime_r(&tt, &tmd); 113 114 if (odate != NULL) { 115 tmd.tm_year = strtol(odate, NULL, 10) - 1900; 116 tmd.tm_mon = strtol(odate + 5, NULL, 10) - 1; 117 tmd.tm_mday = strtol(odate + 8, NULL, 10); 118 /* Use midnight as the middle of the night */ 119 tmd.tm_hour = 0; 120 tmd.tm_min = 0; 121 tmd.tm_sec = 0; 122 } 123 if (otime != NULL) { 124 tmd.tm_hour = strtol(otime, NULL, 10); 125 tmd.tm_min = strtol(otime + 3, NULL, 10); 126 tmd.tm_sec = strtol(otime + 6, NULL, 10); 127 } 128 tt = mktime(&tmd); 129 } 130 131 gmtime_r(&tt, &GMT); 132 days = (GMT.tm_yday + 1) + ((GMT.tm_hour + 133 (GMT.tm_min / 60.0) + (GMT.tm_sec / 3600.0)) / 24.0); 134 for (cnt = EPOCH; cnt < GMT.tm_year; ++cnt) 135 days += isleap(1900 + cnt) ? 366 : 365; 136 today = potm(days) + .5; 137 (void)printf("The Moon is "); 138 if ((int)today == 100) 139 (void)printf("Full\n"); 140 else if (!(int)today) 141 (void)printf("New\n"); 142 else { 143 tomorrow = potm(days + 1); 144 if ((int)today == 50) 145 (void)printf("%s\n", tomorrow > today ? 146 "at the First Quarter" : "at the Last Quarter"); 147 else { 148 (void)printf("%s ", tomorrow > today ? 149 "Waxing" : "Waning"); 150 if (today > 50) 151 (void)printf("Gibbous (%1.0f%% of Full)\n", 152 today); 153 else if (today < 50) 154 (void)printf("Crescent (%1.0f%% of Full)\n", 155 today); 156 } 157 } 158 159 return 0; 160} 161 162/* 163 * potm -- 164 * return phase of the moon 165 */ 166static double 167potm(double days) 168{ 169 double N, Msol, Ec, LambdaSol, l, Mm, Ev, Ac, A3, Mmprime; 170 double A4, lprime, V, ldprime, D, Nm; 171 172 N = 360 * days / 365.2422; /* sec 42 #3 */ 173 adj360(&N); 174 Msol = N + EPSILONg - RHOg; /* sec 42 #4 */ 175 adj360(&Msol); 176 Ec = 360 / PI * ECCEN * sin(dtor(Msol)); /* sec 42 #5 */ 177 LambdaSol = N + Ec + EPSILONg; /* sec 42 #6 */ 178 adj360(&LambdaSol); 179 l = 13.1763966 * days + lzero; /* sec 61 #4 */ 180 adj360(&l); 181 Mm = l - (0.1114041 * days) - Pzero; /* sec 61 #5 */ 182 adj360(&Mm); 183 Nm = Nzero - (0.0529539 * days); /* sec 61 #6 */ 184 adj360(&Nm); 185 Ev = 1.2739 * sin(dtor(2*(l - LambdaSol) - Mm)); /* sec 61 #7 */ 186 Ac = 0.1858 * sin(dtor(Msol)); /* sec 61 #8 */ 187 A3 = 0.37 * sin(dtor(Msol)); 188 Mmprime = Mm + Ev - Ac - A3; /* sec 61 #9 */ 189 Ec = 6.2886 * sin(dtor(Mmprime)); /* sec 61 #10 */ 190 A4 = 0.214 * sin(dtor(2 * Mmprime)); /* sec 61 #11 */ 191 lprime = l + Ev + Ec - Ac + A4; /* sec 61 #12 */ 192 V = 0.6583 * sin(dtor(2 * (lprime - LambdaSol))); /* sec 61 #13 */ 193 ldprime = lprime + V; /* sec 61 #14 */ 194 D = ldprime - LambdaSol; /* sec 63 #2 */ 195 return(50 * (1 - cos(dtor(D)))); /* sec 63 #3 */ 196} 197 198/* 199 * dtor -- 200 * convert degrees to radians 201 */ 202static double 203dtor(double deg) 204{ 205 206 return(deg * PI / 180); 207} 208 209/* 210 * adj360 -- 211 * adjust value so 0 <= deg <= 360 212 */ 213static void 214adj360(double *deg) 215{ 216 217 for (;;) 218 if (*deg < 0) 219 *deg += 360; 220 else if (*deg > 360) 221 *deg -= 360; 222 else 223 break; 224} 225 226static void 227usage(char *progname) 228{ 229 230 fprintf(stderr, "Usage: %s [-d yyyy.mm.dd] [-t hh:mm:ss]\n", progname); 231 exit(EX_USAGE); 232} 233