2 * Copyright (c) 1989, 1993
3 * The Regents of the University of California. All rights reserved.
5 * This code is derived from software posted to USENET.
7 * Redistribution and use in source and binary forms, with or without
8 * modification, are permitted provided that the following conditions
10 * 1. Redistributions of source code must retain the above copyright
11 * notice, this list of conditions and the following disclaimer.
12 * 2. Redistributions in binary form must reproduce the above copyright
13 * notice, this list of conditions and the following disclaimer in the
14 * documentation and/or other materials provided with the distribution.
15 * 4. Neither the name of the University nor the names of its contributors
16 * may be used to endorse or promote products derived from this software
17 * without specific prior written permission.
19 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
20 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
21 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
22 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
23 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
24 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
25 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
26 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
27 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
28 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
34 static const char copyright[] =
35 "@(#) Copyright (c) 1989, 1993\n\
36 The Regents of the University of California. All rights reserved.\n";
40 static const char sccsid[] = "@(#)pom.c 8.1 (Berkeley) 5/31/93";
43 #include <sys/cdefs.h>
44 __FBSDID("$FreeBSD$");
47 * Phase of the Moon. Calculates the current phase of the moon.
48 * Based on routines from `Practical Astronomy with Your Calculator',
49 * by Duffett-Smith. Comments give the section from the book that
50 * particular piece of code was adapted from.
52 * -- Keith E. Brandt VIII 1984
67 #define PI 3.14159265358979323846
70 #define EPSILONg 279.611371 /* solar ecliptic long at EPOCH */
71 #define RHOg 282.680403 /* solar ecliptic long of perigee at EPOCH */
72 #define ECCEN 0.01671542 /* solar orbit eccentricity */
73 #define lzero 18.251907 /* lunar mean long at EPOCH */
74 #define Pzero 192.917585 /* lunar mean long of perigee at EPOCH */
75 #define Nzero 55.204723 /* lunar mean long of node at EPOCH */
76 #define isleap(y) ((((y) % 4) == 0 && ((y) % 100) != 0) || ((y) % 400) == 0)
78 static void adj360(double *);
79 static double dtor(double);
80 static double potm(double onday);
81 static double potm_minute(double onday, int olddir);
84 pom(int year, double utcoffset, int *fms, int *nms)
86 double ffms[MAXMOONS];
87 double fnms[MAXMOONS];
90 fpom(year, utcoffset, ffms, fnms);
93 for (i = 0; ffms[i] != 0; i++)
94 fms[j++] = round(ffms[i]);
96 for (i = 0; fnms[i] != 0; i++)
97 nms[i] = round(fnms[i]);
102 fpom(int year, double utcoffset, double *ffms, double *fnms)
105 struct tm GMT, tmd_today, tmd_tomorrow;
106 double days_today, days_tomorrow, today, tomorrow;
110 double *pfnms, *pffms, t;
116 * We take the phase of the moon one second before and one second
119 memset(&tmd_today, 0, sizeof(tmd_today));
120 tmd_today.tm_year = year - 1900;
121 tmd_today.tm_mon = 0;
122 tmd_today.tm_mday = -1; /* 31 December */
123 tmd_today.tm_hour = 23;
124 tmd_today.tm_min = 59;
125 tmd_today.tm_sec = 59;
126 memset(&tmd_tomorrow, 0, sizeof(tmd_tomorrow));
127 tmd_tomorrow.tm_year = year - 1900;
128 tmd_tomorrow.tm_mon = 0;
129 tmd_tomorrow.tm_mday = 0; /* 01 January */
130 tmd_tomorrow.tm_hour = 0;
131 tmd_tomorrow.tm_min = 0;
132 tmd_tomorrow.tm_sec = 1;
134 tt = mktime(&tmd_today);
137 for (cnt = EPOCH; cnt < GMT.tm_year; ++cnt)
138 yeardays += isleap(1900 + cnt) ? DAYSPERLEAPYEAR : DAYSPERYEAR;
139 days_today = (GMT.tm_yday + 1) + ((GMT.tm_hour +
140 (GMT.tm_min / FSECSPERMINUTE) + (GMT.tm_sec / FSECSPERHOUR)) /
142 days_today += yeardays;
144 tt = mktime(&tmd_tomorrow);
147 for (cnt = EPOCH; cnt < GMT.tm_year; ++cnt)
148 yeardays += isleap(1900 + cnt) ? DAYSPERLEAPYEAR : DAYSPERYEAR;
149 days_tomorrow = (GMT.tm_yday + 1) + ((GMT.tm_hour +
150 (GMT.tm_min / FSECSPERMINUTE) + (GMT.tm_sec / FSECSPERHOUR)) /
152 days_tomorrow += yeardays;
154 today = potm(days_today); /* 30 December 23:59:59 */
155 tomorrow = potm(days_tomorrow); /* 31 December 00:00:01 */
156 olddir = today > tomorrow ? -1 : +1;
158 yeardays = 1 + isleap(year) ? DAYSPERLEAPYEAR : DAYSPERYEAR; /* reuse */
159 for (d = 0; d <= yeardays; d++) {
160 today = potm(days_today);
161 tomorrow = potm(days_tomorrow);
162 newdir = today > tomorrow ? -1 : +1;
163 if (olddir != newdir) {
164 t = potm_minute(days_today - 1, olddir) +
165 utcoffset / FHOURSPERDAY;
166 if (olddir == -1 && newdir == +1) {
169 } else if (olddir == +1 && newdir == -1) {
183 potm_minute(double onday, int olddir) {
184 double period = FSECSPERDAY / 2.0;
186 double before, after;
189 // printf("---> days:%g olddir:%d\n", days, olddir);
191 p1 = onday + (period / SECSPERDAY);
194 while (period > 30) { /* half a minute */
195 // printf("period:%g - p1:%g - ", period, p1);
196 p2 = p1 + (2.0 / SECSPERDAY);
199 // printf("before:%10.10g - after:%10.10g\n", before, after);
200 newdir = before < after ? -1 : +1;
201 if (olddir != newdir)
202 p1 += (period / SECSPERDAY);
204 p1 -= (period / SECSPERDAY);
206 // printf("newdir:%d - p1:%10.10f - period:%g\n",
207 // newdir, p1, period);
216 * return phase of the moon, as a percentage [0 ... 100]
221 double N, Msol, Ec, LambdaSol, l, Mm, Ev, Ac, A3, Mmprime;
222 double A4, lprime, V, ldprime, D, Nm;
224 N = 360 * onday / 365.2422; /* sec 42 #3 */
226 Msol = N + EPSILONg - RHOg; /* sec 42 #4 */
228 Ec = 360 / PI * ECCEN * sin(dtor(Msol)); /* sec 42 #5 */
229 LambdaSol = N + Ec + EPSILONg; /* sec 42 #6 */
231 l = 13.1763966 * onday + lzero; /* sec 61 #4 */
233 Mm = l - (0.1114041 * onday) - Pzero; /* sec 61 #5 */
235 Nm = Nzero - (0.0529539 * onday); /* sec 61 #6 */
237 Ev = 1.2739 * sin(dtor(2*(l - LambdaSol) - Mm)); /* sec 61 #7 */
238 Ac = 0.1858 * sin(dtor(Msol)); /* sec 61 #8 */
239 A3 = 0.37 * sin(dtor(Msol));
240 Mmprime = Mm + Ev - Ac - A3; /* sec 61 #9 */
241 Ec = 6.2886 * sin(dtor(Mmprime)); /* sec 61 #10 */
242 A4 = 0.214 * sin(dtor(2 * Mmprime)); /* sec 61 #11 */
243 lprime = l + Ev + Ec - Ac + A4; /* sec 61 #12 */
244 V = 0.6583 * sin(dtor(2 * (lprime - LambdaSol))); /* sec 61 #13 */
245 ldprime = lprime + V; /* sec 61 #14 */
246 D = ldprime - LambdaSol; /* sec 63 #2 */
247 return(50 * (1 - cos(dtor(D)))); /* sec 63 #3 */
252 * convert degrees to radians
258 return(deg * PI / 180);
263 * adjust value so 0 <= deg <= 360