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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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