//  GLCDIV.C
    //  Program to calculate GLC divergence
    //  Author: Peter Meyer
    //  Last mod.: 1998-11-17
    
    #include <STDIO.H>
    #include <MATH.H>
    
    #include <L_DATE.H>
    
    #define SECONDS_PER_DAY    (86400L)        //  Length of fixed-length day in atomic seconds
    #define JDN_2000           (2451545L)      //  Julian Day Number of 2000-01-01 G
    #define JULIAN_CENTURY     (36525L)        //  Julian century in days
    #define SOLAR_YEAR         (365.2422)      //  Length of solar year
    #define SYNODIC_MONTH_2000 (29.5305888531) //  Synodic month at 2000-01-01 G
    #define AV_GLC_MONTH       (29.53058880)   //  Average length of GLC month
    #define C1                 (2.1621E-7)
    #define C2                 (3.64E-10)
    #define K                  (3.42238E-7)    //  The rate at which the Earth is slowing
                                               //  in its rotation in revolutions/second^2
    /*-----------*/
    void main(void)
    {
    int year;
    long jd, yrs;
    double s, t, d, rsm, diff;
    Date dt;
    
    dt.day = dt.month = 1;
    dt.calendar = 'G';
    
    for ( year=-1000; year <=5000; year+=500 )
        {
        //  Get JDN for 1 January of year
        dt.year = year;
        jd = euro_date_to_jdn(&dt);
        //  S = 29.5305888531 + 0.00000021621 T - 3.64E-10 T^2 (days)
        //  where S is the synodic month and
        //  T = (JD - 2,451,545) / 36,525
        t = ( (double)jd - JDN_2000 ) / JULIAN_CENTURY;
        s = SYNODIC_MONTH_2000 + C1*t + C2*t*t;
    
        //  D = ( K * T * 36,525 + s ) / s
        d = ( K * t * JULIAN_CENTURY + SECONDS_PER_DAY ) / SECONDS_PER_DAY;
    
        rsm = s / d;
        diff = rsm - AV_GLC_MONTH;
        yrs = (long)( ( (1/fabs(diff)) * AV_GLC_MONTH ) / SOLAR_YEAR );
    
        printf("%6d%13.8f%13.9f%13.8f%13.8f%9ld\n",
            year,s,d,rsm,diff,yrs);
        }
    }