SUBROUTINE HCPVD( P, T, PV, IERR )
C
C     CALCULATE THE TRUE HELIOCENTRIC POSITION AND VELOCITY OF PLANET P
C     AT TIME T, USING THE VSOP87D THEORY (P BRETAGNON AND G FRANCOU,
C     ASTRON. ASTROPHYS., 202 (1988), 309-315).
C
C     B. KNAPP, 1998-05-07, 1999-05-11, 2001-04-24 (i*4)
C
C     P IS 0,...,9 FOR SUN,...,PLUTO.  T IS DYNAMIC TIME IN JULIAN
C     MILLENIA ELAPSED SINCE J2000 (I.E., T = (JDE-2451545)/365250).
C     PV(1) AND PV(2), THE HELIOCENTRIC LONGITUDE ("L") AND LATITUDE
C     ("B"), RESPECTIVELY, ARE RETURNED AS RADIANS, AND PV(3), THE
C     HELIOCENTRIC RADIUS ("R"), IS RETURNED AS AU.  THE RESPECTIVE
C     VELOCITIES (PV(4) TO PV(6)) ARE RETURNED AS RADIANS/DAY OR AU/DAY.
C
C     NOTE: THIS ROUTINE MUST BE COMPILED WITH A FORTRAN/77 OR LATER
C     COMPILER, AS THE DO 2 LOOP BELOW HAS DEGENERATE CASES
C
C
C     RCS DATA
C     
C     $Header$
C     
C     $Log$
C
C
      IMPLICIT NONE
C
C     INPUT
      INTEGER*4 P
      REAL*8 T
C
C     OUTPUT
      REAL*8 PV(6)
      INTEGER*4 IERR
C
C     CONSTANTS
      REAL*8 TWOPI, JM
      PARAMETER (TWOPI = 2.D0*3.14159265358979324D0)
      PARAMETER (JM=365250.D0)
C
C     GET COSINE COEFFICIENTS FROM BLOCK DATA SUBPROGRAM VSOP87
      INTEGER*4 NDATA
      PARAMETER (NDATA=31577)
      INTEGER*4 OFFSET(0:5,3,8), NTERMS(0:5,3,8)
      REAL*8 COEFFS(3,NDATA)
      COMMON /VSOP87_COEFFS/ OFFSET, NTERMS, COEFFS
C
C     LOCAL VARIABLES
      INTEGER*4 I,J,K
      REAL*8 Q,R,U,V,W,X,Y,CU,SU
C
C     EXTERNAL FUNCTIONS, SUBROUTINES
      EXTERNAL VSOP87, PLUTO
C
C     INITIALIZE (HELIOCENTRIC PV OF SUN)
      DO 1 J=1,6
         PV(J) = 0.0
    1 CONTINUE
      IF ( P .LT. 0 .OR. 9 .LT. P ) THEN
         IERR = 1
         RETURN
      ENDIF
      IERR = 0
      IF ( P .EQ. 0 ) RETURN
      IF ( P .EQ. 9 ) THEN
         CALL PLUTO( T, PV, IERR )
         RETURN
      ENDIF
C
C     COMPUTE THE POSITION AND VELOCITY
      DO 4 J=1,3
         V = 0.D0
         W = 0.D0
         Y = 0.D0
         DO 3 I=5,0,-1
             X = 0.D0
             Q = 0.D0
             R = 0.D0
             DO 2 K=OFFSET(I,J,P)+NTERMS(I,J,P)-1,OFFSET(I,J,P),-1
                U = MOD( COEFFS(2,K)+COEFFS(3,K)*T, TWOPI )
                CU = COS( U )
                SU = SIN( U )
                X = X + COEFFS(1,K)*CU
                Q = Q + DBLE(I)*COEFFS(1,K)*CU
                R = R - COEFFS(1,K)*COEFFS(3,K)*SU
    2        CONTINUE
             Y = Y*T + X
             IF ( I .GT. 0 ) THEN
                V = V*T + Q
             ENDIF
             W = W*T + R
    3    CONTINUE
         PV(J) = Y
         PV(J+3) = (V + W)/JM
    4 CONTINUE
C
C     PUT LONGITUDE IN RANGE 0..2*PI
      PV(1) = MOD( PV(1), TWOPI )
      IF ( PV(1) .LT. 0.D0 ) PV(1) = PV(1)+TWOPI
C
      RETURN
      END