def Lieskeprecangles(jd1, jd2):
    """
----------------------------------------------------------------------
Purpose:   Calculate IAU 1976 precession angles for a precession
           of epoch corresponding to Julian date jd1 to epoch corresponding
	   to Julian date jd2.
Input:     Julian date jd1 for start epoch
           Julian date jd2 for end epoch
Returns:   Angles zeta, z, theta in degrees
Reference: Lieske,J.H., 1979. Astron.Astrophys.,73,282.
           equations (6) & (7), p283.
Notes:     The ES (Explanatory Supplement to the Astronomical Almanac)
           lists for a IAU1976 precession from 1984, January 1d0h to J2000 
           the angles:  xi_a=368.9985, ze_a=369.0188 and th_a=320.7279
           Using the functions in this module, this can be calculated 
           by applying:

           jd1 = JD(1984,1,1) and 
           jd2 = JD(2000,1,1.5)
           print Lieskeprecangles(jd1, jd2) 

           which returns the equivalent values in degrees of:
           368.998509562 369.018811234 320.727934238 in arcsec.

           Lieske's fit belongs to the so called Quasi-Linear Types
           Precision (according to IAU SOFA):
             1960AD to 2040AD: < 0.1"
             1640AD to 2360AD: < 1"
             500BC to 3000AD: < 3"
             1200BC to 3900AD: > 10"
             < 4200BC or > 5600AD: > 100"
             < 6800BC or > 8200AD: > 1000"
----------------------------------------------------------------------
   """
    T = (jd1 - 2451545.0) / 36525.0
    t = (jd2 - jd1) / 36525.0
    d1 = 2306.2181
    d2 = 1.39656
    d3 = -0.000139
    d4 = 0.30188
    d5 = 0.000344
    d6 = 0.017998
    D1 = d1 + T * (d2 + T * d3)
    zeta_a = t * (D1 + t * (d4 + d5 * T + t * d6))
    d4 = 1.09468
    d5 = -6.6e-05
    d6 = 0.018203
    z_a = t * (D1 + t * (d4 + d5 * T + t * d6))
    d1 = 2004.3109
    d2 = -0.8533
    d3 = -0.000217
    d4 = -0.42665
    d5 = -0.000217
    d6 = -0.041833
    D1 = d1 + T * (d2 + T * d3)
    theta_a = t * (D1 + t * (d4 + d5 * T + t * d6))
    return zeta_a / 3600.0, z_a / 3600.0, theta_a / 3600.0
