#ifndef __DT_ASTRO_SOLAR_C__
#define __DT_ASTRO_SOLAR_C__
#include "dt_astro.h"

/* [1] Table 12.1 p.183 (zero-padded to align veritcally) */
#define SOLAR_LONGITUDE_ARGS_SIZE 49
static const double SOLAR_LONGITUDE_ARGS[59][3] = {
    /*      X          Y             Z        */
    /* left side of table 12.1                */
    { 403406, 270.54861,      0.9287892 },
    { 119433,  63.91854,  35999.4089666 },
    {   3891, 317.84300,  71998.2026100 },
    {   1721, 240.05200,  36000.3572600 },
    {    350, 247.23000,  32964.4678000 },
    {    314, 297.82000, 445267.1117000 },
    {    242, 166.79000,      3.1008000 },
    {    158,   3.50000,    -19.9739000 },
    {    129, 182.95000,   9038.0293000 },
    {     99,  29.80000,  33718.1480000 },
    {     86, 249.20000,  -2280.7730000 },
    {     72, 257.80000,  31556.4930000 },
    {     64,  69.90000,   9037.7500000 },
    {     38, 197.10000,  -4444.1760000 },
    {     32,  65.30000,  67555.3160000 },
    {     28, 341.50000,  -4561.5400000 },
    {     27,  98.50000,   1221.6550000 },
    {     24, 110.00000,  31437.3690000 },
    {     21, 342.60000, -31931.7570000 },
    {     18, 256.10000,   1221.9990000 },
    {     14, 242.90000,  -4442.0390000 },
    {     13, 151.80000,    119.0660000 },
    {     12,  53.30000,     -4.5780000 },
    {     10, 205.70000,    -39.1270000 },
    {     10, 146.10000,  90073.7780000 },
    /* right side of table 12.1               */
    { 195207, 340.19128,  35999.1376958 },
    { 112392, 331.26220,  35998.7287385 },
    {   2819,  86.63100,  71998.4403000 },
    {    660, 310.26000,  71997.4812000 },
    {    334, 260.87000,    -19.4410000 },
    {    268, 343.14000,  45036.8840000 },
    {    234,  81.53000,  22518.4434000 },
    {    132, 132.75000,  65928.9345000 },
    {    114, 162.03000,   3034.7684000 },
    {     93, 266.40000,   3034.4480000 },
    {     78, 157.60000,  29929.9920000 },
    {     68, 185.10000,    149.5880000 },
    {     46,   8.00000, 107997.4050000 },
    {     37, 250.40000,    151.7710000 },
    {     29, 162.70000,  31556.0800000 },
    {     27, 291.60000, 107996.7060000 },
    {     25, 146.70000,  62894.1670000 },
    {     21,   5.20000,  14578.2980000 },
    {     20, 230.90000,  34777.2430000 },
    {     17,  45.30000,  62894.5110000 },
    {     13, 115.20000, 107997.9090000 },
    {     13, 285.30000,  16859.0710000 },
    {     10, 126.60000,  26895.2920000 },
    {     10,  85.90000,  12297.5360000 }
};

