      SUBROUTINE PLUTO( T, PV, IERR )
C
C     GIVEN DYNAMIC TIME T (JULIAN MILLENIA SINCE J2000.0) THIS
C     ROUTINE RETURNS THE TRUE HELIOCENTRIC POSITION AND VELOCITY
C     (L, B, R, dL/dT, dB/dT, dR/dT; RADIANS, RADIANS/DAY, AU,
C     AU/DAY) OF THE PLANET PLUTO, IN THE TIME INTERVAL 1885.0 TO
C     2100.0 (IERR IS RETURNED NON-ZERO IF T IS OUT OF RANGE).
C
C     B. KNAPP, 2000-08-31, 2001-04-24 (i*4)
C
C     REFERENCE: J. MEEUS, ASTRONOMICAL ALGORITHMS, 2ND ED,
C     CHAPTER 37 "PLUTO", WILLMANN-BELL 1998.
C
C
C     RCS DATA
C     
C     $Header$
C     
C     $Log$
C
C
      IMPLICIT NONE
C
C     INPUT
      REAL*8 T
C
C     OUTPUT
      REAL*8 PV(6)
      INTEGER*4 IERR
C
C     LOCAL
      REAL*8 JT, ST, PT, ALPHA, dALPHAdT, S, C
      INTEGER*4 M
C
C     CONSTANTS
      REAL*8 D2R, JM, TMIN, TMAX
      PARAMETER (D2R = 3.14159265358979324D0/180.0D0)
      PARAMETER (JM = 365250.0D0)
      PARAMETER (TMIN = (2409542.5D0-2451545.0D0)/JM)
      PARAMETER (TMAX = (2488069.5D0-2451545.0D0)/JM)
      REAL*8 J0, dJdT, S0, dSdT, P0, dPdT
      PARAMETER (J0= 34.35D0, dJdT=30349.057D0)
      PARAMETER (S0= 50.08D0, dSdT=12221.138D0)
      PARAMETER (P0=238.96D0, dPdT= 1449.600D0)
C
C     COEFFICIENTS
      INTEGER*4 N
      PARAMETER (N = 43)
      INTEGER*4 I(N)
     &   / 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1,
     &     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2,
     &     2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3/
      INTEGER*4 J(N)
     &   / 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 2, 2, 2,-1,-1,
     &     0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0,
     &     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0/
      INTEGER*4 K(N)
     &   / 1, 2, 3, 4, 5, 6,-1, 0, 1, 2, 3,-2,-1, 0, 0, 1,
     &    -3,-2,-1, 0, 1, 2, 3, 4,-3,-2,-1, 0, 1, 3,-6,-5,
     &    -4,-3,-2,-1, 0, 1, 2, 3,-2,-1, 0/
      REAL*8 L_A(N)
     &/-19799805,  897144,  611149, -341243,  129287,  -38164,   20442,
     &     -4063,   -6016,   -3956,    -667,    1276,    1152,     630,
     &      2571,     899,   -1016,   -2343,    7042,    1199,     418,
     &       120,     -60,     -82,     -36,     -40,     -14,       4,
     &         5,      -1,       2,      -4,       4,      14,     -49,
     &       163,       9,      -4,      -3,       1,      -3,       5,
     &         0/
      REAL*8 L_B(N)
     &/ 19850055,-4954829, 1211027, -189585,  -34992,   30893,   -9987,
     &     -5071,   -3336,    3039,    3572,     501,    -917,   -1277,
     &      -459,   -1449,    1043,   -1012,     788,    -338,     -67,
     &      -274,    -159,     -29,     -29,       7,      22,      13,
     &         2,       0,       0,       5,      -7,      24,     -34,
     &       -48,     -24,       1,       1,       3,      -1,      -3,
     &         0/
      REAL*8 B_A(N)
     &/ -5452852, 3527812,-1050748,  178690,   18650,  -30697,    4878,
     &       226,    2030,      69,    -247,     -57,    -122,     -49,
     &      -197,     -25,     589,    -269,     185,     315,    -130,
     &         5,       2,       2,       2,       3,       2,       1,
     &         0,       0,       0,       2,      -7,      10,      -3,
     &         6,      14,      -2,       0,       0,       0,       0,
     &         1/
      REAL*8 B_B(N)
     &/-14974862, 1672790,  327647, -292153,  100340,  -25823,   11248,
     &       -64,    -836,    -604,    -567,       1,     175,    -164,
     &       199,     217,    -248,     711,     193,     807,     -43,
     &         3,      17,       5,       3,       1,      -1,      -1,
     &        -1,       0,      -2,       2,       0,      -8,      20,
     &         5,      17,       0,       0,       0,       1,       0,
     &         0/
      REAL*8 R_A(N)
     &/66865439,-11827535, 1593179,  -18444,  -65977,   31174,   -5794,
     &      4601,   -1729,    -415,     239,      67,    1034,    -129,
     &       480,       2,   -3359,    7856,      36,    8663,    -809,
     &       263,    -126,     -35,     -19,     -15,      -4,       5,
     &         3,       6,       2,      -2,      14,     -63,     136,
     &       273,     251,     -25,       9,      -8,       2,      19,
     &        10/
      REAL*8 R_B(N)
     &/ 68951812, -332538,-1438890,  483220,  -85431,   -6032,   22161,
     &      4032,     234,     702,     723,     -67,    -451,     504,
     &      -231,    -441,     265,   -7832,   45763,    8547,    -769,
     &      -144,      32,     -16,      -4,       8,      12,       6,
     &         1,      -2,       2,      -2,      13,      13,    -236,
     &      1065,     149,      -9,      -2,       7,     -10,      35,
     &         3/
