      SUBROUTINE PRECESS_MATRIX(J0,J1,P)
C
C     OBTAIN THE 3 X 3 ROTATION MATRIX WHICH WILL PRECESS AN OBJECT
C     (WHOSE POSITION IS EXPRESSED IN RECTANGULAR EQUATORIAL
C     COORDINATES) FROM FROM JULIAN INITIAL EPOCH J0 TO JULIAN FINAL
C     EPOCH J1 (IN THE FK5 SYSTEM).  EPOCHS ARE GREGORIAN YEAR AND
C     FRACTION.
C
C     B. KNAPP, 1993-06-22
C
C     BASED ON THE ALGORITHM OF J. H. LIESKE ET AL, ASTRONOMY AND
C     ASTROPHYSICS, VOL 58, PP 1-16 (1977), AS PRESENTED IN SKY AND
C     TELESCOPE, VOL 82, PP 408-410 (OCT 1991).  NOTE: COEFFICIENT
C     H5 CHANGED FROM 3.45 TO 3.44--SEE J. H. LIESKE, ASTRONOMY AND
C     ASTROPHYSICS, VOL 73, PP 282-284 (1979).
C
C
C
C     RCS DATA
C     
C     $Header$
C     
C     $Log$
C
C
      IMPLICIT NONE
C
C     INPUT
      REAL*8 J0,J1
C
C     OUTPUT
      REAL*8 P(3,3)
C
C     LOCAL
C
      REAL*8 PI,D2R
      PARAMETER (PI=3.14159265358979324D0,D2R=PI/180.D0)
C
      REAL*8 H1,H2,H3,H4,H5,H6,K1,K2,K3,L1,L2,L3,L4,L5,L6
      PARAMETER (H1=2306.2181D0, H2=1.39656D0, H3=-1.39D-4,
     &   H4=0.30188D0, H5=-3.44D-4, H6=1.7998D-2, K1=1.09468D0,
     &   K2=6.6D-5, K3=1.8203D-2, L1=2004.3109D0, L2=-0.8533D0,
     &   L3=-2.17D-4, L4=-0.42665D0, L5=-2.17D-4, L6=-4.1833D-2)
C
      REAL*8 T0,T1,W,ZT,ZD,TH,S1,S2,S3,C1,C2,C3
C
C
      T0 = (J0-2000.D0)/100.D0
      T1 = (J1-J0)/100.D0
      W = H1+(H2+H3*T0)*T0
      ZT = (W+(H4+H5*T0+H6*T1)*T1)*T1/3600.D0*D2R
      ZD = (W+(K1+K2*T0+K3*T1)*T1)*T1/3600.D0*D2R
      W = L1+(L2+L3*T0)*T0
      TH = (W+(L4+L5*T0+L6*T1)*T1)*T1/3600.D0*D2R
C
      S1 = SIN(ZT)
      C1 = COS(ZT)
      S2 = SIN(ZD)
      C2 = COS(ZD)
      S3 = SIN(TH)
      C3 = COS(TH)
C
C     PRECESSION MATRIX:
      P(1,1) =  C1*C3*C2-S1*S2
      P(1,2) = -S1*C3*C2-C1*S2
      P(1,3) = -S3*C2
      P(2,1) =  C1*C3*S2+S1*C2
      P(2,2) = -S1*C3*S2+C1*C2
      P(2,3) = -S3*S2
      P(3,1) =  C1*S3
      P(3,2) = -S1*S3
      P(3,3) =  C3
C
      RETURN
      END
