view release on metacpan or  search on metacpan

xs/dt_astro_lunar.inc  view on Meta::CPAN

            mpfr_clear(i);
            mpfr_clear(j);
            mpfr_clear(l);
        }
    }

#ifdef ANNOYING_DEBUG
#if (ANNOYING_DEBUG)
mpfr_fprintf(stderr,
    "correction = %.10RNf\nextra = %.10RNf\nadditional = %.10RNf\n", correction, extra, additional );
#endif
#endif
    mpfr_set(*result, approx, GMP_RNDN);
    mpfr_add(*result, *result, correction, GMP_RNDN);
    mpfr_add(*result, *result, extra, GMP_RNDN);
    mpfr_add(*result, *result, additional, GMP_RNDN);

    adjust_lunar_phase_to_zero( result );

    mpfr_clear(n);
    mpfr_clear(k);
    mpfr_clear(C);
    mpfr_clear(approx);
    mpfr_clear(E);
    mpfr_clear(solar_anomaly);
    mpfr_clear(lunar_anomaly);
    mpfr_clear(moon_argument);
    mpfr_clear(omega);
    mpfr_clear(extra);
    mpfr_clear(correction);
    mpfr_clear(additional);


    if (dt_astro_global_cache.cache_size == 0) {
        dt_astro_global_cache.cache_size = 200000;
        Newxz( dt_astro_global_cache.cache, dt_astro_global_cache.cache_size, mpfr_t * );
    }

    if (dt_astro_global_cache.cache_size < n_int + 1) {
        int new_size = dt_astro_global_cache.cache_size * 2;
        Renew( dt_astro_global_cache.cache, new_size, mpfr_t *);
        dt_astro_global_cache.cache_size = new_size;
    }

    {
        mpfr_t *new_data;
        Newxz( new_data, 1, mpfr_t );
        mpfr_init( *new_data );
        mpfr_set( *new_data, *result, GMP_RNDN );
        dt_astro_global_cache.cache[n_int] = new_data;
#if (0)
        PerlIO_printf(PerlIO_stderr(), "CACHE set for %d\n", n_int);
#endif
    }

    return 1;
}

/* TODO: Check out errata 269 on 
    http://emr.cs.uiuc.edu/home/reingold/calendar-book/second-edition/errata.pdf
*/
STATIC_INLINE
int
new_moon_before_from_moment(mpfr_t *result, mpfr_t *o_moment) {
    mpfr_t t0, phi, big_n, radian, delta;
    mpfr_t fullangle;
    long n;

    mpfr_init(t0);
    nth_new_moon(&t0, 0);

    mpfr_init(phi);
    lunar_phase( &phi, o_moment );

    mpfr_init_set(delta, *o_moment, GMP_RNDN);
    mpfr_sub(delta, delta, t0, GMP_RNDN);
    mpfr_div_d(delta, delta, MEAN_SYNODIC_MONTH, GMP_RNDN);

    mpfr_init_set_ui( fullangle, 360, GMP_RNDN );
    mpfr_init_set(radian, phi, GMP_RNDN);
    mpfr_div(radian, radian, fullangle, GMP_RNDN );

    mpfr_init_set(big_n, delta, GMP_RNDN);
    mpfr_sub(big_n, big_n, radian, GMP_RNDN);
    mpfr_round(big_n, big_n);

    n = mpfr_get_si(big_n, GMP_RNDN);

    mpfr_clear(t0);
    mpfr_clear(phi);
    mpfr_clear(big_n);
    mpfr_clear(radian);
    mpfr_clear(delta);
    mpfr_clear(fullangle);

    {
        int increment = 1;
        mpfr_t last_good, tmp;
        mpfr_init(tmp);
        mpfr_init(last_good);

        n = n - 1;
        nth_new_moon( &last_good, n );
        while (increment) {
            nth_new_moon( &tmp, n );
            if (mpfr_cmp( tmp, *o_moment ) < 0) {
                n = n + 1;
                mpfr_set( last_good, tmp, GMP_RNDN );
            } else {
                increment = 0;
                mpfr_set( *result, last_good, GMP_RNDN );
            }
        }
                
        mpfr_clear(tmp);
        mpfr_clear(last_good);
    }

    return 1;
}