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