C
      SAVE I, J, K, L_A, L_B, B_A, B_B, R_A, R_B
C
C     INITIALIZE
      DO 1 M=1,6
         PV(M) = 0.0D0
 1    CONTINUE
C
C     IS T IN RANGE?
      IF ( T .LT. TMIN .OR. TMAX .LE. T ) THEN
         IERR = 1
      ENDIF
      IERR = 0
C
      JT = J0 + dJdT*T
      ST = S0 + dSdT*T
      PT = P0 + dPdT*T
C
C     NOW SUM UP THE CONTRIBUTIONS FOR LONGITUDE, LATITUDE, AND R
      DO 2 M=N,1,-1
         ALPHA = (I(M)*JT + J(M)*ST + K(M)*PT)*D2R
         dALPHAdT = (I(M)*dJdT + J(M)*dSdT + K(M)*dPdT)*D2R
         S = SIN(ALPHA)
         C = COS(ALPHA)
         PV(1) = PV(1) + (L_A(M)*S+L_B(M)*C)
         PV(2) = PV(2) + (B_A(M)*S+B_B(M)*C)
         PV(3) = PV(3) + (R_A(M)*S+R_B(M)*C)
         PV(4) = PV(4) + (L_A(M)*C-L_B(M)*S)*dALPHAdT
         PV(5) = PV(5) + (B_A(M)*C-B_B(M)*S)*dALPHAdT
         PV(6) = PV(6) + (R_A(M)*C-R_B(M)*S)*dALPHAdT
 2    CONTINUE
C
C     CONVERT TO RADIANS, RADIANS/DAY
      PV(1) = (238.958116D0 + 1449.6D0*T + PV(1)*1.0D-6)*D2R
      PV(2) = (-3.908239D0 + PV(2)*1.0D-6)*D2R
      PV(3) = (40.7241346D0 + PV(3)*1.0D-7)
      PV(4) = (1449.6D0 + PV(4)*1.0D-6)*D2R/JM
      PV(5) = (PV(5)*1.0D-6)*D2R/JM
      PV(6) = (PV(6)*1.0D-7)/JM
C
      RETURN
      END
