/**
 * A super quick and dirty implementation of conversion routines between Julian 
 * Day and Gregorian Calendar Dates.  The main() function demonstrates the 
 * conversion between the two is consistent.  Further validation was done be 
 * comparing against a bunch of javascript conversion utilities found online.
 *
 * The algorithms themselves are uncommented and mysterious.  They were taken 
 * from a website by Peter Baum:
 * http://www.vsg.cape.com/~pbaum/date/date0.htm
 *
 * This could falls into the do-whatever-you-want-with-it category.
 * 
 * Cheers,
 *     Niek Sanders
 *     http://www.cis.rit.edu/~njs8030/
 * 
 **/
// Required by the conversion routines
#include <cmath>

// Used in the demo main() program and its helpers
#include <iostream>
#include <ctime>
#include <cstdlib>
#include <cassert>


/**
 * Get the Julian day corresponding to the supplied Gregorian date.  Note that 
 * the supplied date should be in UT.  For instance, 2 PM on the 14th is 
 * 14+(14/24)=14.583
 *
 * \param D          The (fractional) day.  First day is 1.
 * \param M          The month.  January is 1.
 * \param Y          The year.
 *
 * \return           The Julian day.
 */
double getJulianDay( double D, int M, int Y ) {

    if ( M < 3 ) {
        M += 12;
        Y -= 1;
    }

    double F = -122 + static_cast<int>( 30.6 * ( M + 1 ));
    return D + F + 365 * Y + std::floor( Y / 4.0 ) - std::floor( Y / 100.0 ) + 
           std::floor( Y / 400.0 ) + 1721118.5;
}


/**
 * Get the modified Julian day (MJD).  The supplied date should be in UT.  Note 
 * that MJD 0.0 corresponds to Nov 17.0, 1858.
 *
 * \param D          The (fractional) day.  First day is 1.
 * \param M          The month.  January is 1.
 * \param Y          The year.
 *
 * \return           The modified Julian day.
 */
double getModifiedJulianDay( double D, int M, int Y ) {

    return getJulianDay( D, M, Y ) - 2400000.5;
}


/**
 * Get the Gregorian calendar date corresponding to a given Julian Day.
 *
 * \param J          The Julian Day.
 * \param D          The computed (fractional) day.  First day is 1.
 * \param M          The computed month.  January is 1.
 * \param Y          The computed year.
 */
void getGregorianDate( double J, double& D, int& M, int& Y ) {

    double Z = std::floor( J - 1721118.5 );
    double R = J - 1721118.5 - Z;
    double G = Z - 0.25;
    double A = std::floor( G / 36524.25 );
    double B = A - std::floor( A / 4.0 );
    Y = static_cast<int>( std::floor( (B + G) / 365.25 ) );
    double C = B + Z - std::floor( 365.25 * Y );
    M = static_cast<int>(( 5.0 * C + 456.0 ) / 153.0 );
    D = C - static_cast<int>( 30.6 * M - 91.4 ) + R;
    if ( M > 12 ) {
        Y += 1;
        M -= 12;
    }
}




double drnd( void ) {
    double aa = (double)::random();
    return aa / (double) 0x7FFFFFFFU;
}

int main( int argc, char* argv[] ) {

    // Seed the random number generators from the system time
    srand48( time( 0 ) );
    srand( time( 0 ) );

    // Generate ten million random Julian Days between year 1700 and 2200.
    // Check if we get same value after converting to Gregorian date and back 
    // again.
    static const double jd1700 = 2341972;
    static const double jd2200 = 2524593;

    for ( int i = 0; i < 10000000; i++ ) {
        double testJ = jd1700 + (drnd() * (jd2200 - jd1700));

        double D = 0.0;
        int M = 0, Y = 0;

        getGregorianDate( testJ, D, M, Y );
        double resultJ = getJulianDay( D, M, Y );

        if ( std::fabs( testJ - resultJ ) > 1E-6 ) {
            std::cerr << "ERROR!" << std::endl;
            std::cerr << "    iteration -->" << i << std::endl;
            std::cerr.precision( 12 );
            std::cerr << "   testing JD -->" << testJ << std::endl;
            std::cerr.precision( 12 );
            std::cerr << "  Gregorian D -->" << D << std::endl;
            std::cerr << "  Gregorian M -->" << M << std::endl;
            std::cerr << "  Gregorian Y -->" << Y << std::endl;
            std::cerr.precision( 12 );
            std::cerr << "  returned JD -->" << resultJ << std::endl;
            exit( 1 );
        }
    }

    return 0;
}