def Newcombprecangles(epoch1, epoch2):
   """
----------------------------------------------------------------------
Purpose:    Calculate precession angles for a precession in FK4, using 
            Newcombs method (Woolard and Clemence angles)
Input:      Besselian start epoch1 and Besselian end epoch2
Returns:    Angles zeta, z, theta in degrees
Reference:  ES 3.214 p.106
Notes:      Newcomb's precession angles for old catalogs (FK4),
            see ES 3.214 p.106
            Input are Besselian epochs!
            Adopted accumulated precession angles from equator
            and equinox at B1950 to 1984 January 1d 0h according
            to ES (table 3.214.1, p 107) are:
            zeta=783.7092, z=783.8009 and theta=681.3883
            The Woolard and Clemence angles (derived in this routine)
            are:
            zeta=783.70925, z=783.80093 and theta=681.38830 (same ES table)
            This routine found (in seconds of arc):
            zeta,z,theta =  783.709246271 783.800934641 681.388298284
            for t1 = 0.1 and t2 = 0.133999566814
            using the lines in the example 
Examples:
            >>> b1 = 1950.0
            >>> b2 = celestial.epochs("F1984-01-01")[0]
            >>> print [x*3600 for x in celestial.Newcombprecangles(be1, be2)]
                [783.70924627097793, 783.80093464073127, 681.38829828393466]
----------------------------------------------------------------------
   """
   t1 = (epoch1-1850.0)/1000.0    #1000 tropical years
   t2 = (epoch2-1850.0)/1000.0
   tau = t2 - t1

   d0 = 23035.545; d1 = 139.720; d2 = 0.060; d3 = 30.240; d4 = -0.27; d5 = 17.995
   a0 = d0 + t1*(d1+d2*t1); a1 = d3 + d4*t1; a2 = d5
   zeta_a = tau*(a0+tau*(a1+tau*a2))

   d0 = 23035.545; d1 = 139.720; d2 = 0.060; d3 = 109.480; d4 = 0.39; d5 = 18.325
   a0 = d0 + t1*(d1+d2*t1); a1 = d3 + d4*t1; a2 = d5
   z_a = tau*(a0+tau*(a1+tau*a2))

   d0 = 20051.12; d1 = -85.29; d2 = -0.37; d3 = -42.65; d4 = -0.37; d5 = -41.80
   a0 = d0 + t1*(d1+d2*t1); a1 = d3 + d4*t1; a2 = d5
   theta_a = tau*(a0+tau*(a1+tau*a2))
   # Return values in degrees
   return zeta_a/3600.0, z_a/3600.0, theta_a/3600.0