1 /*-
2 * SPDX-License-Identifier: BSD-3-Clause
3 *
4 * Copyright (c) 2019-2020 The DragonFly Project. All rights reserved.
5 *
6 * This code is derived from software contributed to The DragonFly Project
7 * by Aaron LI <aly@aaronly.me>
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 *
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
17 * the documentation and/or other materials provided with the
18 * distribution.
19 * 3. Neither the name of The DragonFly Project nor the names of its
20 * contributors may be used to endorse or promote products derived
21 * from this software without specific, prior written permission.
22 *
23 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
24 * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
25 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
26 * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
27 * COPYRIGHT HOLDERS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
28 * INCIDENTAL, SPECIAL, EXEMPLARY OR CONSEQUENTIAL DAMAGES (INCLUDING,
29 * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
30 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
31 * AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
32 * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
33 * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
34 * SUCH DAMAGE.
35 *
36 * Reference:
37 * Calendrical Calculations, The Ultimate Edition (4th Edition)
38 * by Edward M. Reingold and Nachum Dershowitz
39 * 2018, Cambridge University Press
40 */
41
42 #include <math.h>
43 #include <stdbool.h>
44 #include <stddef.h>
45 #include <stdio.h>
46 #include <stdlib.h>
47
48 #include "basics.h"
49 #include "gregorian.h"
50 #include "moon.h"
51 #include "sun.h"
52 #include "utils.h"
53
54
55 /*
56 * Time from new moon to new moon (i.e., a lunation).
57 * Ref: Sec.(14.6), Eq.(14.44)
58 */
59 const double mean_synodic_month = 29.530588861;
60
61
62 /*
63 * Argument data 1 used by 'nth_new_moon()'.
64 * Ref: Sec.(14.6), Table(14.3)
65 */
66 static const struct nth_new_moon_arg1 {
67 double v;
68 int w;
69 int x;
70 int y;
71 int z;
72 } nth_new_moon_data1[] = {
73 { -0.40720, 0, 0, 1, 0 },
74 { 0.17241, 1, 1, 0, 0 },
75 { 0.01608, 0, 0, 2, 0 },
76 { 0.01039, 0, 0, 0, 2 },
77 { 0.00739, 1, -1, 1, 0 },
78 { -0.00514, 1, 1, 1, 0 },
79 { 0.00208, 2, 2, 0, 0 },
80 { -0.00111, 0, 0, 1, -2 },
81 { -0.00057, 0, 0, 1, 2 },
82 { 0.00056, 1, 1, 2, 0 },
83 { -0.00042, 0, 0, 3, 0 },
84 { 0.00042, 1, 1, 0, 2 },
85 { 0.00038, 1, 1, 0, -2 },
86 { -0.00024, 1, -1, 2, 0 },
87 { -0.00007, 0, 2, 1, 0 },
88 { 0.00004, 0, 0, 2, -2 },
89 { 0.00004, 0, 3, 0, 0 },
90 { 0.00003, 0, 1, 1, -2 },
91 { 0.00003, 0, 0, 2, 2 },
92 { -0.00003, 0, 1, 1, 2 },
93 { 0.00003, 0, -1, 1, 2 },
94 { -0.00002, 0, -1, 1, -2 },
95 { -0.00002, 0, 1, 3, 0 },
96 { 0.00002, 0, 0, 4, 0 },
97 };
98
99 /*
100 * Argument data 2 used by 'nth_new_moon()'.
101 * Ref: Sec.(14.6), Table(14.4)
102 */
103 static const struct nth_new_moon_arg2 {
104 double i;
105 double j;
106 double l;
107 } nth_new_moon_data2[] = {
108 { 251.88, 0.016321, 0.000165 },
109 { 251.83, 26.651886, 0.000164 },
110 { 349.42, 36.412478, 0.000126 },
111 { 84.66, 18.206239, 0.000110 },
112 { 141.74, 53.303771, 0.000062 },
113 { 207.14, 2.453732, 0.000060 },
114 { 154.84, 7.306860, 0.000056 },
115 { 34.52, 27.261239, 0.000047 },
116 { 207.19, 0.121824, 0.000042 },
117 { 291.34, 1.844379, 0.000040 },
118 { 161.72, 24.198154, 0.000037 },
119 { 239.56, 25.513099, 0.000035 },
120 { 331.55, 3.592518, 0.000023 },
121 };
122
123 /*
124 * Calculate the moment of the $n-th new moon after the new moon of
125 * January 11, 1 (Gregorian), which is the first new moon after RD 0.
126 * This function is centered upon January, 2000.
127 * Ref: Sec.(14.6), Eq.(14.45)
128 */
129 double
nth_new_moon(int n)130 nth_new_moon(int n)
131 {
132 int n0 = 24724; /* Months from RD 0 until j2000 */
133 int k = n - n0; /* Months since j2000 */
134 double nm = 1236.85; /* Months per century */
135 double c = k / nm; /* Julian centuries */
136 double j2000 = 0.5 + gregorian_new_year(2000);
137
138 double coef_ap[] = { 5.09766, mean_synodic_month * nm,
139 0.00015437, -0.000000150, 0.00000000073 };
140 double approx = j2000 + poly(c, coef_ap, nitems(coef_ap));
141
142 double extra = 0.000325 * sin_deg(299.77 + 132.8475848 * c -
143 0.009173 * c*c);
144
145 double coef_sa[] = { 2.5534, 29.10535670 * nm,
146 -0.0000014, -0.00000011 };
147 double coef_la[] = { 201.5643, 385.81693528 * nm,
148 0.0107582, 0.00001238, -0.000000058 };
149 double coef_ma[] = { 160.7108, 390.67050284 * nm,
150 -0.0016118, -0.00000227, 0.000000011 };
151 double coef_om[] = { 124.7746, -1.56375588 * nm,
152 0.0020672, 0.00000215 };
153 double solar_anomaly = poly(c, coef_sa, nitems(coef_sa));
154 double lunar_anomaly = poly(c, coef_la, nitems(coef_la));
155 double moon_argument = poly(c, coef_ma, nitems(coef_ma));
156 double omega = poly(c, coef_om, nitems(coef_om));
157 double E = 1.0 - 0.002516 * c - 0.0000074 * c*c;
158
159 double sum_c = 0.0;
160 const struct nth_new_moon_arg1 *arg1;
161 for (size_t i = 0; i < nitems(nth_new_moon_data1); i++) {
162 arg1 = &nth_new_moon_data1[i];
163 sum_c += arg1->v * pow(E, arg1->w) * sin_deg(
164 arg1->x * solar_anomaly +
165 arg1->y * lunar_anomaly +
166 arg1->z * moon_argument);
167 }
168 double correction = -0.00017 * sin_deg(omega) + sum_c;
169
170 double additional = 0.0;
171 const struct nth_new_moon_arg2 *arg2;
172 for (size_t i = 0; i < nitems(nth_new_moon_data2); i++) {
173 arg2 = &nth_new_moon_data2[i];
174 additional += arg2->l * sin_deg(arg2->i + arg2->j * k);
175 }
176
177 double dt = approx + correction + extra + additional;
178 return universal_from_dynamical(dt);
179 }
180
181 /*
182 * Mean longitude of moon at moment given in Julian centuries $c.
183 * Ref: Sec.(14.6), Eq.(14.49)
184 */
185 static double
lunar_longitude_mean(double c)186 lunar_longitude_mean(double c)
187 {
188 double coef[] = { 218.3164477, 481267.88123421, -0.0015786,
189 1.0 / 538841.0, -1.0 / 65194000.0 };
190 return mod_f(poly(c, coef, nitems(coef)), 360);
191 }
192
193 /*
194 * Elongation of moon (angular distance from Sun) at moment given in Julian
195 * centuries $c.
196 * Ref: Sec.(14.6), Eq.(14.50)
197 */
198 static double
lunar_elongation(double c)199 lunar_elongation(double c)
200 {
201 double coef[] = { 297.8501921, 445267.1114034, -0.0018819,
202 1.0 / 545868.0, -1.0 / 113065000.0 };
203 return mod_f(poly(c, coef, nitems(coef)), 360);
204 }
205
206 /*
207 * Mean anomaly of Sun (angular distance from perihelion) at moment given in
208 * Julian centuries $c.
209 * Ref: Sec.(14.6), Eq.(14.51)
210 */
211 static double
solar_anomaly(double c)212 solar_anomaly(double c)
213 {
214 double coef[] = { 357.5291092, 35999.0502909, -0.0001536,
215 1.0 / 24490000.0 };
216 return mod_f(poly(c, coef, nitems(coef)), 360);
217 }
218
219 /*
220 * Mean anomaly of moon (angular distance from perigee) at moment given in
221 * Julian centuries $c.
222 * Ref: Sec.(14.6), Eq.(14.52)
223 */
224 static double
lunar_anomaly(double c)225 lunar_anomaly(double c)
226 {
227 double coef[] = { 134.9633964, 477198.8675055, 0.0087414,
228 1.0 / 69699.0, -1.0 / 14712000.0 };
229 return mod_f(poly(c, coef, nitems(coef)), 360);
230 }
231
232 /*
233 * Moon's argument of latitude (the distance from the moon's node) at moment
234 * given in Julian centuries $c.
235 * Ref: Sec.(14.6), Eq.(14.53)
236 */
237 static double
moon_node(double c)238 moon_node(double c)
239 {
240 double coef[] = { 93.2720950, 483202.0175233, -0.0036539,
241 -1.0 / 3526000.0, 1.0 / 863310000.0 };
242 return mod_f(poly(c, coef, nitems(coef)), 360);
243 }
244
245 /*
246 * Argument data used by 'lunar_longitude()'.
247 * Ref: Sec.(14.6), Table(14.5)
248 */
249 static const struct lunar_longitude_arg {
250 int v;
251 int w;
252 int x;
253 int y;
254 int z;
255 } lunar_longitude_data[] = {
256 { 6288774, 0, 0, 1, 0 },
257 { 1274027, 2, 0, -1, 0 },
258 { 658314, 2, 0, 0, 0 },
259 { 213618, 0, 0, 2, 0 },
260 { -185116, 0, 1, 0, 0 },
261 { -114332, 0, 0, 0, 2 },
262 { 58793, 2, 0, -2, 0 },
263 { 57066, 2, -1, -1, 0 },
264 { 53322, 2, 0, 1, 0 },
265 { 45758, 2, -1, 0, 0 },
266 { -40923, 0, 1, -1, 0 },
267 { -34720, 1, 0, 0, 0 },
268 { -30383, 0, 1, 1, 0 },
269 { 15327, 2, 0, 0, -2 },
270 { -12528, 0, 0, 1, 2 },
271 { 10980, 0, 0, 1, -2 },
272 { 10675, 4, 0, -1, 0 },
273 { 10034, 0, 0, 3, 0 },
274 { 8548, 4, 0, -2, 0 },
275 { -7888, 2, 1, -1, 0 },
276 { -6766, 2, 1, 0, 0 },
277 { -5163, 1, 0, -1, 0 },
278 { 4987, 1, 1, 0, 0 },
279 { 4036, 2, -1, 1, 0 },
280 { 3994, 2, 0, 2, 0 },
281 { 3861, 4, 0, 0, 0 },
282 { 3665, 2, 0, -3, 0 },
283 { -2689, 0, 1, -2, 0 },
284 { -2602, 2, 0, -1, 2 },
285 { 2390, 2, -1, -2, 0 },
286 { -2348, 1, 0, 1, 0 },
287 { 2236, 2, -2, 0, 0 },
288 { -2120, 0, 1, 2, 0 },
289 { -2069, 0, 2, 0, 0 },
290 { 2048, 2, -2, -1, 0 },
291 { -1773, 2, 0, 1, -2 },
292 { -1595, 2, 0, 0, 2 },
293 { 1215, 4, -1, -1, 0 },
294 { -1110, 0, 0, 2, 2 },
295 { -892, 3, 0, -1, 0 },
296 { -810, 2, 1, 1, 0 },
297 { 759, 4, -1, -2, 0 },
298 { -713, 0, 2, -1, 0 },
299 { -700, 2, 2, -1, 0 },
300 { 691, 2, 1, -2, 0 },
301 { 596, 2, -1, 0, -2 },
302 { 549, 4, 0, 1, 0 },
303 { 537, 0, 0, 4, 0 },
304 { 520, 4, -1, 0, 0 },
305 { -487, 1, 0, -2, 0 },
306 { -399, 2, 1, 0, -2 },
307 { -381, 0, 0, 2, -2 },
308 { 351, 1, 1, 1, 0 },
309 { -340, 3, 0, -2, 0 },
310 { 330, 4, 0, -3, 0 },
311 { 327, 2, -1, 2, 0 },
312 { -323, 0, 2, 1, 0 },
313 { 299, 1, 1, -1, 0 },
314 { 294, 2, 0, 3, 0 },
315 };
316
317 /*
318 * Calculate the geocentric longitude of moon (in degrees) at moment $t.
319 * Ref: Sec.(14.6), Eq.(14.48)
320 */
321 double
lunar_longitude(double t)322 lunar_longitude(double t)
323 {
324 double c = julian_centuries(t);
325 double nu = nutation(t);
326
327 double L_prime = lunar_longitude_mean(c);
328 double D = lunar_elongation(c);
329 double M = solar_anomaly(c);
330 double M_prime = lunar_anomaly(c);
331 double F = moon_node(c);
332 double E = 1.0 - 0.002516 * c - 0.0000074 * c*c;
333
334 double sum = 0.0;
335 const struct lunar_longitude_arg *arg;
336 for (size_t i = 0; i < nitems(lunar_longitude_data); i++) {
337 arg = &lunar_longitude_data[i];
338 sum += arg->v * pow(E, abs(arg->x)) * sin_deg(
339 arg->w * D +
340 arg->x * M +
341 arg->y * M_prime +
342 arg->z * F);
343 }
344 double correction = sum / 1e6;
345
346 double venus = sin_deg(119.75 + 131.849 * c) * 3958 / 1e6;
347 double jupiter = sin_deg(53.09 + 479264.29 * c) * 318 / 1e6;
348 double flat_earth = sin_deg(L_prime - F) * 1962 / 1e6;
349
350 return mod_f((L_prime + correction + venus + jupiter +
351 flat_earth + nu), 360);
352 }
353
354 /*
355 * Argument data used by 'lunar_latitude()'.
356 * Ref: Sec.(14.6), Table(14.6)
357 */
358 static const struct lunar_latitude_arg {
359 int v;
360 int w;
361 int x;
362 int y;
363 int z;
364 } lunar_latitude_data[] = {
365 { 5128122, 0, 0, 0, 1 },
366 { 280602, 0, 0, 1, 1 },
367 { 277693, 0, 0, 1, -1 },
368 { 173237, 2, 0, 0, -1 },
369 { 55413, 2, 0, -1, 1 },
370 { 46271, 2, 0, -1, -1 },
371 { 32573, 2, 0, 0, 1 },
372 { 17198, 0, 0, 2, 1 },
373 { 9266, 2, 0, 1, -1 },
374 { 8822, 0, 0, 2, -1 },
375 { 8216, 2, -1, 0, -1 },
376 { 4324, 2, 0, -2, -1 },
377 { 4200, 2, 0, 1, 1 },
378 { -3359, 2, 1, 0, -1 },
379 { 2463, 2, -1, -1, 1 },
380 { 2211, 2, -1, 0, 1 },
381 { 2065, 2, -1, -1, -1 },
382 { -1870, 0, 1, -1, -1 },
383 { 1828, 4, 0, -1, -1 },
384 { -1794, 0, 1, 0, 1 },
385 { -1749, 0, 0, 0, 3 },
386 { -1565, 0, 1, -1, 1 },
387 { -1491, 1, 0, 0, 1 },
388 { -1475, 0, 1, 1, 1 },
389 { -1410, 0, 1, 1, -1 },
390 { -1344, 0, 1, 0, -1 },
391 { -1335, 1, 0, 0, -1 },
392 { 1107, 0, 0, 3, 1 },
393 { 1021, 4, 0, 0, -1 },
394 { 833, 4, 0, -1, 1 },
395 { 777, 0, 0, 1, -3 },
396 { 671, 4, 0, -2, 1 },
397 { 607, 2, 0, 0, -3 },
398 { 596, 2, 0, 2, -1 },
399 { 491, 2, -1, 1, -1 },
400 { -451, 2, 0, -2, 1 },
401 { 439, 0, 0, 3, -1 },
402 { 422, 2, 0, 2, 1 },
403 { 421, 2, 0, -3, -1 },
404 { -366, 2, 1, -1, 1 },
405 { -351, 2, 1, 0, 1 },
406 { 331, 4, 0, 0, 1 },
407 { 315, 2, -1, 1, 1 },
408 { 302, 2, -2, 0, -1 },
409 { -283, 0, 0, 1, 3 },
410 { -229, 2, 1, 1, -1 },
411 { 223, 1, 1, 0, -1 },
412 { 223, 1, 1, 0, 1 },
413 { -220, 0, 1, -2, -1 },
414 { -220, 2, 1, -1, -1 },
415 { -185, 1, 0, 1, 1 },
416 { 181, 2, -1, -2, -1 },
417 { -177, 0, 1, 2, 1 },
418 { 176, 4, 0, -2, -1 },
419 { 166, 4, -1, -1, -1 },
420 { -164, 1, 0, 1, -1 },
421 { 132, 4, 0, 1, -1 },
422 { -119, 1, 0, -1, -1 },
423 { 115, 4, -1, 0, -1 },
424 { 107, 2, -2, 0, 1 },
425 };
426
427 /*
428 * Calculate the geocentric latitude of moon (in degrees) at moment $t.
429 * Lunar latitude ranges from about -6 to 6 degress.
430 * Ref: Sec.(14.6), Eq.(14.63)
431 */
432 double
lunar_latitude(double t)433 lunar_latitude(double t)
434 {
435 double c = julian_centuries(t);
436
437 double L_prime = lunar_longitude_mean(c);
438 double D = lunar_elongation(c);
439 double M = solar_anomaly(c);
440 double M_prime = lunar_anomaly(c);
441 double F = moon_node(c);
442 double E = 1.0 - 0.002516 * c - 0.0000074 * c*c;
443
444 double sum = 0.0;
445 const struct lunar_latitude_arg *arg;
446 for (size_t i = 0; i < nitems(lunar_latitude_data); i++) {
447 arg = &lunar_latitude_data[i];
448 sum += arg->v * pow(E, abs(arg->x)) * sin_deg(
449 arg->w * D +
450 arg->x * M +
451 arg->y * M_prime +
452 arg->z * F);
453 }
454 double beta = sum / 1e6;
455
456 double venus = (sin_deg(119.75 + 131.849 * c + F) +
457 sin_deg(119.75 + 131.849 * c - F)) * 175 / 1e6;
458 double flat_earth = (-2235 * sin_deg(L_prime) +
459 127 * sin_deg(L_prime - M_prime) -
460 115 * sin_deg(L_prime + M_prime)) / 1e6;
461 double extra = sin_deg(313.45 + 481266.484 * c) * 382 / 1e6;
462
463 return (beta + venus + flat_earth + extra);
464 }
465
466 /*
467 * Argument data used by 'lunar_distance()'.
468 * Ref: Sec.(14.6), Table(14.7)
469 */
470 static const struct lunar_distance_arg {
471 int v;
472 int w;
473 int x;
474 int y;
475 int z;
476 } lunar_distance_data[] = {
477 { -20905355, 0, 0, 1, 0 },
478 { -3699111, 2, 0, -1, 0 },
479 { -2955968, 2, 0, 0, 0 },
480 { -569925, 0, 0, 2, 0 },
481 { 48888, 0, 1, 0, 0 },
482 { -3149, 0, 0, 0, 2 },
483 { 246158, 2, 0, -2, 0 },
484 { -152138, 2, -1, -1, 0 },
485 { -170733, 2, 0, 1, 0 },
486 { -204586, 2, -1, 0, 0 },
487 { -129620, 0, 1, -1, 0 },
488 { 108743, 1, 0, 0, 0 },
489 { 104755, 0, 1, 1, 0 },
490 { 10321, 2, 0, 0, -2 },
491 { 0, 0, 0, 1, 2 },
492 { 79661, 0, 0, 1, -2 },
493 { -34782, 4, 0, -1, 0 },
494 { -23210, 0, 0, 3, 0 },
495 { -21636, 4, 0, -2, 0 },
496 { 24208, 2, 1, -1, 0 },
497 { 30824, 2, 1, 0, 0 },
498 { -8379, 1, 0, -1, 0 },
499 { -16675, 1, 1, 0, 0 },
500 { -12831, 2, -1, 1, 0 },
501 { -10445, 2, 0, 2, 0 },
502 { -11650, 4, 0, 0, 0 },
503 { 14403, 2, 0, -3, 0 },
504 { -7003, 0, 1, -2, 0 },
505 { 0, 2, 0, -1, 2 },
506 { 10056, 2, -1, -2, 0 },
507 { 6322, 1, 0, 1, 0 },
508 { -9884, 2, -2, 0, 0 },
509 { 5751, 0, 1, 2, 0 },
510 { 0, 0, 2, 0, 0 },
511 { -4950, 2, -2, -1, 0 },
512 { 4130, 2, 0, 1, -2 },
513 { 0, 2, 0, 0, 2 },
514 { -3958, 4, -1, -1, 0 },
515 { 0, 0, 0, 2, 2 },
516 { 3258, 3, 0, -1, 0 },
517 { 2616, 2, 1, 1, 0 },
518 { -1897, 4, -1, -2, 0 },
519 { -2117, 0, 2, -1, 0 },
520 { 2354, 2, 2, -1, 0 },
521 { 0, 2, 1, -2, 0 },
522 { 0, 2, -1, 0, -2 },
523 { -1423, 4, 0, 1, 0 },
524 { -1117, 0, 0, 4, 0 },
525 { -1571, 4, -1, 0, 0 },
526 { -1739, 1, 0, -2, 0 },
527 { 0, 2, 1, 0, -2 },
528 { -4421, 0, 0, 2, -2 },
529 { 0, 1, 1, 1, 0 },
530 { 0, 3, 0, -2, 0 },
531 { 0, 4, 0, -3, 0 },
532 { 0, 2, -1, 2, 0 },
533 { 1165, 0, 2, 1, 0 },
534 { 0, 1, 1, -1, 0 },
535 { 0, 2, 0, 3, 0 },
536 { 8752, 2, 0, -1, -2 },
537 };
538
539 /*
540 * Calculate the distance to moon (in meters) at moment $t.
541 * Ref: Sec.(14.6), Eq.(14.65)
542 */
543 double
lunar_distance(double t)544 lunar_distance(double t)
545 {
546 double c = julian_centuries(t);
547
548 double D = lunar_elongation(c);
549 double M = solar_anomaly(c);
550 double M_prime = lunar_anomaly(c);
551 double F = moon_node(c);
552 double E = 1.0 - 0.002516 * c - 0.0000074 * c*c;
553
554 double correction = 0.0;
555 const struct lunar_distance_arg *arg;
556 for (size_t i = 0; i < nitems(lunar_distance_data); i++) {
557 arg = &lunar_distance_data[i];
558 correction += arg->v * pow(E, abs(arg->x)) * cos_deg(
559 arg->w * D +
560 arg->x * M +
561 arg->y * M_prime +
562 arg->z * F);
563 }
564
565 return 385000560.0 + correction;
566 }
567
568 /*
569 * Calculate the altitude of moon (in degrees) above the horizon at
570 * location ($latitude, $longitude) and moment $t, ignoring parallax
571 * and refraction.
572 * NOTE: This calculates the geocentric altitude viewed from the center
573 * of the Earth.
574 * Ref: Sec.(14.6), Eq.(14.64)
575 */
576 double
lunar_altitude(double t,double latitude,double longitude)577 lunar_altitude(double t, double latitude, double longitude)
578 {
579 double lambda = lunar_longitude(t);
580 double beta = lunar_latitude(t);
581 double alpha = right_ascension(t, beta, lambda);
582 double delta = declination(t, beta, lambda);
583 double theta = sidereal_from_moment(t);
584 double H = mod_f(theta + longitude - alpha, 360);
585
586 double v = (sin_deg(latitude) * sin_deg(delta) +
587 cos_deg(latitude) * cos_deg(delta) * cos_deg(H));
588 return mod3_f(arcsin_deg(v), -180, 180);
589 }
590
591 /*
592 * Parallax of moon at moment $t and location ($latitude, $longitude).
593 * Ref: Sec.(14.6), Eq.(14.66)
594 */
595 static double
lunar_parallax(double t,double latitude,double longitude)596 lunar_parallax(double t, double latitude, double longitude)
597 {
598 double geo = lunar_altitude(t, latitude, longitude);
599 double distance = lunar_distance(t);
600 /* Equatorial horizontal parallax of the moon */
601 double sin_pi = 6378140.0 / distance;
602 return arcsin_deg(sin_pi * cos_deg(geo));
603 }
604
605 /*
606 * Calculate the topocentric altitude of moon viewed from the Earth surface
607 * at location ($latitude, $longitude) and at moment $t.
608 * Ref: Sec.(14.6), Eq.(14.67)
609 */
610 static double
lunar_altitude_topocentric(double t,double latitude,double longitude)611 lunar_altitude_topocentric(double t, double latitude, double longitude)
612 {
613 return (lunar_altitude(t, latitude, longitude) -
614 lunar_parallax(t, latitude, longitude));
615 }
616
617 /*
618 * Calculate the observed altitude of the upper limb of moon at location
619 * ($latitude, $longitude) and moment $t.
620 * Ref: Sec.(14.7), Eq.(14.82)
621 */
622 double
lunar_altitude_observed(double t,const struct location * loc)623 lunar_altitude_observed(double t, const struct location *loc)
624 {
625 double moon_radius = 16.0 / 60.0; /* 16 arcminutes */
626 return (lunar_altitude_topocentric(t, loc->latitude, loc->longitude) +
627 refraction(loc->elevation) + moon_radius);
628 }
629
630 /*
631 * Calculate the phase of the moon, defined as the difference in
632 * longitudes of the Sun and moon, at the given moment $t.
633 * Ref: Sec.(14.6), Eq.(14.56)
634 */
635 double
lunar_phase(double t)636 lunar_phase(double t)
637 {
638 double phi = mod_f(lunar_longitude(t) - solar_longitude(t), 360);
639
640 /*
641 * To check whether the above result conflicts with the time of
642 * new moon as calculated by the more precise 'nth_new_moon()'
643 */
644 double t0 = nth_new_moon(0);
645 int n = (int)lround((t - t0) / mean_synodic_month);
646 double tn = nth_new_moon(n);
647 double phi2 = mod_f((t - tn) / mean_synodic_month, 1) * 360;
648
649 /* prefer the approximation based on the 'nth_new_moon()' moment */
650 return (fabs(phi - phi2) > 180) ? phi2 : phi;
651 }
652
653 /*
654 * Calculate the moment of the next time at or after the given moment $t
655 * when the phase of the moon is $phi degree.
656 * Ref: Sec.(14.6), Eq.(14.58)
657 */
658 double
lunar_phase_atafter(double phi,double t)659 lunar_phase_atafter(double phi, double t)
660 {
661 double rate = mean_synodic_month / 360.0;
662 double phase = lunar_phase(t);
663 double tau = t + rate * mod_f(phi - phase, 360);
664
665 /* estimate range (within 2 days) */
666 double a = (t > tau - 2) ? t : tau - 2;
667 double b = tau + 2;
668
669 return invert_angular(lunar_phase, phi, a, b);
670 }
671
672 /*
673 * Calculate the moment of the new moon before the given moment $t.
674 * Ref: Sec.(14.6), Eq.(14.46)
675 */
676 double
new_moon_before(double t)677 new_moon_before(double t)
678 {
679 double t0 = nth_new_moon(0);
680 double phi = lunar_phase(t);
681 int n = (int)lround((t - t0) / mean_synodic_month - phi / 360.0);
682
683 int k = n - 1;
684 double t1 = nth_new_moon(k);
685 while (t1 < t) {
686 k++;
687 t1 = nth_new_moon(k);
688 }
689
690 return nth_new_moon(k-1);
691 }
692
693 /*
694 * Calculate the moment of the new moon at or after the given moment $t.
695 * Ref: Sec.(14.6), Eq.(14.47)
696 */
697 double
new_moon_atafter(double t)698 new_moon_atafter(double t)
699 {
700 double t0 = nth_new_moon(0);
701 double phi = lunar_phase(t);
702 int n = (int)lround((t - t0) / mean_synodic_month - phi / 360.0);
703
704 double t1 = nth_new_moon(n);
705 while (t1 < t) {
706 n++;
707 t1 = nth_new_moon(n);
708 }
709
710 return t1;
711 }
712
713 /*
714 * Calculate the moment of moonrise in standard time on fixed date $rd
715 * at location $loc.
716 * NOTE: Return an NaN if no moonrise.
717 * Ref: Sec.(14.7), Eq.(14.83)
718 */
719 double
moonrise(int rd,const struct location * loc)720 moonrise(int rd, const struct location *loc)
721 {
722 double t = (double)rd - loc->zone; /* universal time */
723 bool waning = lunar_phase(t) > 180.0;
724 /* lunar altitude at midnight */
725 double alt = lunar_altitude_observed(t, loc);
726 double offset = alt / (4.0 * (90.0 - fabs(loc->latitude)));
727
728 /* approximate rising time */
729 double t_approx = t;
730 if (waning) {
731 t_approx -= offset;
732 if (offset > 0)
733 t_approx += 1;
734 } else {
735 t_approx += 0.5 + offset;
736 }
737
738 /* binary search to determine the rising time */
739 const double eps = 30.0 / 3600 / 24; /* accuracy of 30 seconds */
740 double a = t_approx - 6.0/24.0;
741 double b = t_approx + 6.0/24.0;
742 double t_rise;
743 do {
744 t_rise = (a + b) / 2.0;
745 if (lunar_altitude_observed(t_rise, loc) > 0)
746 b = t_rise;
747 else
748 a = t_rise;
749 } while (fabs(a - b) >= eps);
750
751 if (t_rise < t + 1) {
752 t_rise += loc->zone; /* standard time */
753 /* may be just before to midnight */
754 return (t_rise > (double)rd) ? t_rise : (double)rd;
755 } else {
756 /* no moonrise on this day */
757 return NAN;
758 }
759 }
760
761 /*
762 * Calculate the moment of moonset in standard time on fixed date $rd
763 * at location $loc.
764 * NOTE: Return an NaN if no moonset.
765 * Ref: Sec.(14.7), Eq.(14.84)
766 */
767 double
moonset(int rd,const struct location * loc)768 moonset(int rd, const struct location *loc)
769 {
770 double t = (double)rd - loc->zone; /* universal time */
771 bool waxing = lunar_phase(t) < 180.0;
772 /* lunar altitude at midnight */
773 double alt = lunar_altitude_observed(t, loc);
774 double offset = alt / (4.0 * (90.0 - fabs(loc->latitude)));
775
776 /* approximate setting time */
777 double t_approx = t;
778 if (waxing) {
779 t_approx += offset;
780 if (offset <= 0)
781 t_approx += 1;
782 } else {
783 t_approx += 0.5 - offset;
784 }
785
786 /* binary search to determine the setting time */
787 const double eps = 30.0 / 3600 / 24; /* accuracy of 30 seconds */
788 double a = t_approx - 6.0/24.0;
789 double b = t_approx + 6.0/24.0;
790 double t_set;
791 do {
792 t_set = (a + b) / 2.0;
793 if (lunar_altitude_observed(t_set, loc) < 0)
794 b = t_set;
795 else
796 a = t_set;
797 } while (fabs(a - b) >= eps);
798
799 if (t_set < t + 1) {
800 t_set += loc->zone; /* standard time */
801 /* may be just before to midnight */
802 return (t_set > (double)rd) ? t_set : (double)rd;
803 } else {
804 /* no moonrise on this day */
805 return NAN;
806 }
807 }
808
809 /**************************************************************************/
810
811 /*
812 * Print moon information at the given moment $t (in standard time)
813 * and events in the year.
814 */
815 void
show_moon_info(double t,const struct location * loc)816 show_moon_info(double t, const struct location *loc)
817 {
818 char buf[64];
819 int rd = (int)floor(t);
820 double t_u = t - loc->zone; /* universal time */
821
822 /*
823 * Lunar phase
824 */
825
826 double eps = 1e-5;
827 double phi = lunar_phase(t_u);
828 bool waxing = (phi <= 180);
829 bool crescent = (phi <= 90 || phi > 270);
830 const char *phase_name;
831 if (fabs(phi) < eps || fabs(phi - 360) < eps)
832 phase_name = "New Moon";
833 else if (fabs(phi - 90) < eps)
834 phase_name = "First Quarter";
835 else if (fabs(phi - 180) < eps)
836 phase_name = "Full Moon";
837 else if (fabs(phi - 270) < eps)
838 phase_name = "Last Quarter";
839 else
840 phase_name = NULL;
841
842 if (phase_name) {
843 snprintf(buf, sizeof(buf), "%s", phase_name);
844 } else {
845 snprintf(buf, sizeof(buf), "%s %s",
846 (waxing ? "Waxing" : "Waning"),
847 (crescent ? "Crescent" : "Gibbous"));
848 }
849 printf("Moon phase: %.2lf° (%s)\n", phi, buf);
850
851 /*
852 * Moon position
853 */
854 double lon = lunar_longitude(t_u);
855 double alt = lunar_altitude_observed(t_u, loc);
856 printf("Moon position: %.4lf° (longitude), %.4lf° (altitude)\n",
857 lon, alt);
858
859 /*
860 * Moon rise and set
861 */
862
863 double moments[2] = { moonrise(rd, loc), moonset(rd, loc) };
864 const char *names[2] = { "Moonrise", "Moonset" };
865 if (!isnan(moments[0]) && !isnan(moments[1]) &&
866 moments[0] > moments[1]) {
867 double t_tmp = moments[0];
868 moments[0] = moments[1];
869 moments[1] = t_tmp;
870 const char *p = names[0];
871 names[0] = names[1];
872 names[1] = p;
873 }
874 for (size_t i = 0; i < nitems(moments); i++) {
875 if (isnan(moments[i]))
876 snprintf(buf, sizeof(buf), "(null)");
877 else
878 format_time(buf, sizeof(buf), moments[i]);
879 printf("%-8s: %s\n", names[i], buf);
880 }
881
882 /*
883 * Moon phases in the year
884 */
885
886 int year = gregorian_year_from_fixed(rd);
887 struct date date = { year, 1, 1 };
888 double t_begin = fixed_from_gregorian(&date) - loc->zone;
889 date.year++;
890 double t_end = fixed_from_gregorian(&date) - loc->zone;
891
892 printf("\nLunar events in year %d:\n", year);
893 printf("%19s %19s %19s %19s\n",
894 "New Moon", "First Quarter", "Full Moon", "Last Quarter");
895
896 double t_newmoon = t_begin;
897 while ((t_newmoon = new_moon_atafter(t_newmoon)) < t_end) {
898 t_newmoon += loc->zone; /* to standard time */
899 gregorian_from_fixed((int)floor(t_newmoon), &date);
900 format_time(buf, sizeof(buf), t_newmoon);
901 printf("%d-%02d-%02d %s",
902 date.year, date.month, date.day, buf);
903
904 /*
905 * first quarter, full moon, last quarter
906 */
907 double t_event = t_newmoon;
908 int phi_events[] = { 90, 180, 270 };
909 for (size_t i = 0; i < nitems(phi_events); i++) {
910 t_event = lunar_phase_atafter(phi_events[i], t_event);
911 t_event += loc->zone;
912 gregorian_from_fixed((int)floor(t_event), &date);
913 format_time(buf, sizeof(buf), t_event);
914 printf(" %d-%02d-%02d %s",
915 date.year, date.month, date.day, buf);
916 }
917 printf("\n");
918
919 /* go to the next new moon */
920 t_newmoon += 28;
921 }
922 }
923