```    //  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);
}
}
```