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 = (Current epoch - 1 jan, 2000, 12h noon)
   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))
   # d1 = 2306.2181
   # d2 = 1.39656
   # d3 = -0.000139
   d4 = 1.09468
   d5 = -0.000066
   d6 = 0.018203
   z_a = t*(D1 + t*((d4+d5*T) + t*d6))
   d1 = 2004.3109
   d2 = -0.85330
   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 values in degrees
   return zeta_a/3600.0, z_a/3600.0, theta_a/3600.0