SUBROUTINE EPHEM( JD, P, POSNID, FRAMID, COORID, PV, IERR )
C
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)
      INTEGER*4 IERR
C
C     CONSTANTS
      INTEGER*4 SUN, EARTH, MOON, TO_SPH, TO_XYZ
      PARAMETER (SUN=0, EARTH=3, MOON=10, TO_SPH=-1, TO_XYZ=+1)
      REAL*8 PI, TWOPI, D2R, C, AU, JM, AU2JM
      PARAMETER (PI=3.14159265358979324D0, TWOPI=2.D0*PI, D2R=PI/180.D0)
      PARAMETER (C=299792.458D0, AU=149597870.66D0, JM=365250.D0,
     &   AU2JM = AU/C/86400.D0/JM)
C
C     LOCAL VARIABLES
      REAL*8 JDE, T0, T, TAU, DPSI, DEPS, DDPSI, DDEPS, EPS,
     &   PVE(6), PVM(6)
      INTEGER*4 PLANET, IERR_P, IERR_E
      INTEGER*4 I, J, MODE
C
C     EXTERNAL FUNCTIONS, SUBROUTINES
      REAL*8 UT2TD, OBLIQ
      EXTERNAL UT2TD, HCPVD, FK5, NUTATE, EQUECL, OBLIQ, MOONPV, SPHXYZ
C
C
C     INITIALIZE
      DO 1, J=1,6
         PV(J) = 0.0D0
 1    CONTINUE
C
      IF ( P      .LT. 0 .OR. 10 .LT. P      .OR.
     &     POSNID .LT. 1 .OR.  2 .LT. POSNID .OR.
     &     FRAMID .LT. 1 .OR.  3 .LT. FRAMID .OR.
     &     COORID .LT. 1 .OR.  2 .LT. COORID ) THEN
         IERR = 1
         RETURN
      ENDIF
      IERR = 0
C
C     THE RESULT IS EXACTLY ZERO IN TWO CASES
      IF ( P .EQ. SUN   .AND. FRAMID .EQ. 3 .OR.
     &     P .EQ. EARTH .AND. FRAMID .NE. 3 ) RETURN
C
C     GET DYNAMIC TIME (IN JULIAN MILLENIA SINCE J2000)
      JDE = UT2TD(JD)
      T0 = (JDE-2451545.D0)/365250.D0
C
C     IF P IS THE MOON THEN "PLANET" IS THE EARTH
      IF ( P .EQ. MOON ) THEN
         PLANET = EARTH
      ELSE
         PLANET = P
      ENDIF
C
      TAU = 0.0
      DO 5 I=1,POSNID
         T = T0 - TAU
         IF ( P .EQ. MOON .AND. FRAMID .NE. 3 ) THEN
            DO 2 J=1,6
               PV(J) = 0.D0
 2          CONTINUE
         ELSE
            CALL HCPVD( PLANET, T, PV, IERR_P )
            IERR = MAX( IERR, IERR_P )
            PV(1) = MOD( PV(1), TWOPI )
            PV(2) = MOD( PV(2), TWOPI )
            CALL SPHXYZ( TO_XYZ, PV, PV )
            IF ( FRAMID .NE. 3 ) THEN
               CALL HCPVD( EARTH, T, PVE, IERR_E )
               IERR = MAX( IERR, IERR_E )
               PVE(1) = MOD( PVE(1), TWOPI )
               PVE(2) = MOD( PVE(2), TWOPI )
               CALL SPHXYZ( TO_XYZ, PVE, PVE )
               DO 3 J=1,6
                  PV(J) = PV(J)-PVE(J)
 3             CONTINUE
            ENDIF
         ENDIF
         IF ( P .EQ. MOON ) THEN
            CALL MOONPV( T, PVM )
            PVM(1) = MOD( PVM(1), TWOPI )
            PVM(2) = MOD( PVM(2), TWOPI )
            PVM(3) = PVM(3)/AU
            PVM(6) = PVM(6)/AU
            CALL SPHXYZ( TO_XYZ, PVM, PVM )
            DO 4 J=1,6
               PV(J) = PV(J)+PVM(J)
 4          CONTINUE
         ENDIF
C
C        LIGHT TRAVEL TIME
         TAU = SQRT( PV(1)**2+PV(2)**2+PV(3)**2 )*AU2JM
 5    CONTINUE
C
C     APPLY FK5 CORRECTIONS
      CALL SPHXYZ( TO_SPH, PV, PV )
      CALL FK5( T, PV )
C
C     APPLY LONGITUDE NUTATION CORRECTION
      CALL NUTATE( T, DPSI, DEPS, DDPSI, DDEPS )
      PV(1) = PV(1) + DPSI
      PV(4) = PV(4) + DDPSI/JM
C
C     CONVERT TO EQUATORIAL FRAME?
      IF ( FRAMID .EQ. 1 ) THEN
C
C        GET OBLIQUITY OF THE ECLIPTIC (WITH NUTATION CORRECTION)
         EPS = OBLIQ( T ) + DEPS
         MODE = -1
         CALL EQUECL( EPS, MODE, PV, PV )
      ENDIF
C
      IF ( COORID .EQ. 2 ) THEN
C
C        CONVERT TO RECTANGULAR COORDINATES
         CALL SPHXYZ( TO_XYZ, PV, PV )
      ELSE
C
C        CONVERT TO DEGREES
         PV(1) = PV(1)/D2R
         PV(2) = PV(2)/D2R
         PV(4) = PV(4)/D2R
         PV(5) = PV(5)/D2R
      ENDIF
C
      RETURN
      END