1#ifndef __DT_ASTRO_LUNAR_C__ 2#define __DT_ASTRO_LUNAR_C__ 3#include "dt_astro.h" 4#define NV_1E6 1000000.0 5 6STATIC_INLINE 7int 8lunar_longitude( mpfr_t *result, mpfr_t *moment ) { 9 10 mpfr_t C, mean_moon, elongation, solar_anomaly, lunar_anomaly, moon_node, E, correction, venus, jupiter, flat_earth, N, fullangle; 11 12 mpfr_init(C); 13 julian_centuries( &C, moment ); 14 15 { 16 mpfr_t a, b, c, d, e; 17 18 mpfr_init(mean_moon); 19 mpfr_init_set_d(a, 218.316591, GMP_RNDN); 20 mpfr_init_set_d(b, 481267.88134236, GMP_RNDN); 21 mpfr_init_set_d(c, -0.0013268, GMP_RNDN); 22 mpfr_init_set_ui(d, 1, GMP_RNDN); 23 mpfr_div_ui(d, d, 538841, GMP_RNDN); 24 mpfr_init_set_si(e, -1, GMP_RNDN); 25 mpfr_div_ui(e, e, 65194000, GMP_RNDN); 26 27 polynomial( &mean_moon, &C, 5, &a, &b, &c, &d, &e ); 28 mpfr_clear(a); 29 mpfr_clear(b); 30 mpfr_clear(c); 31 mpfr_clear(d); 32 mpfr_clear(e); 33 } 34 35 { 36 mpfr_t a, b, c, d, e; 37 mpfr_init(elongation); 38 39 mpfr_init_set_d(a, 297.8502042, GMP_RNDN); 40 mpfr_init_set_d(b, 445267.1115168, GMP_RNDN); 41 mpfr_init_set_d(c, -0.00163, GMP_RNDN); 42 mpfr_init_set_ui(d, 1, GMP_RNDN); 43 mpfr_div_ui(d, d, 545868, GMP_RNDN); 44 mpfr_init_set_si(e, -1, GMP_RNDN); 45 mpfr_div_ui(e, e, 113065000, GMP_RNDN); 46 polynomial( &elongation, &C, 5, &a, &b, &c, &d, &e ); 47 mpfr_clear(a); 48 mpfr_clear(b); 49 mpfr_clear(c); 50 mpfr_clear(d); 51 mpfr_clear(e); 52 } 53 54 { 55 mpfr_t a, b, c, d; 56 mpfr_init(solar_anomaly); 57 mpfr_init_set_d(a, 357.5291092, GMP_RNDN); 58 mpfr_init_set_d(b, 35999.0502909, GMP_RNDN); 59 mpfr_init_set_d(c, -0.0001536, GMP_RNDN); 60 mpfr_init_set_ui(d, 1, GMP_RNDN); 61 mpfr_div_ui(d, d, 24490000, GMP_RNDN); 62 polynomial( &solar_anomaly, &C, 4, &a, &b, &c, &d ); 63 mpfr_clear(a); 64 mpfr_clear(b); 65 mpfr_clear(c); 66 mpfr_clear(d); 67 } 68 69 { 70 mpfr_t a, b, c, d, e; 71 mpfr_init(lunar_anomaly); 72 73 mpfr_init_set_d(a, 134.9634114, GMP_RNDN); 74 mpfr_init_set_d(b, 477198.8676313, GMP_RNDN); 75 mpfr_init_set_d(c, 0.0008997, GMP_RNDN); 76 mpfr_init_set_ui(d, 1, GMP_RNDN); 77 mpfr_div_ui(d, d, 69699, GMP_RNDN); 78 mpfr_init_set_si(e, -1, GMP_RNDN); 79 mpfr_div_ui(e, e, 14712000, GMP_RNDN); 80 polynomial( &lunar_anomaly, &C, 5, &a, &b, &c, &d, &e); 81 mpfr_clear(a); 82 mpfr_clear(b); 83 mpfr_clear(c); 84 mpfr_clear(d); 85 mpfr_clear(e); 86 } 87 88 { 89 mpfr_t a, b, c, d, e; 90 mpfr_init(moon_node); 91 mpfr_init_set_d(a, 93.2720993, GMP_RNDN); 92 mpfr_init_set_d(b, 483202.0175273, GMP_RNDN); 93 mpfr_init_set_d(c, -0.0034029, GMP_RNDN); 94 mpfr_init_set_si(d, -1, GMP_RNDN); 95 mpfr_div_ui(d, d, 3526000, GMP_RNDN); 96 mpfr_init_set_ui(e, 1, GMP_RNDN); 97 mpfr_div_ui(e, e, 863310000, GMP_RNDN); 98 polynomial(&moon_node, &C, 5, &a, &b, &c, &d, &e); 99 mpfr_clear(a); 100 mpfr_clear(b); 101 mpfr_clear(c); 102 mpfr_clear(d); 103 mpfr_clear(e); 104 } 105 106 { 107 mpfr_t a, b, c; 108 mpfr_init(E); 109 mpfr_init_set_ui(a, 1, GMP_RNDN); 110 mpfr_init_set_d(b, -0.002516, GMP_RNDN); 111 mpfr_init_set_d(c, -0.0000074, GMP_RNDN); 112 polynomial( &E, &C, 3, &a, &b, &c ); 113 mpfr_clear(a); 114 mpfr_clear(b); 115 mpfr_clear(c); 116 } 117 118 { 119 int i; 120 mpfr_t fugly; 121 mpfr_init_set_ui(fugly, 0, GMP_RNDN); 122 123 for(i = 0; i < LUNAR_LONGITUDE_ARGS_SIZE; i++) { 124 mpfr_t a, b, v, w, x, y, z; 125 mpfr_init_set_d( v, LUNAR_LONGITUDE_ARGS[i][0], GMP_RNDN ); 126 mpfr_init_set_d( w, LUNAR_LONGITUDE_ARGS[i][1], GMP_RNDN ); 127 mpfr_init_set_d( x, LUNAR_LONGITUDE_ARGS[i][2], GMP_RNDN ); 128 mpfr_init_set_d( y, LUNAR_LONGITUDE_ARGS[i][3], GMP_RNDN ); 129 mpfr_init_set_d( z, LUNAR_LONGITUDE_ARGS[i][4], GMP_RNDN ); 130 131 mpfr_init(b); 132 mpfr_pow(b, E, x, GMP_RNDN); 133 134 mpfr_mul(w, w, elongation, GMP_RNDN); 135 mpfr_mul(x, x, solar_anomaly, GMP_RNDN); 136 mpfr_mul(y, y, lunar_anomaly, GMP_RNDN); 137 mpfr_mul(z, z, moon_node, GMP_RNDN); 138 139 mpfr_init_set(a, w, GMP_RNDN); 140 mpfr_add(a, a, x, GMP_RNDN); 141 mpfr_add(a, a, y, GMP_RNDN); 142 mpfr_add(a, a, z, GMP_RNDN); 143 dt_astro_sin(&a, &a); 144 145 mpfr_mul(a, a, v, GMP_RNDN); 146 mpfr_mul(a, a, b, GMP_RNDN); 147 mpfr_add(fugly, fugly, a, GMP_RNDN); 148 149 mpfr_clear(a); 150 mpfr_clear(b); 151 mpfr_clear(v); 152 mpfr_clear(w); 153 mpfr_clear(x); 154 mpfr_clear(y); 155 mpfr_clear(z); 156 } 157 158 mpfr_init_set_d( correction, 0.000001, GMP_RNDN ); 159 mpfr_mul( correction, correction, fugly, GMP_RNDN); 160 mpfr_clear(fugly); 161 } 162 163 { 164 mpfr_t a, b; 165 mpfr_init(venus); 166 mpfr_init_set_d(a, 119.75, GMP_RNDN); 167 mpfr_init_set(b, C, GMP_RNDN); 168 mpfr_mul_d(b, b, 131.849, GMP_RNDN); 169 170 mpfr_add(a, a, b, GMP_RNDN); 171 dt_astro_sin(&a, &a); 172 mpfr_mul_d(venus, a, 0.003957, GMP_RNDN ); 173 mpfr_clear(a); 174 mpfr_clear(b); 175 } 176 177 { 178 mpfr_t a, b; 179 mpfr_init(jupiter); 180 mpfr_init_set_d(a, 53.09, GMP_RNDN); 181 mpfr_init_set(b, C, GMP_RNDN); 182 mpfr_mul_d(b, b, 479264.29, GMP_RNDN); 183 184 mpfr_add(a, a, b, GMP_RNDN); 185 dt_astro_sin(&a, &a); 186 mpfr_mul_d(jupiter, a, 0.000318, GMP_RNDN ); 187 mpfr_clear(a); 188 mpfr_clear(b); 189 } 190 191 { 192 mpfr_t a; 193 mpfr_init(flat_earth); 194 mpfr_init_set(a, mean_moon, GMP_RNDN); 195 mpfr_sub(a, a, moon_node, GMP_RNDN); 196 dt_astro_sin(&a, &a); 197 mpfr_mul_d(flat_earth, a, 0.001962, GMP_RNDN); 198 mpfr_clear(a); 199 } 200 201 mpfr_set(*result, mean_moon, GMP_RNDN); 202 mpfr_add(*result, *result, correction, GMP_RNDN); 203 mpfr_add(*result, *result, venus, GMP_RNDN); 204 mpfr_add(*result, *result, jupiter, GMP_RNDN); 205 mpfr_add(*result, *result, flat_earth, GMP_RNDN); 206 207#ifdef ANNOYING_DEBUG 208#if (ANNOYING_DEBUG) 209mpfr_fprintf(stderr, 210 "mean_moon = %.10RNf\ncorrection = %.10RNf\nvenus = %.10RNf\njupiter = %.10RNf\nflat_earth = %.10RNf\n", 211 mean_moon, 212 correction, 213 venus, 214 jupiter, 215 flat_earth); 216#endif 217#endif 218 219 mpfr_init(N); 220 nutation(&N, moment); 221 mpfr_add(*result, *result, N, GMP_RNDN); 222 223 mpfr_init_set_ui(fullangle, 360, GMP_RNDN); 224 225#ifdef ANNOYING_DEBUG 226#if (ANNOYING_DEBUG) 227mpfr_fprintf(stderr, "lunar = mod(%.10RNf) = ", *result ); 228#endif 229#endif 230 dt_astro_mod(result, result, &fullangle); 231#ifdef ANNOYING_DEBUG 232#if (ANNOYING_DEBUG) 233mpfr_fprintf(stderr, "%.10RNf\n", *result ); 234#endif 235#endif 236 237 mpfr_clear(C); 238 mpfr_clear(mean_moon); 239 mpfr_clear(elongation); 240 mpfr_clear(solar_anomaly); 241 mpfr_clear(lunar_anomaly); 242 mpfr_clear(moon_node); 243 mpfr_clear(E); 244 mpfr_clear(correction); 245 mpfr_clear(venus); 246 mpfr_clear(jupiter); 247 mpfr_clear(flat_earth); 248 mpfr_clear(N); 249 mpfr_clear(fullangle); 250 return 1; 251} 252 253STATIC_INLINE 254int 255lunar_phase( mpfr_t *result, mpfr_t *moment ) { 256 mpfr_t sl, ll, fullangle; 257 mpfr_init(sl); 258 mpfr_init(ll); 259 mpfr_init_set_ui(fullangle, 360, GMP_RNDN); 260 261 solar_longitude( &sl, moment ); 262 lunar_longitude( &ll, moment ); 263 mpfr_sub(*result, ll, sl, GMP_RNDN ); 264 dt_astro_mod(result, result, &fullangle); 265 266 mpfr_clear(sl); 267 mpfr_clear(ll); 268 mpfr_clear(fullangle); 269 return 1; 270} 271 272STATIC_INLINE 273void 274adjust_lunar_phase_to_zero(mpfr_t *result) { 275 mpfr_t ll, delta; 276 int mode = -1; 277 int loop = 1; 278 int count = 0; 279 /* Adjust values so that it's as close as possible to 0 degrees. 280 * if we have a delta of 1 degree, then we're about 281 * 1 / ( 360 / MEAN_SYNODIC_MONTH ) 282 * days apart 283 */ 284 285 mpfr_init(ll); 286 mpfr_init_set_d(delta, 0.0001, GMP_RNDN); 287 288 while (loop) { 289 int flipped = mode; 290 mpfr_t new_moment; 291 count++; 292 mpfr_init(new_moment); 293 lunar_phase(&ll, result); 294#if (TRACE) 295mpfr_fprintf(stderr, 296 "Adjusting ll from (%.30RNf) moment is %.5RNf delta is %.30RNf\n", ll, *result, delta); 297#endif 298 /* longitude was greater than 180, so we're looking to add a few 299 * degrees to make it close to 360 ( 0 ) 300 */ 301 if (mpfr_cmp_ui( ll, 180 ) > 0) { 302 mode = 1; 303 mpfr_sub_ui(delta, ll, 360, GMP_RNDN); 304 mpfr_div_d(delta, delta, 360 / MEAN_SYNODIC_MONTH, GMP_RNDN); 305 mpfr_add( new_moment, *result, delta, GMP_RNDN ); 306#if (TRACE) 307mpfr_fprintf(stderr, "add %.30RNf -> %.30RNf\n", *result, new_moment); 308#endif 309 mpfr_set(*result, new_moment, GMP_RNDN); 310 if (mpfr_cmp(new_moment, *result) == 0) { 311 loop = 0; 312 } 313 } else if (mpfr_cmp_ui( ll, 180 ) < 0 ) { 314 if ( mpfr_cmp_d( ll, 0.000000000000000000001 ) < 0) { 315 loop = 0; 316 } else { 317 mode = 0; 318 mpfr_sub_ui(delta, ll, 0, GMP_RNDN); 319 mpfr_div_d(delta, delta, 360 / MEAN_SYNODIC_MONTH, GMP_RNDN); 320 mpfr_sub( new_moment, *result, delta, GMP_RNDN ); 321#if (TRACE) 322mpfr_fprintf(stderr, "sub %.120RNf -> %.120RNf\n", *result, new_moment); 323#endif 324 if (mpfr_cmp(new_moment, *result) == 0) { 325 loop = 0; 326 } 327 mpfr_set(*result, new_moment, GMP_RNDN); 328 } 329 } else { 330 loop = 0; 331 } 332 if (flipped != -1 && flipped != mode) { 333 mpfr_div_d(delta, delta, 1.1, GMP_RNDN); 334 } 335 mpfr_clear(new_moment); 336 } 337 mpfr_clear(delta); 338 mpfr_clear(ll); 339} 340 341STATIC_INLINE 342int 343nth_new_moon( mpfr_t *result, int n_int ) { 344 mpfr_t n, k, C, approx, E, solar_anomaly, lunar_anomaly, moon_argument, omega, extra, correction, additional; 345 346#if(0) 347PerlIO_printf(PerlIO_stderr(), "nth_new_moon = %d\n", n_int ); 348#endif 349 if ( dt_astro_global_cache.cache_size > n_int ) { 350 mpfr_t *cached = dt_astro_global_cache.cache[n_int]; 351 if (cached != NULL) { 352#if(0) 353 PerlIO_printf(PerlIO_stderr(), "Cache HIT for %d\n", n_int); 354#endif 355 mpfr_set( *result, *cached, GMP_RNDN ); 356 return 1; 357 } 358 } 359 360 mpfr_init_set_ui( n, n_int, GMP_RNDN ); 361 362 /* k = n - 24724 */ 363 mpfr_init_set(k, n, GMP_RNDN); 364 mpfr_sub_ui(k, k, 24724, GMP_RNDN ); 365 366 /* c = k / 1236.85 */ 367 mpfr_init_set(C, k, GMP_RNDN ); 368 mpfr_div_d(C, C, 1236.85, GMP_RNDN); 369 370 { 371 mpfr_t a, b, c, d, e; 372 mpfr_init(approx); 373 mpfr_init_set_d(a, 730125.59765, GMP_RNDN ); 374 mpfr_init_set_d(b, MEAN_SYNODIC_MONTH * 1236.85, GMP_RNDN ); 375 mpfr_init_set_d(c, 0.0001337, GMP_RNDN ); 376 mpfr_init_set_d(d, -0.000000150, GMP_RNDN ); 377 mpfr_init_set_d(e, 0.00000000073, GMP_RNDN ); 378 polynomial( &approx, &C, 5, &a, &b, &c, &d, &e ); 379 mpfr_clear(a); 380 mpfr_clear(b); 381 mpfr_clear(c); 382 mpfr_clear(d); 383 mpfr_clear(e); 384#ifdef ANNOYING_DEBUG 385#if (ANNOYING_DEBUG) 386mpfr_fprintf(stderr, 387 "approx = %.10RNf\n", approx); 388#endif 389#endif 390 } 391 392 { 393 mpfr_t a, b, c; 394 mpfr_init(E); 395 mpfr_init_set_ui(a, 1, GMP_RNDN); 396 mpfr_init_set_d(b, -0.002516, GMP_RNDN ); 397 mpfr_init_set_d(c, -0.0000074, GMP_RNDN ); 398 polynomial( &E, &C, 3, &a, &b, &c ); 399 mpfr_clear(a); 400 mpfr_clear(b); 401 mpfr_clear(c); 402 } 403 404 { 405 mpfr_t a, b, c, d; 406 mpfr_init(solar_anomaly); 407 mpfr_init_set_d(a, 2.5534, GMP_RNDN); 408 mpfr_init_set_d(b, 1236.85, GMP_RNDN); 409 mpfr_mul_d(b, b, 29.10535669, GMP_RNDN); 410 mpfr_init_set_d(c, -0.0000218, GMP_RNDN ); 411 mpfr_init_set_d(d, -0.00000011, GMP_RNDN ); 412 polynomial( &solar_anomaly, &C, 4, &a, &b, &c, &d); 413 mpfr_clear(a); 414 mpfr_clear(b); 415 mpfr_clear(c); 416 mpfr_clear(d); 417 } 418 419 { 420 mpfr_t a, b, c, d, e; 421 mpfr_init(lunar_anomaly); 422 mpfr_init_set_d(a, 201.5643, GMP_RNDN); 423 mpfr_init_set_d(b, 385.81693528 * 1236.85, GMP_RNDN); 424 mpfr_init_set_d(c, 0.0107438, GMP_RNDN); 425 mpfr_init_set_d(d, 0.00001239, GMP_RNDN); 426 mpfr_init_set_d(e, -0.000000058, GMP_RNDN); 427 polynomial( &lunar_anomaly, &C, 5, &a, &b, &c, &d, &e); 428 mpfr_clear(a); 429 mpfr_clear(b); 430 mpfr_clear(c); 431 mpfr_clear(d); 432 mpfr_clear(e); 433 } 434 435 { 436 mpfr_t a, b, c, d, e; 437 mpfr_init(moon_argument); 438 mpfr_init_set_d(a, 160.7108, GMP_RNDN); 439 mpfr_init_set_d(b, 390.67050274 * 1236.85, GMP_RNDN); 440 mpfr_init_set_d(c, -0.0016431, GMP_RNDN); 441 mpfr_init_set_d(d, -0.00000227, GMP_RNDN); 442 mpfr_init_set_d(e, 0.000000011, GMP_RNDN); 443 polynomial( &moon_argument, &C, 5, &a, &b, &c, &d, &e); 444 mpfr_clear(a); 445 mpfr_clear(b); 446 mpfr_clear(c); 447 mpfr_clear(d); 448 mpfr_clear(e); 449 } 450 451 { 452 mpfr_t a, b, c, d; 453 mpfr_init(omega); 454 mpfr_init_set_d(a, 124.7746, GMP_RNDN); 455 mpfr_init_set_d(b, -1.56375580 * 1236.85, GMP_RNDN); 456 mpfr_init_set_d(c, 0.0020691, GMP_RNDN); 457 mpfr_init_set_d(d, 0.00000215, GMP_RNDN); 458 polynomial( &omega, &C, 4, &a, &b, &c, &d); 459 mpfr_clear(a); 460 mpfr_clear(b); 461 mpfr_clear(c); 462 mpfr_clear(d); 463 } 464 465 { 466 mpfr_t a, b, c; 467 mpfr_init(extra); 468 mpfr_init_set_d(a, 299.77, GMP_RNDN); 469 mpfr_init_set_d(b, 132.8475848, GMP_RNDN); 470 mpfr_init_set_d(c, -0.009173, GMP_RNDN); 471 polynomial(&extra, &C, 3, &a, &b, &c); 472 dt_astro_sin(&extra, &extra); 473 mpfr_mul_d(extra, extra, 0.000325, GMP_RNDN); 474 mpfr_clear(a); 475 mpfr_clear(b); 476 mpfr_clear(c); 477 } 478 479 mpfr_init(correction); 480 dt_astro_sin(&correction, &omega); 481 mpfr_mul_d(correction, correction, -0.00017, GMP_RNDN); 482 483 { 484 int i; 485 for( i = 0; i < NTH_NEW_MOON_CORRECTION_ARGS_SIZE; i++ ) { 486 mpfr_t a, v, w, x, y, z; 487 mpfr_init_set_d(v, NTH_NEW_MOON_CORRECTION_ARGS[i][0], GMP_RNDN); 488 mpfr_init_set_d(w, NTH_NEW_MOON_CORRECTION_ARGS[i][1], GMP_RNDN); 489 mpfr_init_set_d(x, NTH_NEW_MOON_CORRECTION_ARGS[i][2], GMP_RNDN); 490 mpfr_init_set_d(y, NTH_NEW_MOON_CORRECTION_ARGS[i][3], GMP_RNDN); 491 mpfr_init_set_d(z, NTH_NEW_MOON_CORRECTION_ARGS[i][4], GMP_RNDN); 492 493 mpfr_mul(x, x, solar_anomaly, GMP_RNDN); 494 mpfr_mul(y, y, lunar_anomaly, GMP_RNDN); 495 mpfr_mul(z, z, moon_argument, GMP_RNDN); 496 497 mpfr_add(x, x, y, GMP_RNDN); 498 mpfr_add(x, x, z, GMP_RNDN); 499 dt_astro_sin(&x, &x); 500 501 mpfr_init(a); 502 mpfr_pow(a, E, w, GMP_RNDN); 503 504 mpfr_mul(a, a, v, GMP_RNDN); 505 mpfr_mul(a, a, x, GMP_RNDN); 506 mpfr_add( correction, correction, a, GMP_RNDN ); 507 508 mpfr_clear(a); 509 mpfr_clear(v); 510 mpfr_clear(w); 511 mpfr_clear(x); 512 mpfr_clear(y); 513 mpfr_clear(z); 514 } 515 } 516 517 { 518 int z; 519 mpfr_init_set_ui(additional, 0, GMP_RNDN); 520 for (z = 0; z < NTH_NEW_MOON_ADDITIONAL_ARGS_SIZE; z++) { 521 mpfr_t i, j, l; 522 mpfr_init_set_d(i, NTH_NEW_MOON_ADDITIONAL_ARGS[z][0], GMP_RNDN); 523 mpfr_init_set_d(j, NTH_NEW_MOON_ADDITIONAL_ARGS[z][1], GMP_RNDN); 524 mpfr_init_set_d(l, NTH_NEW_MOON_ADDITIONAL_ARGS[z][2], GMP_RNDN); 525 526 mpfr_mul(j, j, k, GMP_RNDN); 527 mpfr_add(j, j, i, GMP_RNDN); 528 dt_astro_sin(&j, &j); 529 mpfr_mul(l, l, j, GMP_RNDN); 530 531 mpfr_add(additional, additional, l, GMP_RNDN); 532 533 mpfr_clear(i); 534 mpfr_clear(j); 535 mpfr_clear(l); 536 } 537 } 538 539#ifdef ANNOYING_DEBUG 540#if (ANNOYING_DEBUG) 541mpfr_fprintf(stderr, 542 "correction = %.10RNf\nextra = %.10RNf\nadditional = %.10RNf\n", correction, extra, additional ); 543#endif 544#endif 545 mpfr_set(*result, approx, GMP_RNDN); 546 mpfr_add(*result, *result, correction, GMP_RNDN); 547 mpfr_add(*result, *result, extra, GMP_RNDN); 548 mpfr_add(*result, *result, additional, GMP_RNDN); 549 550 adjust_lunar_phase_to_zero( result ); 551 552 mpfr_clear(n); 553 mpfr_clear(k); 554 mpfr_clear(C); 555 mpfr_clear(approx); 556 mpfr_clear(E); 557 mpfr_clear(solar_anomaly); 558 mpfr_clear(lunar_anomaly); 559 mpfr_clear(moon_argument); 560 mpfr_clear(omega); 561 mpfr_clear(extra); 562 mpfr_clear(correction); 563 mpfr_clear(additional); 564 565 566 if (dt_astro_global_cache.cache_size == 0) { 567 dt_astro_global_cache.cache_size = 200000; 568 Newxz( dt_astro_global_cache.cache, dt_astro_global_cache.cache_size, mpfr_t * ); 569 } 570 571 if (dt_astro_global_cache.cache_size < n_int + 1) { 572 int new_size = dt_astro_global_cache.cache_size * 2; 573 Renew( dt_astro_global_cache.cache, new_size, mpfr_t *); 574 dt_astro_global_cache.cache_size = new_size; 575 } 576 577 { 578 mpfr_t *new_data; 579 Newxz( new_data, 1, mpfr_t ); 580 mpfr_init( *new_data ); 581 mpfr_set( *new_data, *result, GMP_RNDN ); 582 dt_astro_global_cache.cache[n_int] = new_data; 583#if (0) 584 PerlIO_printf(PerlIO_stderr(), "CACHE set for %d\n", n_int); 585#endif 586 } 587 588 return 1; 589} 590 591/* TODO: Check out errata 269 on 592 http://emr.cs.uiuc.edu/home/reingold/calendar-book/second-edition/errata.pdf 593*/ 594STATIC_INLINE 595int 596new_moon_before_from_moment(mpfr_t *result, mpfr_t *o_moment) { 597 mpfr_t t0, phi, big_n, radian, delta; 598 mpfr_t fullangle; 599 long n; 600 601 mpfr_init(t0); 602 nth_new_moon(&t0, 0); 603 604 mpfr_init(phi); 605 lunar_phase( &phi, o_moment ); 606 607 mpfr_init_set(delta, *o_moment, GMP_RNDN); 608 mpfr_sub(delta, delta, t0, GMP_RNDN); 609 mpfr_div_d(delta, delta, MEAN_SYNODIC_MONTH, GMP_RNDN); 610 611 mpfr_init_set_ui( fullangle, 360, GMP_RNDN ); 612 mpfr_init_set(radian, phi, GMP_RNDN); 613 mpfr_div(radian, radian, fullangle, GMP_RNDN ); 614 615 mpfr_init_set(big_n, delta, GMP_RNDN); 616 mpfr_sub(big_n, big_n, radian, GMP_RNDN); 617 mpfr_round(big_n, big_n); 618 619 n = mpfr_get_si(big_n, GMP_RNDN); 620 621 mpfr_clear(t0); 622 mpfr_clear(phi); 623 mpfr_clear(big_n); 624 mpfr_clear(radian); 625 mpfr_clear(delta); 626 mpfr_clear(fullangle); 627 628 { 629 int increment = 1; 630 mpfr_t last_good, tmp; 631 mpfr_init(tmp); 632 mpfr_init(last_good); 633 634 n = n - 1; 635 nth_new_moon( &last_good, n ); 636 while (increment) { 637 nth_new_moon( &tmp, n ); 638 if (mpfr_cmp( tmp, *o_moment ) < 0) { 639 n = n + 1; 640 mpfr_set( last_good, tmp, GMP_RNDN ); 641 } else { 642 increment = 0; 643 mpfr_set( *result, last_good, GMP_RNDN ); 644 } 645 } 646 647 mpfr_clear(tmp); 648 mpfr_clear(last_good); 649 } 650 651 return 1; 652} 653 654STATIC_INLINE 655int 656new_moon_after_from_moment(mpfr_t *result, mpfr_t *o_moment) { 657 mpfr_t t0, phi, big_n, radian, delta; 658 mpfr_t fullangle; 659 long n; 660 661 mpfr_init(t0); 662 nth_new_moon(&t0, 0); 663 664 mpfr_init(phi); 665 666 lunar_phase( &phi, o_moment ); 667 668 mpfr_init_set(delta, *o_moment, GMP_RNDN); 669 mpfr_sub(delta, delta, t0, GMP_RNDN); 670 mpfr_div_d(delta, delta, MEAN_SYNODIC_MONTH, GMP_RNDN); 671 672 mpfr_init_set_ui( fullangle, 360, GMP_RNDN ); 673 mpfr_init_set(radian, phi, GMP_RNDN); 674 mpfr_div(radian, radian, fullangle, GMP_RNDN ); 675 676 mpfr_init_set(big_n, delta, GMP_RNDN); 677 mpfr_sub(big_n, big_n, radian, GMP_RNDN); 678 mpfr_round(big_n, big_n); 679 680 n = mpfr_get_si(big_n, GMP_RNDN); 681 682 mpfr_clear(t0); 683 mpfr_clear(phi); 684 mpfr_clear(big_n); 685 mpfr_clear(radian); 686 mpfr_clear(delta); 687 mpfr_clear(fullangle); 688 689 nth_new_moon( result, n ); 690#if (0) 691mpfr_fprintf(stderr, 692 "got result = %.10RNf against %.10RNf\n", 693 result, 694 o_moment); 695#endif 696 697 { 698 /* if the result and moment are too close to each other, we 699 wrongly match the same time. 700 fractional parts below 0.0001 doesn't make any change to 701 the calculation, so we make that exception, and make sure 702 to think of those two as the same 703 */ 704 int loop_count = 0; 705 int cmp_result = -1; 706 mpfr_t delta; 707 mpfr_init(delta); 708 709 while (cmp_result <= 0) { 710 loop_count++; 711 cmp_result = mpfr_cmp( *result, *o_moment ); 712 if (cmp_result > 0) { 713 mpfr_dim(delta, *result, *o_moment, GMP_RNDN); 714 if (mpfr_cmp_d(delta, 0.001) <= 0) { 715 cmp_result = 0; 716 } 717 } 718 719 if (cmp_result <= 0) { 720 n = n + 1; 721 nth_new_moon( result, n ); 722 } 723 } 724 mpfr_clear(delta); 725 } 726 727 return 1; 728} 729 730#endif /** __DT_ASTRO_LUNAR_C__ **/ 731 732