1 /* 2 * Copyright (c) 1989 The Regents of the University of California. 3 * All rights reserved. 4 * 5 * This code is derived from software posted to USENET. 6 * 7 * Redistribution and use in source and binary forms are permitted 8 * provided that the above copyright notice and this paragraph are 9 * duplicated in all such forms and that any documentation, 10 * advertising materials, and other materials related to such 11 * distribution and use acknowledge that the software was developed 12 * by the University of California, Berkeley. The name of the 13 * University may not be used to endorse or promote products derived 14 * from this software without specific prior written permission. 15 * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR 16 * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED 17 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. 18 */ 19 20 #ifndef lint 21 char copyright[] = 22 "@(#) Copyright (c) 1989 The Regents of the University of California.\n\ 23 All rights reserved.\n"; 24 #endif /* not lint */ 25 26 #ifndef lint 27 static char sccsid[] = "@(#)pom.c 5.1 (Berkeley) 09/07/89"; 28 #endif /* not lint */ 29 30 /* 31 * Phase of the Moon. Calculates the current phase of the moon. 32 * Based on routines from `Practical Astronomy with Your Calculator', 33 * by Duffett-Smith. Comments give the section from the book that 34 * particular piece of code was adapted from. 35 * 36 * -- Keith E. Brandt VIII 1984 37 * 38 */ 39 40 #include <sys/time.h> 41 #include <stdio.h> 42 #include <tzfile.h> 43 #include <math.h> 44 45 #define PI 3.141592654 46 #define EPOCH 85 47 #define EPSILONg 279.611371 /* solar ecliptic long at EPOCH */ 48 #define RHOg 282.680403 /* solar ecliptic long of perigee at EPOCH */ 49 #define ECCEN 0.01671542 /* solar orbit eccentricity */ 50 #define lzero 18.251907 /* lunar mean long at EPOCH */ 51 #define Pzero 192.917585 /* lunar mean long of perigee at EPOCH */ 52 #define Nzero 55.204723 /* lunar mean long of node at EPOCH */ 53 54 main() 55 { 56 extern int errno; 57 struct timeval tp; 58 struct timezone tzp; 59 struct tm *GMT, *gmtime(); 60 double days, today, tomorrow, dtor(), adj360(), potm(); 61 int cnt; 62 char *strerror(); 63 64 if (gettimeofday(&tp,&tzp)) { 65 (void)fprintf(stderr, "pom: %s\n", strerror(errno)); 66 exit(1); 67 } 68 GMT = gmtime(&tp.tv_sec); 69 days = (GMT->tm_yday + 1) + ((GMT->tm_hour + 70 (GMT->tm_min / 60.0) + (GMT->tm_sec / 3600.0)) / 24.0); 71 for (cnt = EPOCH; cnt < GMT->tm_year; ++cnt) 72 days += isleap(cnt) ? 366 : 365; 73 today = potm(days) + .5; 74 (void)printf("The Moon is "); 75 if ((int)today == 100) 76 (void)printf("Full\n"); 77 else if (!(int)today) 78 (void)printf("New\n"); 79 else { 80 tomorrow = potm(days + 1); 81 if ((int)today == 50) 82 (void)printf("%s\n", tomorrow > today ? 83 "at the First Quarter" : "at the Last Quarter"); 84 else { 85 (void)printf("%s ", tomorrow > today ? 86 "Waxing" : "Waning"); 87 if (today > 50) 88 (void)printf("Gibbous (%1.0f%% of Full)\n", 89 today); 90 else if (today < 50) 91 (void)printf("Crescent (%1.0f%% of Full)\n", 92 today); 93 } 94 } 95 } 96 97 /* 98 * potm -- 99 * return phase of the moon 100 */ 101 double 102 potm(days) 103 double days; 104 { 105 double N, Msol, Ec, LambdaSol, l, Mm, Ev, Ac, A3, Mmprime; 106 double A4, lprime, V, ldprime, D, Nm; 107 108 N = 360 * days / 365.2422; /* sec 42 #3 */ 109 adj360(&N); 110 Msol = N + EPSILONg - RHOg; /* sec 42 #4 */ 111 adj360(&Msol); 112 Ec = 360 / PI * ECCEN * sin(dtor(Msol)); /* sec 42 #5 */ 113 LambdaSol = N + Ec + EPSILONg; /* sec 42 #6 */ 114 adj360(&LambdaSol); 115 l = 13.1763966 * days + lzero; /* sec 61 #4 */ 116 adj360(&l); 117 Mm = l - (0.1114041 * days) - Pzero; /* sec 61 #5 */ 118 adj360(&Mm); 119 Nm = Nzero - (0.0529539 * days); /* sec 61 #6 */ 120 adj360(&Nm); 121 Ev = 1.2739 * sin(dtor(2*(l - LambdaSol) - Mm)); /* sec 61 #7 */ 122 Ac = 0.1858 * sin(dtor(Msol)); /* sec 61 #8 */ 123 A3 = 0.37 * sin(dtor(Msol)); 124 Mmprime = Mm + Ev - Ac - A3; /* sec 61 #9 */ 125 Ec = 6.2886 * sin(dtor(Mmprime)); /* sec 61 #10 */ 126 A4 = 0.214 * sin(dtor(2 * Mmprime)); /* sec 61 #11 */ 127 lprime = l + Ev + Ec - Ac + A4; /* sec 61 #12 */ 128 V = 0.6583 * sin(dtor(2 * (lprime - LambdaSol))); /* sec 61 #13 */ 129 ldprime = lprime + V; /* sec 61 #14 */ 130 D = ldprime - LambdaSol; /* sec 63 #2 */ 131 return(50 * (1 - cos(dtor(D)))); /* sec 63 #3 */ 132 } 133 134 /* 135 * dtor -- 136 * convert degrees to radians 137 */ 138 double 139 dtor(deg) 140 double deg; 141 { 142 return(deg * PI / 180); 143 } 144 145 /* 146 * adj360 -- 147 * adjust value so 0 <= deg <= 360 148 */ 149 double 150 adj360(deg) 151 double *deg; 152 { 153 for (;;) 154 if (*deg < 0) 155 *deg += 360; 156 else if (*deg > 360) 157 *deg -= 360; 158 else 159 break; 160 } 161