STATIC_INLINE
int
solar_longitude( mpfr_t *result, mpfr_t *moment ) {
    mpfr_t C, fugly, A, N, longitude, fullangle;

    mpfr_init(C);
    julian_centuries(&C, moment);

#ifdef ANNOYING_DEBUG
#if (ANNOYING_DEBUG)
mpfr_fprintf(stderr,
    "julian_centuries = %.10RNf\n", C);
#endif
#endif
    {
        int i;
        mpfr_init_set_ui(fugly, 0, GMP_RNDN);
        for( i = 0; i < SOLAR_LONGITUDE_ARGS_SIZE; i++ ) {
            mpfr_t a, b, c;
            mpfr_init_set_d( a, SOLAR_LONGITUDE_ARGS[i][0], GMP_RNDN );
            mpfr_init_set_d( b, SOLAR_LONGITUDE_ARGS[i][1], GMP_RNDN );
            mpfr_init_set_d( c, SOLAR_LONGITUDE_ARGS[i][2], GMP_RNDN );
            mpfr_mul(c, c, C, GMP_RNDN);
            mpfr_add(b, b, c, GMP_RNDN);
            dt_astro_sin(&b, &b);
            mpfr_mul(a, a, b, GMP_RNDN);
            mpfr_add(fugly, fugly, a, GMP_RNDN);
            mpfr_clear(a);
            mpfr_clear(b);
            mpfr_clear(c);
        }
    }

    {
        mpfr_t b, c;
        mpfr_init_set_d(longitude, 282.7771834, GMP_RNDN);
        mpfr_init_set_d(b, 36000.76953744, GMP_RNDN);
        mpfr_mul(b, b, C, GMP_RNDN);
        mpfr_init_set_d(c, 0.000005729577951308232, GMP_RNDN);
        mpfr_mul(c, c, fugly, GMP_RNDN);
        mpfr_add(longitude, longitude, b, GMP_RNDN);
        mpfr_add(longitude, longitude, c, GMP_RNDN);
        mpfr_clear(b);
        mpfr_clear(c);
    }

    mpfr_init(A);
    aberration(&A, moment);

    mpfr_init(N);
    nutation(&N, moment);

#ifdef ANNOYING_DEBUG
#if (ANNOYING_DEBUG)
mpfr_fprintf(stderr,
    "longitude = %.10RNf\naberration = %.10RNf\nnutation = %.10RNf\n",
    longitude,
    A,
    N);
#endif
#endif

    mpfr_set( *result, longitude, GMP_RNDN);
    mpfr_add( *result, *result, A, GMP_RNDN );
    mpfr_add( *result, *result, N, GMP_RNDN );

    mpfr_init_set_ui( fullangle, 360, GMP_RNDN );
#ifdef ANNOYING_DEBUG
#if (ANNOYING_DEBUG)
mpfr_fprintf(stderr, "(solar) mod(%.10RNf) = ", *result );
#endif
#endif
    dt_astro_mod( result, result, &fullangle );
#ifdef ANNOYING_DEBUG
#if (ANNOYING_DEBUG)
mpfr_fprintf(stderr, "%.10RNf\n", *result );
#endif
#endif

    mpfr_clear(A);
    mpfr_clear(N);
    mpfr_clear(C);
    mpfr_clear(longitude);
    mpfr_clear(fullangle);
    mpfr_clear(fugly);

    return 1;
}

STATIC_INLINE
int
estimate_prior_solar_longitude(mpfr_t *result, mpfr_t *moment, mpfr_t *phi) {
    mpfr_t tau, delta;

    mpfr_init_set(tau, *moment, GMP_RNDN);
    mpfr_init(delta);

    {
        mpfr_t tmp;
        mpfr_t fullangle;
        mpfr_init(tmp);
        mpfr_init_set_ui(fullangle, 360, GMP_RNDN);

        solar_longitude( &tmp, moment );
        mpfr_sub( tmp, tmp, *phi, GMP_RNDN );
        dt_astro_mod( &tmp, &tmp, &fullangle );
        mpfr_mul_d( tmp, tmp, SOLAR_YEAR_RATE, GMP_RNDN );

        mpfr_sub( tau, tau, tmp, GMP_RNDN );

        {
            mpfr_t tau_lon;

            mpfr_init(tau_lon);
            /* tau_lon = solar_longitude(tau_lon_arg) */
            solar_longitude( &tau_lon, &tau );

            mpfr_sub( tau_lon, tau_lon, *phi, GMP_RNDN );
            mpfr_add_ui( tau_lon, tau_lon, 180, GMP_RNDN );

            dt_astro_mod( &delta, &tau_lon, &fullangle );

            mpfr_clear(tau_lon);
        }

        mpfr_sub_ui( delta, delta, 180, GMP_RNDN );

        mpfr_clear(tmp);
        mpfr_clear(fullangle);
    }

    mpfr_mul_d( delta, delta, SOLAR_YEAR_RATE, GMP_RNDN );
    mpfr_sub( tau, tau, delta, GMP_RNDN );
    if (mpfr_cmp( *moment, tau ) > 0) {
        mpfr_set( *result, tau, GMP_RNDN );
    } else {
        mpfr_set( *result, *moment, GMP_RNDN);
    }
        
    mpfr_clear(tau);
    mpfr_clear(delta);

    return 1;
}

