subroutine pnmatrix_j2000( jde, np )
c
c     Return the orthogonal rotation matrix which applies first
c     the precession, then the nutation from the mean equinox and
c     equator of J2000.0 to the true equinox and equator for a
c     given Julian Ephemeris (Dynamical)  date.
c
c     B. Knapp, 2000-09-18, 2001-04-24 (i*4)
c
C
C     RCS DATA
C     
C     $Header$
C     
C     $Log$
C
C
      implicit none
c
c     Input
      real*8 jde
c
c     Output
      real*8 np(3,3)
c
c     Constants
      real*8 J2000, JDE_J2000
      parameter (J2000=2000.0d0, JDE_J2000=2451545.d0)

c     Local
      integer*4 i, j, k
      real*8 ny, td, acc, p(3,3), n(3,3)
c
c     Externals
      external precess_matrix, nutate_matrix
c
      NY = (JDE-JDE_J2000)/365.25D0
      CALL PRECESS_MATRIX(J2000,J2000+NY,P)
      TD = (JDE-JDE_J2000)/365250.D0
      CALL NUTATE_MATRIX(TD,N)
      DO I=1,3
         DO J=1,3
            ACC = 0.D0
            DO K=1,3
               ACC = ACC+N(I,K)*P(K,J)
            ENDDO
            NP(I,J) = ACC
         ENDDO
      ENDDO
c
      return
      end