1 /* ---------------------------------------------------------------------------
2
3 moonphas.c
4
5 Notes:
6
7 This file contains routines for performing moon phase calculations.
8
9 The following suite of routines are to calculate the phase of the moon
10 for a given month, day, year. They compute the phase of the moon for
11 noon (UT) on the day requested, the start of the Julian day.
12
13 Revision history:
14
15 4.11.0
16 B.Marr 2007-12-15
17
18 Fix long-standing bug that only manifested itself in a DOS
19 build environment whereby the phases of the moon were being
20 erroneously detected on adjacent days.
21
22 Remove long-obsolete external 'moon file' concept. Now, we
23 depend solely on the algorithmic determination of moon phases.
24
25 Rename some variables, structures, and/or routines to be
26 clearer about their purpose and/or to allow easier searching
27 with fewer "false positives".
28
29 4.10.0
30 B.Marr 2006-07-19
31
32 Rename macro 'CALC_PHASE()' to 'CALC_MOON_PHASE()' to avoid
33 search collisions with routine 'calc_phase()'.
34
35 Reformatted comments and code to match my standards.
36
37 B.Marr 2006-07-12
38
39 Rename old 'getline()' routine to 'get_pcal_line()' to avoid a
40 compile-time namespace collision with the standard C library,
41 which seems to manifest itself only in Windows+Cygwin build
42 environments.
43
44 Fix longstanding bugs in the timezone offset ('-z' option)
45 calculations.
46
47 Use explicit units as part of variable name ('utc_offset'
48 becomes 'utc_offset_days') to avoid confusion.
49
50 Get rid of all the '#ifdef PROTOS' checks, which are pretty
51 much obsolete these days and just needlessly clutter up the
52 code.
53
54 Use common definition of PI already provided by 'math.h'.
55
56 4.8.0
57 B.Marr 2004-11-15
58
59 Remove Ctl-L (page eject) characters from source file. Remove
60 spaces embedded within tab fields.
61
62 4.7 AWR 02/24/1998 rename "template" (C++ reserved word)
63 as "templt" (cf. make_moonfile())
64
65 12/21/1997 clean up gcc warnings in -Wall mode
66
67 07/27/1997 delete obsolete FPR and PRT macros
68
69 4.5 AWR 11/06/1993 accept "opt -[AE]" in the moon file
70
71 06/25/1993 revise for pre-ANSI math libraries
72 without fmod()
73
74 04/28/1993 restructure function definitions so
75 function name appears in first column
76 (to facilitate searching for definition
77 by name)
78
79 04/22/1992 eliminated some unused variables and
80 calculations
81
82 4.4 AWR 01/20/1992 use alternate timezone spec (-z)
83
84 12/16/1991 Revise find_moonfile() for efficiency
85
86 4.3 AWR 12/05/1991 Search for moon file in Pcal's path
87 as well as in calendar file's path
88
89 10/25/1991 Many changes for support of moon
90 phase wildcards and -Z flag
91
92 4.11 AWR 08/23/1991 Revise is_quarter() to eliminate
93 occasional missing or duplicate
94 quarter-moons in "-m" mode; add
95 gen_phases()
96
97 4.1 AWR 08/02/1991 Fix bug in julday() where result is
98 calculated incorrectly if floor() has
99 no prototype
100
101 4.01 RLD 03/19/91 Upgraded calc_phase() to accurately
102 calculate the lunar phase for any
103 day without needing to resort to
104 a moon file.
105
106 4.0 AWR 03/07/1991 Add find_moonfile()
107
108 01/15/1991 Author: translated PostScript
109 routines to C and added moon
110 file routines
111 */
112
113 /* ---------------------------------------------------------------------------
114
115 Header Files
116
117 */
118
119 #include <stdio.h>
120 #include <string.h>
121 #include <ctype.h>
122 #include <math.h>
123
124 #include "pcaldefs.h"
125 #include "pcallang.h"
126 #include "protos.h"
127
128 /* ---------------------------------------------------------------------------
129
130 Type, Struct, & Enum Declarations
131
132 */
133
134 /* ---------------------------------------------------------------------------
135
136 Constant Declarations
137
138 */
139
140 /* Astronomical constants. */
141
142 #define epoch 2444238.5 /* 1980 January 0.0 */
143
144 /* Constants defining the Sun's apparent orbit. */
145
146 #define elonge 278.833540 /* ecliptic longitude of the Sun at epoch 1980.0 */
147 #define elongp 282.596403 /* ecliptic longitude of the Sun at perigee */
148 #define eccent 0.016718 /* eccentricity of Earth's orbit */
149
150 /* Elements of the Moon's orbit, epoch 1980.0. */
151
152 #define mmlong 64.975464 /* moon's mean lonigitude at the epoch */
153 #define mmlongp 349.383063 /* mean longitude of the perigee at the epoch */
154 #define mlnode 151.950429 /* mean longitude of the node at the epoch */
155 #define synmonth 29.53058868 /* synodic month (new Moon to new Moon) */
156
157 /* ---------------------------------------------------------------------------
158
159 Macro Definitions
160
161 */
162
163 /* Handy mathematical functions. */
164
165 #define sgn(x) (((x) < 0) ? -1 : ((x) > 0 ? 1 : 0)) /* extract sign */
166 #ifndef abs
167 #define abs(x) ((x) < 0 ? (-(x)) : (x)) /* absolute val */
168 #endif
169 #define fixangle(a) ((a) - 360.0 * (floor((a) / 360.0))) /* fix angle */
170 #define torad(d) ((d) * (M_PI / 180.0)) /* deg->rad */
171 #define todeg(d) ((d) * (180.0 / M_PI)) /* rad->deg */
172 #define FNITG(x) (sgn (x) * floor (abs (x)))
173
174 /* ---------------------------------------------------------------------------
175
176 Data Declarations (including externals)
177
178 */
179
180 /* ---------------------------------------------------------------------------
181
182 External Routine References & Function Prototypes
183
184 */
185
186 /* ---------------------------------------------------------------------------
187
188 * Routines to accurately calculate the phase of the moon
189 *
190 * Originally adapted from "moontool.c" by John Walker, Release 2.0.
191 *
192 * This routine (calc_phase) and its support routines were adapted from
193 * phase.c (v 1.2 88/08/26 22:29:42 jef) in the program "xphoon"
194 * (v 1.9 88/08/26 22:29:47 jef) by Jef Poskanzer and Craig Leres.
195 * The necessary notice follows...
196 *
197 * Copyright (C) 1988 by Jef Poskanzer and Craig Leres.
198 *
199 * Permission to use, copy, modify, and distribute this software and its
200 * documentation for any purpose and without fee is hereby granted,
201 * provided that the above copyright notice appear in all copies and that
202 * both that copyright notice and this permission notice appear in
203 * supporting documentation. This software is provided "as is" without
204 * express or implied warranty.
205 *
206 * These were added to "pcal" by RLD on 19-MAR-1991
207 */
208
209 /* ---------------------------------------------------------------------------
210
211 * julday -- calculate the julian date from input month, day, year
212 * N.B. - The Julian date is computed for noon UT.
213 *
214 * Adopted from Peter Duffett-Smith's book `Astronomy With Your
215 * Personal Computer' by Rick Dyson 18-MAR-1991
216 */
julday(int month,int day,int year)217 static double julday (int month, int day, int year)
218 {
219 int mn1, yr1;
220 double a, b, c, d, djd;
221
222 mn1 = month;
223 yr1 = year;
224 if ( yr1 < 0 ) yr1 = yr1 + 1;
225 if ( month < 3 ) {
226 mn1 = month + 12;
227 yr1 = yr1 - 1;
228 }
229 if (( year < 1582 ) ||
230 ( year == 1582 && month < 10 ) ||
231 ( year == 1582 && month == 10 && day < 15.0 )) {
232 b = 0;
233 }
234 else {
235 a = floor (yr1 / 100.0);
236 b = 2 - a + floor (a / 4);
237 }
238
239 if ( yr1 >= 0 ) c = floor (365.25 * yr1) - 694025;
240 else c = FNITG ((365.25 * yr1) - 0.75) - 694025;
241
242 d = floor (30.6001 * (mn1 + 1));
243 djd = b + c + d + day + 2415020.0;
244
245 return djd;
246 }
247
248 /* ---------------------------------------------------------------------------
249
250 * kepler - solve the equation of Kepler
251 */
kepler(double m,double ecc)252 static double kepler (double m, double ecc)
253 {
254 double e, delta;
255 #define EPSILON 1E-6
256
257 e = m = torad(m);
258 do {
259 delta = e - ecc * sin(e) - m;
260 e -= delta / (1 - ecc * cos(e));
261 } while (abs(delta) > EPSILON);
262 return e;
263 }
264
265 /* ---------------------------------------------------------------------------
266
267 * calc_phase - calculate phase of moon as a fraction:
268 *
269 * The argument is the time for which the phase is requested,
270 * expressed as the month, day and year. It returns the phase of
271 * the moon (0.0 -> 0.99) with the ordering as New Moon, First Quarter,
272 * Full Moon, and Last Quarter.
273 *
274 * Converted from the subroutine phase.c used by "xphoon.c" (see
275 * above disclaimer) into calc_phase() for use in "moonphas.c"
276 * by Rick Dyson 18-MAR-1991
277 */
calc_phase(int month,int inday,int year)278 double calc_phase (int month, int inday, int year)
279 {
280 double Day, N, M, Ec, Lambdasun, ml, MM;
281 double Ev, Ae, A3, MmP, mEc, A4, lP, V, lPP, MoonAge, pdate, moon_phase;
282 static int first = TRUE;
283 static double utc_offset_days;
284
285 /* Get the UTC offset on the first pass.
286
287 The original code used to normalize the UTC offset to +/- 12 hours. But
288 it was bug-ridden and also failed to take into account that some parts
289 of the world have offsets from UTC greater than 12 hours! Therefore,
290 beginning with v2.0.0, we don't attempt to normalize the user-specified
291 UTC timezone offset at all.
292 */
293 if (first) {
294
295 utc_offset_days = atof(time_zone) / 24.0;
296
297 if (DEBUG(DEBUG_MOON)) {
298 fprintf(stderr, "time_zone='%s' utc_offset_days = %.5lf\n", time_zone, utc_offset_days);
299 }
300
301 first = FALSE;
302 }
303
304 /* need to convert month, day, year into a Julian pdate */
305 pdate = julday(month, inday, year) + utc_offset_days;
306
307 /* Calculation of the Sun's position. */
308
309 Day = pdate - epoch; /* date within epoch */
310 N = fixangle((360 / 365.2422) * Day); /* mean anomaly of the Sun */
311 M = fixangle(N + elonge - elongp); /* convert from perigee
312 co-ordinates to epoch 1980.0 */
313 Ec = kepler(M, eccent); /* solve equation of Kepler */
314 Ec = sqrt((1 + eccent) / (1 - eccent)) * tan(Ec / 2);
315 Ec = 2 * todeg(atan(Ec)); /* true anomaly */
316 Lambdasun = fixangle(Ec + elongp); /* Sun's geocentric ecliptic longitude */
317
318 /* Calculation of the Moon's position. */
319
320 /* Moon's mean longitude. */
321 ml = fixangle(13.1763966 * Day + mmlong);
322
323 /* Moon's mean anomaly. */
324 MM = fixangle(ml - 0.1114041 * Day - mmlongp);
325
326 /* Moon's ascending node mean longitude. */
327 /* Not used -- commented out. */
328 /* MN = fixangle(mlnode - 0.0529539 * Day); */
329
330 /* Evection. */
331 Ev = 1.2739 * sin(torad(2 * (ml - Lambdasun) - MM));
332
333 /* Annual equation. */
334 Ae = 0.1858 * sin(torad(M));
335
336 /* Correction term. */
337 A3 = 0.37 * sin(torad(M));
338
339 /* Corrected anomaly. */
340 MmP = MM + Ev - Ae - A3;
341
342 /* Correction for the equation of the centre. */
343 mEc = 6.2886 * sin(torad(MmP));
344
345 /* Another correction term. */
346 A4 = 0.214 * sin(torad(2 * MmP));
347
348 /* Corrected longitude. */
349 lP = ml + Ev + mEc - Ae + A4;
350
351 /* Variation. */
352 V = 0.6583 * sin(torad(2 * (lP - Lambdasun)));
353
354 /* True longitude. */
355 lPP = lP + V;
356
357 /* Calculation of the phase of the Moon. */
358
359 /* Age of the Moon in degrees. */
360 MoonAge = lPP - Lambdasun;
361
362 moon_phase = fixangle(MoonAge) / 360.0;
363 if (moon_phase < 0.0) moon_phase += 1.0;
364
365 /* fprintf(stderr, "Moon phase on %04d-%02d-%02d: %.5lf\n", year, month, inday, moon_phase); */
366
367 return (moon_phase);
368 }
369
370 /* ---------------------------------------------------------------------------
371
372 * is_quarter - if current phase of moon coincides with quarter moon, return
373 * MOON_NM, MOON_1Q, etc.; otherwise return MOON_OTHER;
374 *
375 */
is_quarter(double prev,double curr,double next)376 static int is_quarter (double prev, double curr, double next)
377 {
378 int quarter;
379 double phase, diff;
380
381 /* adjust phases for 1 -> 0 wraparound */
382 if (curr < prev) curr++;
383 if (next < prev) next++;
384
385 /* if a quarter moon has occurred between "prev" and "next", return TRUE if
386 * "curr" is closer to it than "prev" or "next" is
387 */
388 for (quarter = 1; quarter <= 4; quarter++) {
389 if (prev < (phase = quarter/4.0) && next > phase &&
390 (diff = fabs(curr - phase)) < phase - prev && diff < next - phase) {
391 return quarter % 4; /* MOON_NM == 0 */
392 }
393 }
394
395 return MOON_OTHER;
396 }
397
398 /* ---------------------------------------------------------------------------
399
400 * gen_phases - fill array with moon phases for all days in month/year (plus
401 * extra entries at beginning and end for last day of previous month and
402 * first day of next month, respectively)
403 */
gen_phases(double phase[],int month,int year)404 static void gen_phases (double phase[], int month, int year)
405 {
406 int day, len;
407 date_str date;
408
409 /* start with moon phase for last day of previous month */
410 MAKE_DATE(date, month, 0, year);
411 normalize(&date);
412 phase[0] = calc_phase(date.mm, date.dd, date.yy);
413
414 /* add the moon phases for all days in the current month */
415 for (day = 1, len = LENGTH_OF(month, year); day <= len; day++) {
416 phase[day] = calc_phase(month, day, year);
417 }
418
419 /* append the moon phase for the first day of next month */
420 MAKE_DATE(date, month, len + 1, year);
421 normalize(&date);
422 phase[len + 1] = calc_phase(date.mm, date.dd, date.yy);
423 }
424
425 /* ---------------------------------------------------------------------------
426
427 * find_phase - calculate phase of moon using calc_phase() above. Sets
428 * *pquarter to MOON_NM, MOON_1Q, etc. if quarter moon, MOON_OTHER if not
429 */
find_phase(int month,int day,int year,int * pquarter)430 double find_phase (int month, int day, int year, int *pquarter)
431 {
432 static int sv_year = 0, sv_month = 0;
433 static double mphase[33]; /* 31 days + 2 dummies */
434 double phase;
435
436 /* calculate moon phase */
437
438 /* new month? fill mphase[] with moon phases */
439 if (month != sv_month || year != sv_year) {
440 gen_phases(mphase, month, year);
441 sv_month = month;
442 sv_year = year;
443 }
444
445 phase = mphase[day];
446 *pquarter = is_quarter(mphase[day-1], phase, mphase[day+1]);
447
448 return phase;
449 }
450