1 /*-
2 * SPDX-License-Identifier: BSD-3-Clause
3 *
4 * Copyright (c) 1989, 1993
5 * The Regents of the University of California. All rights reserved.
6 *
7 * This code is derived from software posted to USENET.
8 *
9 * Redistribution and use in source and binary forms, with or without
10 * modification, are permitted provided that the following conditions
11 * are met:
12 * 1. Redistributions of source code must retain the above copyright
13 * notice, this list of conditions and the following disclaimer.
14 * 2. Redistributions in binary form must reproduce the above copyright
15 * notice, this list of conditions and the following disclaimer in the
16 * documentation and/or other materials provided with the distribution.
17 * 3. Neither the name of the University nor the names of its contributors
18 * may be used to endorse or promote products derived from this software
19 * without specific prior written permission.
20 *
21 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
22 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
23 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
24 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
25 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
26 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
27 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
28 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
29 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
30 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
31 * SUCH DAMAGE.
32 */
33
34 /*
35 * Phase of the Moon. Calculates the current phase of the moon.
36 * Based on routines from `Practical Astronomy with Your Calculator',
37 * by Duffett-Smith. Comments give the section from the book that
38 * particular piece of code was adapted from.
39 *
40 * -- Keith E. Brandt VIII 1984
41 *
42 */
43
44 #include <stdio.h>
45 #include <stdlib.h>
46 #include <math.h>
47 #include <string.h>
48 #include <sysexits.h>
49 #include <time.h>
50 #include <unistd.h>
51
52 #include "calendar.h"
53
54 #ifndef PI
55 #define PI 3.14159265358979323846
56 #endif
57 #define EPOCH 85
58 #define EPSILONg 279.611371 /* solar ecliptic long at EPOCH */
59 #define RHOg 282.680403 /* solar ecliptic long of perigee at EPOCH */
60 #define ECCEN 0.01671542 /* solar orbit eccentricity */
61 #define lzero 18.251907 /* lunar mean long at EPOCH */
62 #define Pzero 192.917585 /* lunar mean long of perigee at EPOCH */
63 #define Nzero 55.204723 /* lunar mean long of node at EPOCH */
64 #define isleap(y) ((((y) % 4) == 0 && ((y) % 100) != 0) || ((y) % 400) == 0)
65
66 static void adj360(double *);
67 static double dtor(double);
68 static double potm(double onday);
69 static double potm_minute(double onday, int olddir);
70
71 void
pom(int year,double utcoffset,int * fms,int * nms)72 pom(int year, double utcoffset, int *fms, int *nms)
73 {
74 double ffms[MAXMOONS];
75 double fnms[MAXMOONS];
76 int i, j;
77
78 fpom(year, utcoffset, ffms, fnms);
79
80 j = 0;
81 for (i = 0; ffms[i] != 0; i++)
82 fms[j++] = round(ffms[i]);
83 fms[i] = -1;
84 for (i = 0; fnms[i] != 0; i++)
85 nms[i] = round(fnms[i]);
86 nms[i] = -1;
87 }
88
89 void
fpom(int year,double utcoffset,double * ffms,double * fnms)90 fpom(int year, double utcoffset, double *ffms, double *fnms)
91 {
92 time_t tt;
93 struct tm GMT, tmd_today, tmd_tomorrow;
94 double days_today, days_tomorrow, today, tomorrow;
95 int cnt, d;
96 int yeardays;
97 int olddir, newdir;
98 double *pfnms, *pffms, t;
99
100 pfnms = fnms;
101 pffms = ffms;
102
103 /*
104 * We take the phase of the moon one second before and one second
105 * after midnight.
106 */
107 memset(&tmd_today, 0, sizeof(tmd_today));
108 tmd_today.tm_year = year - 1900;
109 tmd_today.tm_mon = 0;
110 tmd_today.tm_mday = -1; /* 31 December */
111 tmd_today.tm_hour = 23;
112 tmd_today.tm_min = 59;
113 tmd_today.tm_sec = 59;
114 memset(&tmd_tomorrow, 0, sizeof(tmd_tomorrow));
115 tmd_tomorrow.tm_year = year - 1900;
116 tmd_tomorrow.tm_mon = 0;
117 tmd_tomorrow.tm_mday = 0; /* 01 January */
118 tmd_tomorrow.tm_hour = 0;
119 tmd_tomorrow.tm_min = 0;
120 tmd_tomorrow.tm_sec = 1;
121
122 tt = mktime(&tmd_today);
123 gmtime_r(&tt, &GMT);
124 yeardays = 0;
125 for (cnt = EPOCH; cnt < GMT.tm_year; ++cnt)
126 yeardays += isleap(1900 + cnt) ? DAYSPERLEAPYEAR : DAYSPERYEAR;
127 days_today = (GMT.tm_yday + 1) + ((GMT.tm_hour +
128 (GMT.tm_min / FSECSPERMINUTE) + (GMT.tm_sec / FSECSPERHOUR)) /
129 FHOURSPERDAY);
130 days_today += yeardays;
131
132 tt = mktime(&tmd_tomorrow);
133 gmtime_r(&tt, &GMT);
134 yeardays = 0;
135 for (cnt = EPOCH; cnt < GMT.tm_year; ++cnt)
136 yeardays += isleap(1900 + cnt) ? DAYSPERLEAPYEAR : DAYSPERYEAR;
137 days_tomorrow = (GMT.tm_yday + 1) + ((GMT.tm_hour +
138 (GMT.tm_min / FSECSPERMINUTE) + (GMT.tm_sec / FSECSPERHOUR)) /
139 FHOURSPERDAY);
140 days_tomorrow += yeardays;
141
142 today = potm(days_today); /* 30 December 23:59:59 */
143 tomorrow = potm(days_tomorrow); /* 31 December 00:00:01 */
144 olddir = today > tomorrow ? -1 : +1;
145
146 yeardays = 1 + (isleap(year) ? DAYSPERLEAPYEAR : DAYSPERYEAR); /* reuse */
147 for (d = 0; d <= yeardays; d++) {
148 today = potm(days_today);
149 tomorrow = potm(days_tomorrow);
150 newdir = today > tomorrow ? -1 : +1;
151 if (olddir != newdir) {
152 t = potm_minute(days_today - 1, olddir) +
153 utcoffset / FHOURSPERDAY;
154 if (olddir == -1 && newdir == +1) {
155 *pfnms = d - 1 + t;
156 pfnms++;
157 } else if (olddir == +1 && newdir == -1) {
158 *pffms = d - 1 + t;
159 pffms++;
160 }
161 }
162 olddir = newdir;
163 days_today++;
164 days_tomorrow++;
165 }
166 *pffms = -1;
167 *pfnms = -1;
168 }
169
170 static double
potm_minute(double onday,int olddir)171 potm_minute(double onday, int olddir) {
172 double period = FSECSPERDAY / 2.0;
173 double p1, p2;
174 double before, after;
175 int newdir;
176
177 // printf("---> days:%g olddir:%d\n", days, olddir);
178
179 p1 = onday + (period / SECSPERDAY);
180 period /= 2;
181
182 while (period > 30) { /* half a minute */
183 // printf("period:%g - p1:%g - ", period, p1);
184 p2 = p1 + (2.0 / SECSPERDAY);
185 before = potm(p1);
186 after = potm(p2);
187 // printf("before:%10.10g - after:%10.10g\n", before, after);
188 newdir = before < after ? -1 : +1;
189 if (olddir != newdir)
190 p1 += (period / SECSPERDAY);
191 else
192 p1 -= (period / SECSPERDAY);
193 period /= 2;
194 // printf("newdir:%d - p1:%10.10f - period:%g\n",
195 // newdir, p1, period);
196 }
197 p1 -= floor(p1);
198 //exit(0);
199 return (p1);
200 }
201
202 /*
203 * potm --
204 * return phase of the moon, as a percentage [0 ... 100]
205 */
206 static double
potm(double onday)207 potm(double onday)
208 {
209 double N, Msol, Ec, LambdaSol, l, Mm, Ev, Ac, A3, Mmprime;
210 double A4, lprime, V, ldprime, D, Nm;
211
212 N = 360 * onday / 365.2422; /* sec 42 #3 */
213 adj360(&N);
214 Msol = N + EPSILONg - RHOg; /* sec 42 #4 */
215 adj360(&Msol);
216 Ec = 360 / PI * ECCEN * sin(dtor(Msol)); /* sec 42 #5 */
217 LambdaSol = N + Ec + EPSILONg; /* sec 42 #6 */
218 adj360(&LambdaSol);
219 l = 13.1763966 * onday + lzero; /* sec 61 #4 */
220 adj360(&l);
221 Mm = l - (0.1114041 * onday) - Pzero; /* sec 61 #5 */
222 adj360(&Mm);
223 Nm = Nzero - (0.0529539 * onday); /* sec 61 #6 */
224 adj360(&Nm);
225 Ev = 1.2739 * sin(dtor(2*(l - LambdaSol) - Mm)); /* sec 61 #7 */
226 Ac = 0.1858 * sin(dtor(Msol)); /* sec 61 #8 */
227 A3 = 0.37 * sin(dtor(Msol));
228 Mmprime = Mm + Ev - Ac - A3; /* sec 61 #9 */
229 Ec = 6.2886 * sin(dtor(Mmprime)); /* sec 61 #10 */
230 A4 = 0.214 * sin(dtor(2 * Mmprime)); /* sec 61 #11 */
231 lprime = l + Ev + Ec - Ac + A4; /* sec 61 #12 */
232 V = 0.6583 * sin(dtor(2 * (lprime - LambdaSol))); /* sec 61 #13 */
233 ldprime = lprime + V; /* sec 61 #14 */
234 D = ldprime - LambdaSol; /* sec 63 #2 */
235 return(50 * (1 - cos(dtor(D)))); /* sec 63 #3 */
236 }
237
238 /*
239 * dtor --
240 * convert degrees to radians
241 */
242 static double
dtor(double deg)243 dtor(double deg)
244 {
245
246 return(deg * PI / 180);
247 }
248
249 /*
250 * adj360 --
251 * adjust value so 0 <= deg <= 360
252 */
253 static void
adj360(double * deg)254 adj360(double *deg)
255 {
256
257 for (;;)
258 if (*deg < 0)
259 *deg += 360;
260 else if (*deg > 360)
261 *deg -= 360;
262 else
263 break;
264 }
265