#define SOLAR_LONGITUDE_ALLOWED_DELTA 0.00000000000000000000000000000001

STATIC_INLINE
int
__solar_longitude_mu(mpfr_t *lo, mpfr_t *hi) {
    int result = 0;
    mpfr_t delta;
    mpfr_init_set(delta, *lo, GMP_RNDN);
    mpfr_sub(delta, delta, *hi, GMP_RNDN);
    mpfr_abs(delta, delta, GMP_RNDN);

    if (mpfr_cmp_d(delta, SOLAR_LONGITUDE_ALLOWED_DELTA) < 0) {
        result = 1;
    }
    mpfr_clear(delta);
    return result;
}

STATIC_INLINE
int
__solar_longitude_phi(mpfr_t *x, void *args, int n_args) {
    int result = 0;
    mpfr_t phi, lon, fullangle;

    PERL_UNUSED_VAR(n_args);
    
    mpfr_init_set_ui(fullangle, 360, GMP_RNDN);
    
    mpfr_init( lon );
    mpfr_init_set( phi, ((mpfr_t *) args)[0], GMP_RNDN );
    
    solar_longitude( &lon, x );

    mpfr_sub( lon, lon, phi, GMP_RNDN );
    dt_astro_mod( &lon, &lon, &fullangle );

    if ( mpfr_cmp_ui( lon, 180 ) < 0 ) {
        result = 1;
    }

    mpfr_clear(lon);
    mpfr_clear(phi);
    mpfr_clear(fullangle);
    return result;
}

STATIC_INLINE
int
solar_longitude_before( mpfr_t *result, mpfr_t *moment, mpfr_t *phi ) {
    mpfr_t tau, l, u;
    mpfr_init(tau);

    estimate_prior_solar_longitude( &tau, moment, phi );
    mpfr_init_set(l, tau, GMP_RNDN);
    mpfr_sub_ui(l, l, 5, GMP_RNDN);
    mpfr_init_set(u, tau, GMP_RNDN);
    mpfr_add_ui(u, u, 5, GMP_RNDN);
    
    if (mpfr_cmp( *moment, u ) < 0) {
        mpfr_set( u, *moment, GMP_RNDN );
    }

    __binary_search( result, &l, &u, __solar_longitude_phi, (void *) phi, 1, __solar_longitude_mu );

    mpfr_clear(tau);
    mpfr_clear(l);
    mpfr_clear(u);

    return 1;
}

STATIC_INLINE
int
solar_longitude_after( mpfr_t *result, mpfr_t *moment, mpfr_t *phi ) {
    mpfr_t tau, n, l, u;

    mpfr_init_set_d( n, SOLAR_YEAR_RATE, GMP_RNDN );

    {
        mpfr_t lon, x, fullangle;

        mpfr_init_set_ui(fullangle, 360, GMP_RNDN);
        mpfr_init(lon);
        solar_longitude( &lon, moment );
        mpfr_init_set( x, *phi, GMP_RNDN );
        mpfr_sub(x, x, lon, GMP_RNDN );
        dt_astro_mod( &x, &x, &fullangle );
        mpfr_mul( n, n, x, GMP_RNDN );

        mpfr_clear(lon);
        mpfr_clear(x);
        mpfr_clear(fullangle);
    }

    mpfr_init_set( tau, *moment, GMP_RNDN );
    mpfr_add( tau, tau, n, GMP_RNDN );

    {
        mpfr_t x;
        mpfr_init_set( x, tau, GMP_RNDN );
        mpfr_sub_ui( x, x, 5, GMP_RNDN );
        if ( mpfr_cmp (*moment, x) > 0 ) {
            mpfr_init_set( l, *moment, GMP_RNDN );
        } else {
            mpfr_init_set( l, x, GMP_RNDN );
        }

        mpfr_clear(x);
    }

    mpfr_init_set(u, tau, GMP_RNDN);
    mpfr_add_ui(u, u, 5, GMP_RNDN);

    __binary_search( result, &l, &u, __solar_longitude_phi, (void *) phi, 1, __solar_longitude_mu );

    mpfr_clear(tau);
    mpfr_clear(n);
    mpfr_clear(l);
    mpfr_clear(u);
    
    return 1;
}

#endif /** __DT_ASTRO_SOLAR_C__ **/