subroutine ephem_table( jd, p, posnid, framid, coorid, pv, dpv,
     &   ierr)
!
C     GIVEN A JULIAN DATE  JD  AND A SOLAR SYSTEM BODY NUMBER  P
C     (0=SUN, 1=MERCURY, 2=VENUS, 3=EARTH, 4=MARS, 5=JUPITER, 6=SATURN,
C     7=URANUS, 8=NEPTUNE, 9=PLUTO, 10=MOON) THIS ROUTINE COMPUTES
C     THE TRUE GEOMETRIC (POSNID=1) OR APPARENT (POSNID=2) POSITION
C     AND VELOCITY OF THE BODY, IN THE EQUATORIAL (FRAMID=1), ECLIPTIC
C     (FRAMID=2), OR HELIOCENTRIC (FRAMID=3) REFERENCE FRAME, AND IN
C     SPHERICAL (COORID=1) OR RECTANGULAR (COORID=2) COORDINATES.
C
C     B. KNAPP, 1994-06-14, 1999-09-28, 2001-04-24 (i*4)
C
C     THE ELEMENTS OF THE OUTPUT 6-VECTOR PV ARE AS FOLLOWS:
C
C        I   SPHERICAL                   RECTANGULAR
C
C        1   RA OR L (DEG, 0 TO 360)     X (AU)
C        2   DEC OR B (DEG, -90 TO +90)  Y (AU)
C        3   R (AU)                      Z (AU)
C        4   DRA/DT OR DL/DT (DEG/DAY)   DX/DT (AU/DAY)
C        5   DDEC/DT OR DB/DT (DEG/DAY)  DY/DT (AU/DAY)
C        6   DR/DT (AU/DAY)              DZ/DT (AU/DAY)
C
C     IERR IS RETURNED NON-ZERO FOR INVALID INPUT(S)
C
C
C     RCS DATA
C
C     $Header$
C
C     $Log$
C
C
      IMPLICIT NONE
C
C     INPUT
      REAL*8 JD
      INTEGER*4 P, POSNID, FRAMID, COORID
C
C     OUTPUT
      REAL*8 PV(6), DPV(6)
      INTEGER*4 IERR
!
!     Constants
      integer*4 N_SUN, N_MER, N_VEN, N_EAR, N_MAR, 
     &   N_JUP, N_SAT, N_URA, N_NEP, N_PLU, N_MOON, N_SUM
      parameter (N_SUN=24, N_MER=90, N_VEN=36, N_EAR=N_SUN, N_MAR=24, 
     &   N_JUP=12, N_SAT=12, N_URA=12, N_NEP=12, N_PLU=12, N_MOON=300)
      parameter (N_SUM = N_SUN + N_MER + N_VEN + N_EAR + N_MAR + 
     &   N_JUP + N_SAT + N_URA + N_NEP + N_PLU + N_MOON)
!
!     Local
      integer*4 n(0:10) /N_SUN, N_MER, N_VEN, N_EAR, N_MAR, N_JUP,
     &   N_SAT, N_URA, N_NEP, N_PLU, N_MOON/, offset(0:10)
      integer*4 i, j, k, ndx_jk, blksiz
      real*8 dd(0:24*N_SUM-1), y(-1:N_MOON+1,6), pvi(6), jd0, t
      save dd, n, offset, jd0

!     This version works only for apparent equatorial rectangular
!     coordinates
      if (posnid .ne. 2 .or. framid .ne. 1 .or. coorid .ne. 2) then
         ierr = -99
         return
      endif

!     We will build a divided difference table that will cover a
!     period of 30 hours centered on UT noon (an integer JD).

!     Init?
      if (jd-jd0 .ge. 0.625d0 .or. jd0-jd .gt. 0.625d0) then

         jd0 = int(jd+0.5)
         write(*,*) ' Initializing JD ', jd0
!
!        Construct the offset table
         offset(0) = 0
         do k=1,10
            offset(k) = offset(k-1) + 24*n(k-1)
         enddo

!        Build a divided differences table for each body
         do k=0,10

!           Collect the pv data at n+3 points
            do i=-1,n(k)+1
               t = jd0 - 0.625d0 + 1.25d0*dble(i)/dble(n(k))
               call ephem(t, k, posnid, framid, coorid, pvi, ierr)
               if (ierr .ne. 0) return
               do j=1,6
                  y(i,j) = pvi(j)
               enddo
            enddo

!           Build 6 divided difference tables, one each for x, y, z,
!           xdot, ydot, and zdot.
            blksiz = 4*n(k)
            do j=1,6
               ndx_jk = offset(k) + (j-1)*blksiz
               call divdifftable(n(k), y(-1,j), dd(ndx_jk))
            enddo
         enddo
      endif
!
!     Now evaluate the interpolating polynomials for each of
!     the six position elements, for the given JD
      t = (jd - jd0 + 0.625d0)*dble(n(p))/1.25d0
      blksiz = 4*n(p)
      do j=1,6
         ndx_jk = offset(p) + (j-1)*blksiz
         call divdiffeval(n(p), dd(ndx_jk), t, pv(j), dpv(j))
      enddo
!
      return
      end