      SUBROUTINE sla_PLANET (DATE, NP, PV, JSTAT)
*+
*     - - - - - - -
*      P L A N E T
*     - - - - - - -
*
*  Approximate heliocentric position and velocity of a specified
*  major planet.
*
*  Given:
*     DATE      d      Modified Julian Date (JD - 2400000.5)
*     NP        i      planet (1=Mercury, 2=Venus, 3=EMB ... 9=Pluto)
*
*  Returned:
*     PV        d(6)   heliocentric x,y,z,xdot,ydot,zdot, J2000
*                                           equatorial triad (AU,AU/s)
*     JSTAT     i      status: +1 = warning: date out of range
*                               0 = OK
*                              -1 = illegal NP (outside 1-9)
*                              -2 = solution didn't converge
*
*  Called:  sla_PLANEL
*
*  Notes
*
*  1  The epoch, DATE, is in the TDB timescale and is a Modified
*     Julian Date (JD-2400000.5).
*
*  2  The reference frame is equatorial and is with respect to the
*     mean equinox and ecliptic of epoch J2000.
*
*  3  If an NP value outside the range 1-9 is supplied, an error
*     status (JSTAT = -1) is returned and the PV vector set to zeroes.
*
*  4  The algorithm for obtaining the mean elements of the planets
*     from Mercury to Neptune is due to J.L. Simon, P. Bretagnon,
*     J. Chapront, M. Chapront-Touze, G. Francou and J. Laskar
*     (Bureau des Longitudes, Paris).  The (completely different)
*     algorithm for calculating the ecliptic coordinates of Pluto
*     is by Meeus.
*
*  5  Comparisons of the present routine with the JPL DE200 ephemeris
*     give the following RMS errors over the interval 1960-2025:
*
*                      position (km)     speed (metre/sec)
*
*        Mercury            334               0.437
*        Venus             1060               0.855
*        EMB               2010               0.815
*        Mars              7690               1.98
*        Jupiter          71700               7.70
*        Saturn          199000              19.4
*        Uranus          564000              16.4
*        Neptune         158000              14.4
*        Pluto            36400               0.137
*
*     From comparisons with DE102, Simon et al quote the following
*     longitude accuracies over the interval 1800-2200:
*
*        Mercury                 4"
*        Venus                   5"
*        EMB                     6"
*        Mars                   17"
*        Jupiter                71"
*        Saturn                 81"
*        Uranus                 86"
*        Neptune                11"
*
*     In the case of Pluto, Meeus quotes an accuracy of 0.6 arcsec
*     in longitude and 0.2 arcsec in latitude for the period
*     1885-2099.
*
*     For all except Pluto, over the period 1000-3000 the accuracy
*     is better than 1.5 times that over 1800-2200.  Outside the
*     period 1000-3000 the accuracy declines.  For Pluto the
*     accuracy declines rapidly outside the period 1885-2099.
*     Outside these ranges (1885-2099 for Pluto, 1000-3000 for
*     the rest) a "date out of range" warning status (JSTAT=+1)
*     is returned.
*
*  6  The algorithms for (i) Mercury through Neptune and (ii) Pluto
*     are completely independent.  In the Mercury through Neptune
*     case, the present SLALIB implementation differs from the
*     original Simon et al Fortran code in the following respects.
*
*     *  The date is supplied as a Modified Julian Date rather
*        than a Julian Date (MJD = JD - 2400000.5).
*
*     *  The result is returned only in equatorial Cartesian form;
*        the ecliptic longitude, latitude and radius vector are not
*        returned.
*
*     *  The velocity is in AU per second, not AU per day.
*
*     *  Different error/warning status values are used.
*
*     *  Kepler's equation is not solved inline.
*
*     *  Polynomials in T are nested to minimize rounding errors.
*
*     *  Explicit double-precision constants are used to avoid
*        mixed-mode expressions.
*
*     *  There are other, cosmetic, changes to comply with
*        Starlink/SLALIB style guidelines.
*
*     None of the above changes affects the result significantly.
*
*  7  For NP=3 the result is for the Earth-Moon Barycentre.  To
*     obtain the heliocentric position and velocity of the Earth,
*     either use the SLALIB routine sla_EVP or call sla_DMOON and
*     subtract 0.012150581 times the geocentric Moon vector from
*     the EMB vector produced by the present routine.  (The Moon
*     vector should be precessed to J2000 first, but this can
*     be omitted for modern epochs without introducing significant
*     inaccuracy.)
*
*  References:  Simon et al., Astron. Astrophys. 282, 663 (1994).
*               Meeus, Astronomical Algorithms, Willmann-Bell (1991).
*
*  P.T.Wallace   Starlink   27 May 1997
*
*  Copyright (C) 1997 Rutherford Appleton Laboratory
*
*  License:
*    This program is free software; you can redistribute it and/or modify
*    it under the terms of the GNU General Public License as published by
*    the Free Software Foundation; either version 2 of the License, or
*    (at your option) any later version.
*
*    This program is distributed in the hope that it will be useful,
*    but WITHOUT ANY WARRANTY; without even the implied warranty of
*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
*    GNU General Public License for more details.
*
*    You should have received a copy of the GNU General Public License
*    along with this program (see SLA_CONDITIONS); if not, write to the 
*    Free Software Foundation, Inc., 59 Temple Place, Suite 330, 
*    Boston, MA  02111-1307  USA
*
*-

      IMPLICIT NONE

      DOUBLE PRECISION DATE
      INTEGER NP
      DOUBLE PRECISION PV(6)
      INTEGER JSTAT

*  2Pi, deg to radians, arcsec to radians
      DOUBLE PRECISION D2PI,D2R,AS2R
      PARAMETER (D2PI=6.283185307179586476925286766559D0,
     :           D2R=0.017453292519943295769236907684886D0,
     :           AS2R=4.848136811095359935899141023579D-6)

*  Gaussian gravitational constant (exact)
      DOUBLE PRECISION GCON
      PARAMETER (GCON=0.01720209895D0)

*  Seconds per Julian century
      DOUBLE PRECISION SPC
      PARAMETER (SPC=36525D0*86400D0)

*  Sin and cos of J2000 mean obliquity (IAU 1976)
      DOUBLE PRECISION SE,CE
      PARAMETER (SE=0.3977771559319137D0,
     :           CE=0.9174820620691818D0)

      INTEGER I,J,IJSP(3,43)
      DOUBLE PRECISION AMAS(8),A(3,8),DLM(3,8),E(3,8),
     :                 PI(3,8),DINC(3,8),OMEGA(3,8),
     :                 DKP(9,8),CA(9,8),SA(9,8),
     :                 DKQ(10,8),CLO(10,8),SLO(10,8),
     :                 T,DA,DE,DPE,DI,DO,DMU,ARGA,ARGL,DM,
     :                 AB(2,3,43),DJ0,DS0,DP0,DL0,DLD0,DB0,DR0,
     :                 DJ,DS,DP,DJD,DSD,DPD,WLBR(3),WLBRD(3),
     :                 WJ,WS,WP,AL,ALD,SAL,CAL,
     :                 AC,BC,DL,DLD,DB,DBD,DR,DRD,
     :                 SL,CL,SB,CB,SLCB,CLCB,X,Y,Z,XD,YD,ZD

*  -----------------------
*  Mercury through Neptune
*  -----------------------

*  Planetary inverse masses
      DATA AMAS / 6023600D0,408523.5D0,328900.5D0,3098710D0,
     :            1047.355D0,3498.5D0,22869D0,19314D0 /

*
*  Tables giving the mean Keplerian elements, limited to T**2 terms:
*
*         A       semi-major axis (AU)
*         DLM     mean longitude (degree and arcsecond)
*         E       eccentricity
*         PI      longitude of the perihelion (degree and arcsecond)
*         DINC    inclination (degree and arcsecond)
*         OMEGA   longitude of the ascending node (degree and arcsecond)
*
      DATA A /
     :  0.3870983098D0,             0D0,      0D0,
     :  0.7233298200D0,             0D0,      0D0,
     :  1.0000010178D0,             0D0,      0D0,
     :  1.5236793419D0,           3D-10,      0D0,
     :  5.2026032092D0,       19132D-10,  -39D-10,
     :  9.5549091915D0, -0.0000213896D0,  444D-10,
     : 19.2184460618D0,       -3716D-10,  979D-10,
     : 30.1103868694D0,      -16635D-10,  686D-10 /
*
      DATA DLM /
     : 252.25090552D0, 5381016286.88982D0,  -1.92789D0,
     : 181.97980085D0, 2106641364.33548D0,   0.59381D0,
     : 100.46645683D0, 1295977422.83429D0,  -2.04411D0,
     : 355.43299958D0,  689050774.93988D0,   0.94264D0,
     :  34.35151874D0,  109256603.77991D0, -30.60378D0,
     :  50.07744430D0,   43996098.55732D0,  75.61614D0,
     : 314.05500511D0,   15424811.93933D0,  -1.75083D0,
     : 304.34866548D0,    7865503.20744D0,   0.21103D0/
*
      DATA E /
     : 0.2056317526D0,  0.0002040653D0,      -28349D-10,
     : 0.0067719164D0, -0.0004776521D0,       98127D-10,
     : 0.0167086342D0, -0.0004203654D0, -0.0000126734D0,
     : 0.0934006477D0,  0.0009048438D0,      -80641D-10,
     : 0.0484979255D0,  0.0016322542D0, -0.0000471366D0,
     : 0.0555481426D0, -0.0034664062D0, -0.0000643639D0,
     : 0.0463812221D0, -0.0002729293D0,  0.0000078913D0,
     : 0.0094557470D0,  0.0000603263D0,            0D0 /
*
      DATA PI /
     :  77.45611904D0,  5719.11590D0,   -4.83016D0,
     : 131.56370300D0,   175.48640D0, -498.48184D0,
     : 102.93734808D0, 11612.35290D0,   53.27577D0,
     : 336.06023395D0, 15980.45908D0,  -62.32800D0,
     :  14.33120687D0,  7758.75163D0,  259.95938D0,
     :  93.05723748D0, 20395.49439D0,  190.25952D0,
     : 173.00529106D0,  3215.56238D0,  -34.09288D0,
     :  48.12027554D0,  1050.71912D0,   27.39717D0 /
*
      DATA DINC /
     : 7.00498625D0, -214.25629D0,   0.28977D0,
     : 3.39466189D0,  -30.84437D0, -11.67836D0,
     :          0D0,  469.97289D0,  -3.35053D0,
     : 1.84972648D0, -293.31722D0,  -8.11830D0,
     : 1.30326698D0,  -71.55890D0,  11.95297D0,
     : 2.48887878D0,   91.85195D0, -17.66225D0,
     : 0.77319689D0,  -60.72723D0,   1.25759D0,
     : 1.76995259D0,    8.12333D0,   0.08135D0 /
*
      DATA OMEGA /
     :  48.33089304D0,  -4515.21727D0,  -31.79892D0,
     :  76.67992019D0, -10008.48154D0,  -51.32614D0,
     : 174.87317577D0,  -8679.27034D0,   15.34191D0,
     :  49.55809321D0, -10620.90088D0, -230.57416D0,
     : 100.46440702D0,   6362.03561D0,  326.52178D0,
     : 113.66550252D0,  -9240.19942D0,  -66.23743D0,
     :  74.00595701D0,   2669.15033D0,  145.93964D0,
     : 131.78405702D0,   -221.94322D0,   -0.78728D0 /
*
*  Tables for trigonometric terms to be added to the mean elements
*  of the semi-major axes.
*
      DATA DKP /
     : 69613, 75645, 88306, 59899, 15746, 71087, 142173,  3086,    0,
     : 21863, 32794, 26934, 10931, 26250, 43725,  53867, 28939,    0,
     : 16002, 21863, 32004, 10931, 14529, 16368,  15318, 32794,    0,
     : 6345,   7818, 15636,  7077,  8184, 14163,   1107,  4872,    0,
     : 1760,   1454,  1167,   880,   287,  2640,     19,  2047, 1454,
     :  574,      0,   880,   287,    19,  1760,   1167,   306,  574,
     :  204,      0,   177,  1265,     4,   385,    200,   208,  204,
     :    0,    102,   106,     4,    98,  1367,    487,   204,    0 /
*
      DATA CA /
     :      4,    -13,    11,    -9,    -9,    -3,    -1,     4,    0,
     :   -156,     59,   -42,     6,    19,   -20,   -10,   -12,    0,
     :     64,   -152,    62,    -8,    32,   -41,    19,   -11,    0,
     :    124,    621,  -145,   208,    54,   -57,    30,    15,    0,
     : -23437,  -2634,  6601,  6259, -1507, -1821,  2620, -2115,-1489,
     :  62911,-119919, 79336, 17814,-24241, 12068,  8306, -4893, 8902,
     : 389061,-262125,-44088,  8387,-22976, -2093,  -615, -9720, 6633,
     :-412235,-157046,-31430, 37817, -9740,   -13, -7449,  9644,    0 /
*
      DATA SA /
     :     -29,    -1,     9,     6,    -6,     5,     4,     0,    0,
     :     -48,  -125,   -26,   -37,    18,   -13,   -20,    -2,    0,
     :    -150,   -46,    68,    54,    14,    24,   -28,    22,    0,
     :    -621,   532,  -694,   -20,   192,   -94,    71,   -73,    0,
     :  -14614,-19828, -5869,  1881, -4372, -2255,   782,   930,  913,
     :  139737,     0, 24667, 51123, -5102,  7429, -4095, -1976,-9566,
     : -138081,     0, 37205,-49039,-41901,-33872,-27037,-12474,18797,
     :       0, 28492,133236, 69654, 52322,-49577,-26430, -3593,    0 /
*
*  Tables giving the trigonometric terms to be added to the mean
*  elements of the mean longitudes.
*
      DATA DKQ /
     :  3086, 15746, 69613, 59899, 75645, 88306, 12661, 2658,  0,   0,
     : 21863, 32794, 10931,    73,  4387, 26934,  1473, 2157,  0,   0,
     :    10, 16002, 21863, 10931,  1473, 32004,  4387,   73,  0,   0,
     :    10,  6345,  7818,  1107, 15636,  7077,  8184,  532, 10,   0,
     :    19,  1760,  1454,   287,  1167,   880,   574, 2640, 19,1454,
     :    19,   574,   287,   306,  1760,    12,    31,   38, 19, 574,
     :     4,   204,   177,     8,    31,   200,  1265,  102,  4, 204,
     :     4,   102,   106,     8,    98,  1367,   487,  204,  4, 102 /
*
      DATA CLO /
     :     21,   -95, -157,   41,   -5,   42,   23,   30,     0,    0,
     :   -160,  -313, -235,   60,  -74,  -76,  -27,   34,     0,    0,
     :   -325,  -322,  -79,  232,  -52,   97,   55,  -41,     0,    0,
     :   2268,  -979,  802,  602, -668,  -33,  345,  201,   -55,    0,
     :   7610, -4997,-7689,-5841,-2617, 1115, -748, -607,  6074,  354,
     : -18549, 30125,20012, -730,  824,   23, 1289, -352,-14767,-2062,
     :-135245,-14594, 4197,-4030,-5630,-2898, 2540, -306,  2939, 1986,
     :  89948,  2103, 8963, 2695, 3682, 1648,  866, -154, -1963, -283 /
*
      DATA SLO /
     :   -342,   136,  -23,   62,   66,  -52,  -33,   17,     0,    0,
     :    524,  -149,  -35,  117,  151,  122,  -71,  -62,     0,    0,
     :   -105,  -137,  258,   35, -116,  -88, -112,  -80,     0,    0,
     :    854,  -205, -936, -240,  140, -341,  -97, -232,   536,    0,
     : -56980,  8016, 1012, 1448,-3024,-3710,  318,  503,  3767,  577,
     : 138606,-13478,-4964, 1441,-1319,-1482,  427, 1236, -9167,-1918,
     :  71234,-41116, 5334,-4935,-1848,   66,  434,-1748,  3780, -701,
     : -47645, 11647, 2166, 3194,  679,    0, -244, -419, -2531,   48 /

*  -----
*  Pluto
*  -----

*
*  Coefficients for fundamental arguments:  mean longitudes
*  (degrees) and mean rate of change of longitude (degrees per
*  Julian century) for Jupiter, Saturn and Pluto
*
      DATA DJ0, DJD / 34.35D0, 3034.9057D0 /
      DATA DS0, DSD / 50.08D0, 1222.1138D0 /
      DATA DP0, DPD / 238.96D0, 144.9600D0 /

*  Coefficients for latitude, longitude, radius vector
      DATA DL0,DLD0 / 238.956785D0, 144.96D0 /
      DATA DB0 / -3.908202D0 /
      DATA DR0 / 40.7247248D0 /

*
*  Coefficients for periodic terms (Meeus's Table 36.A)
*
*  The coefficients for term n in the series are:
*
*    IJSP(1,n)     J
*    IJSP(2,n)     S
*    IJSP(3,n)     P
*    AB(1,1,n)     longitude sine (degrees)
*    AB(2,1,n)     longitude cosine (degrees)
*    AB(1,2,n)     latitude sine (degrees)
*    AB(2,2,n)     latitude cosine (degrees)
*    AB(1,3,n)     radius vector sine (AU)
*    AB(2,3,n)     radius vector cosine (AU)
*
      DATA (IJSP(I, 1),I=1,3),((AB(J,I, 1),J=1,2),I=1,3) /
     :                             0,  0,  1,
     :            -19798886D-6,  19848454D-6,
     :             -5453098D-6, -14974876D-6,
     :             66867334D-7,  68955876D-7 /
      DATA (IJSP(I, 2),I=1,3),((AB(J,I, 2),J=1,2),I=1,3) /
     :                             0,  0,  2,
     :               897499D-6,  -4955707D-6,
     :              3527363D-6,   1672673D-6,
     :            -11826086D-7,   -333765D-7 /
      DATA (IJSP(I, 3),I=1,3),((AB(J,I, 3),J=1,2),I=1,3) /
     :                             0,  0,  3,
     :               610820D-6,   1210521D-6,
     :             -1050939D-6,    327763D-6,
     :              1593657D-7,  -1439953D-7 /
      DATA (IJSP(I, 4),I=1,3),((AB(J,I, 4),J=1,2),I=1,3) /
     :                             0,  0,  4,
     :              -341639D-6,   -189719D-6,
     :               178691D-6,   -291925D-6,
     :               -18948D-7,    482443D-7 /
      DATA (IJSP(I, 5),I=1,3),((AB(J,I, 5),J=1,2),I=1,3) /
     :                             0,  0,  5,
     :               129027D-6,    -34863D-6,
     :                18763D-6,    100448D-6,
     :               -66634D-7,    -85576D-7 /
      DATA (IJSP(I, 6),I=1,3),((AB(J,I, 6),J=1,2),I=1,3) /
     :                             0,  0,  6,
     :               -38215D-6,     31061D-6,
     :               -30594D-6,    -25838D-6,
     :                30841D-7,     -5765D-7 /
      DATA (IJSP(I, 7),I=1,3),((AB(J,I, 7),J=1,2),I=1,3) /
     :                             0,  1, -1,
     :                20349D-6,     -9886D-6,
     :                 4965D-6,     11263D-6,
     :                -6140D-7,     22254D-7 /
      DATA (IJSP(I, 8),I=1,3),((AB(J,I, 8),J=1,2),I=1,3) /
     :                             0,  1,  0,
     :                -4045D-6,     -4904D-6,
     :                  310D-6,      -132D-6,
     :                 4434D-7,      4443D-7 /
      DATA (IJSP(I, 9),I=1,3),((AB(J,I, 9),J=1,2),I=1,3) /
     :                             0,  1,  1,
     :                -5885D-6,     -3238D-6,
     :                 2036D-6,      -947D-6,
     :                -1518D-7,       641D-7 /
      DATA (IJSP(I,10),I=1,3),((AB(J,I,10),J=1,2),I=1,3) /
     :                             0,  1,  2,
     :                -3812D-6,      3011D-6,
     :                   -2D-6,      -674D-6,
     :                   -5D-7,       792D-7 /
      DATA (IJSP(I,11),I=1,3),((AB(J,I,11),J=1,2),I=1,3) /
     :                             0,  1,  3,
     :                 -601D-6,      3468D-6,
     :                 -329D-6,      -563D-6,
     :                  518D-7,       518D-7 /
      DATA (IJSP(I,12),I=1,3),((AB(J,I,12),J=1,2),I=1,3) /
     :                             0,  2, -2,
     :                 1237D-6,       463D-6,
     :                  -64D-6,        39D-6,
     :                  -13D-7,      -221D-7 /
      DATA (IJSP(I,13),I=1,3),((AB(J,I,13),J=1,2),I=1,3) /
     :                             0,  2, -1,
     :                 1086D-6,      -911D-6,
     :                  -94D-6,       210D-6,
     :                  837D-7,      -494D-7 /
      DATA (IJSP(I,14),I=1,3),((AB(J,I,14),J=1,2),I=1,3) /
     :                             0,  2,  0,
     :                  595D-6,     -1229D-6,
     :                   -8D-6,      -160D-6,
     :                 -281D-7,       616D-7 /
      DATA (IJSP(I,15),I=1,3),((AB(J,I,15),J=1,2),I=1,3) /
     :                             1, -1,  0,
     :                 2484D-6,      -485D-6,
     :                 -177D-6,       259D-6,
     :                  260D-7,      -395D-7 /
      DATA (IJSP(I,16),I=1,3),((AB(J,I,16),J=1,2),I=1,3) /
     :                             1, -1,  1,
     :                  839D-6,     -1414D-6,
     :                   17D-6,       234D-6,
     :                 -191D-7,      -396D-7 /
      DATA (IJSP(I,17),I=1,3),((AB(J,I,17),J=1,2),I=1,3) /
     :                             1,  0, -3,
     :                 -964D-6,      1059D-6,
     :                  582D-6,      -285D-6,
     :                -3218D-7,       370D-7 /
      DATA (IJSP(I,18),I=1,3),((AB(J,I,18),J=1,2),I=1,3) /
     :                             1,  0, -2,
     :                -2303D-6,     -1038D-6,
     :                 -298D-6,       692D-6,
     :                 8019D-7,     -7869D-7 /
      DATA (IJSP(I,19),I=1,3),((AB(J,I,19),J=1,2),I=1,3) /
     :                             1,  0, -1,
     :                 7049D-6,       747D-6,
     :                  157D-6,       201D-6,
     :                  105D-7,     45637D-7 /
      DATA (IJSP(I,20),I=1,3),((AB(J,I,20),J=1,2),I=1,3) /
     :                             1,  0,  0,
     :                 1179D-6,      -358D-6,
     :                  304D-6,       825D-6,
     :                 8623D-7,      8444D-7 /
      DATA (IJSP(I,21),I=1,3),((AB(J,I,21),J=1,2),I=1,3) /
     :                             1,  0,  1,
     :                  393D-6,       -63D-6,
     :                 -124D-6,       -29D-6,
     :                 -896D-7,      -801D-7 /
      DATA (IJSP(I,22),I=1,3),((AB(J,I,22),J=1,2),I=1,3) /
     :                             1,  0,  2,
     :                  111D-6,      -268D-6,
     :                   15D-6,         8D-6,
     :                  208D-7,      -122D-7 /
      DATA (IJSP(I,23),I=1,3),((AB(J,I,23),J=1,2),I=1,3) /
     :                             1,  0,  3,
     :                  -52D-6,      -154D-6,
     :                    7D-6,        15D-6,
     :                 -133D-7,        65D-7 /
      DATA (IJSP(I,24),I=1,3),((AB(J,I,24),J=1,2),I=1,3) /
     :                             1,  0,  4,
     :                  -78D-6,       -30D-6,
     :                    2D-6,         2D-6,
     :                  -16D-7,         1D-7 /
      DATA (IJSP(I,25),I=1,3),((AB(J,I,25),J=1,2),I=1,3) /
     :                             1,  1, -3,
     :                  -34D-6,       -26D-6,
     :                    4D-6,         2D-6,
     :                  -22D-7,         7D-7 /
      DATA (IJSP(I,26),I=1,3),((AB(J,I,26),J=1,2),I=1,3) /
     :                             1,  1, -2,
     :                  -43D-6,         1D-6,
     :                    3D-6,         0D-6,
     :                   -8D-7,        16D-7 /
      DATA (IJSP(I,27),I=1,3),((AB(J,I,27),J=1,2),I=1,3) /
     :                             1,  1, -1,
     :                  -15D-6,        21D-6,
     :                    1D-6,        -1D-6,
     :                    2D-7,         9D-7 /
      DATA (IJSP(I,28),I=1,3),((AB(J,I,28),J=1,2),I=1,3) /
     :                             1,  1,  0,
     :                   -1D-6,        15D-6,
     :                    0D-6,        -2D-6,
     :                   12D-7,         5D-7 /
      DATA (IJSP(I,29),I=1,3),((AB(J,I,29),J=1,2),I=1,3) /
     :                             1,  1,  1,
     :                    4D-6,         7D-6,
     :                    1D-6,         0D-6,
     :                    1D-7,        -3D-7 /
      DATA (IJSP(I,30),I=1,3),((AB(J,I,30),J=1,2),I=1,3) /
     :                             1,  1,  3,
     :                    1D-6,         5D-6,
     :                    1D-6,        -1D-6,
     :                    1D-7,         0D-7 /
      DATA (IJSP(I,31),I=1,3),((AB(J,I,31),J=1,2),I=1,3) /
     :                             2,  0, -6,
     :                    8D-6,         3D-6,
     :                   -2D-6,        -3D-6,
     :                    9D-7,         5D-7 /
      DATA (IJSP(I,32),I=1,3),((AB(J,I,32),J=1,2),I=1,3) /
     :                             2,  0, -5,
     :                   -3D-6,         6D-6,
     :                    1D-6,         2D-6,
     :                    2D-7,        -1D-7 /
      DATA (IJSP(I,33),I=1,3),((AB(J,I,33),J=1,2),I=1,3) /
     :                             2,  0, -4,
     :                    6D-6,       -13D-6,
     :                   -8D-6,         2D-6,
     :                   14D-7,        10D-7 /
      DATA (IJSP(I,34),I=1,3),((AB(J,I,34),J=1,2),I=1,3) /
     :                             2,  0, -3,
     :                   10D-6,        22D-6,
     :                   10D-6,        -7D-6,
     :                  -65D-7,        12D-7 /
      DATA (IJSP(I,35),I=1,3),((AB(J,I,35),J=1,2),I=1,3) /
     :                             2,  0, -2,
     :                  -57D-6,       -32D-6,
     :                    0D-6,        21D-6,
     :                  126D-7,      -233D-7 /
      DATA (IJSP(I,36),I=1,3),((AB(J,I,36),J=1,2),I=1,3) /
     :                             2,  0, -1,
     :                  157D-6,       -46D-6,
     :                    8D-6,         5D-6,
     :                  270D-7,      1068D-7 /
      DATA (IJSP(I,37),I=1,3),((AB(J,I,37),J=1,2),I=1,3) /
     :                             2,  0,  0,
     :                   12D-6,       -18D-6,
     :                   13D-6,        16D-6,
     :                  254D-7,       155D-7 /
      DATA (IJSP(I,38),I=1,3),((AB(J,I,38),J=1,2),I=1,3) /
     :                             2,  0,  1,
     :                   -4D-6,         8D-6,
     :                   -2D-6,        -3D-6,
     :                  -26D-7,        -2D-7 /
      DATA (IJSP(I,39),I=1,3),((AB(J,I,39),J=1,2),I=1,3) /
     :                             2,  0,  2,
     :                   -5D-6,         0D-6,
     :                    0D-6,         0D-6,
     :                    7D-7,         0D-7 /
      DATA (IJSP(I,40),I=1,3),((AB(J,I,40),J=1,2),I=1,3) /
     :                             2,  0,  3,
     :                    3D-6,         4D-6,
     :                    0D-6,         1D-6,
     :                  -11D-7,         4D-7 /
      DATA (IJSP(I,41),I=1,3),((AB(J,I,41),J=1,2),I=1,3) /
     :                             3,  0, -2,
     :                   -1D-6,        -1D-6,
     :                    0D-6,         1D-6,
     :                    4D-7,       -14D-7 /
      DATA (IJSP(I,42),I=1,3),((AB(J,I,42),J=1,2),I=1,3) /
     :                             3,  0, -1,
     :                    6D-6,        -3D-6,
     :                    0D-6,         0D-6,
     :                   18D-7,        35D-7 /
      DATA (IJSP(I,43),I=1,3),((AB(J,I,43),J=1,2),I=1,3) /
     :                             3,  0,  0,
     :                   -1D-6,        -2D-6,
     :                    0D-6,         1D-6,
     :                   13D-7,         3D-7 /


*  Validate the planet number.
      IF (NP.LT.1.OR.NP.GT.9) THEN
         JSTAT=-1
         DO I=1,6
            PV(I)=0D0
         END DO
      ELSE

*     Separate algorithms for Pluto and the rest.
         IF (NP.NE.9) THEN

*        -----------------------
*        Mercury through Neptune
*        -----------------------

*        Time: Julian millennia since J2000.
            T=(DATE-51544.5D0)/365250D0

*        OK status unless remote epoch.
            IF (ABS(T).LE.1D0) THEN
               JSTAT=0
            ELSE
               JSTAT=1
            END IF

*        Compute the mean elements.
            DA=A(1,NP)+(A(2,NP)+A(3,NP)*T)*T
            DL=(3600D0*DLM(1,NP)+(DLM(2,NP)+DLM(3,NP)*T)*T)*AS2R
            DE=E(1,NP)+(E(2,NP)+E(3,NP)*T)*T
            DPE=MOD((3600D0*PI(1,NP)+(PI(2,NP)+PI(3,NP)*T)*T)*AS2R,D2PI)
            DI=(3600D0*DINC(1,NP)+(DINC(2,NP)+DINC(3,NP)*T)*T)*AS2R
            DO=MOD((3600D0*OMEGA(1,NP)
     :                        +(OMEGA(2,NP)+OMEGA(3,NP)*T)*T)*AS2R,D2PI)

*        Apply the trigonometric terms.
            DMU=0.35953620D0*T
            DO J=1,8
               ARGA=DKP(J,NP)*DMU
               ARGL=DKQ(J,NP)*DMU
               DA=DA+(CA(J,NP)*COS(ARGA)+SA(J,NP)*SIN(ARGA))*1D-7
               DL=DL+(CLO(J,NP)*COS(ARGL)+SLO(J,NP)*SIN(ARGL))*1D-7
            END DO
            ARGA=DKP(9,NP)*DMU
            DA=DA+T*(CA(9,NP)*COS(ARGA)+SA(9,NP)*SIN(ARGA))*1D-7
            DO J=9,10
               ARGL=DKQ(J,NP)*DMU
               DL=DL+T*(CLO(J,NP)*COS(ARGL)+SLO(J,NP)*SIN(ARGL))*1D-7
            END DO
            DL=MOD(DL,D2PI)

*        Daily motion.
            DM=GCON*SQRT((1D0+1D0/AMAS(NP))/(DA*DA*DA))

*        Make the prediction.
            CALL sla_PLANEL(DATE,1,DATE,DI,DO,DPE,DA,DE,DL,DM,PV,J)
            IF (J.LT.0) JSTAT=-2

         ELSE

*        -----
*        Pluto
*        -----

*        Time: Julian centuries since J2000.
            T=(DATE-51544.5D0)/36525D0

*        OK status unless remote epoch.
            IF (T.GE.-1.15D0.AND.T.LE.1D0) THEN
               JSTAT=0
            ELSE
               JSTAT=1
            END IF

*        Fundamental arguments (radians).
            DJ=(DJ0+DJD*T)*D2R
            DS=(DS0+DSD*T)*D2R
            DP=(DP0+DPD*T)*D2R

*        Initialize coefficients and derivatives.
            DO I=1,3
               WLBR(I)=0D0
               WLBRD(I)=0D0
            END DO

*        Term by term through Meeus Table 36.A.
            DO J=1,43

*           Argument and derivative (radians, radians per century).
               WJ=DBLE(IJSP(1,J))
               WS=DBLE(IJSP(2,J))
               WP=DBLE(IJSP(3,J))
               AL=WJ*DJ+WS*DS+WP*DP
               ALD=(WJ*DJD+WS*DSD+WP*DPD)*D2R

*           Functions of argument.
               SAL=SIN(AL)
               CAL=COS(AL)

*           Periodic terms in longitude, latitude, radius vector.
               DO I=1,3

*              A and B coefficients (deg, AU).
                  AC=AB(1,I,J)
                  BC=AB(2,I,J)

*              Periodic terms (deg, AU, deg/Jc, AU/Jc).
                  WLBR(I)=WLBR(I)+AC*SAL+BC*CAL
                  WLBRD(I)=WLBRD(I)+(AC*CAL-BC*SAL)*ALD
               END DO
            END DO

*        Heliocentric longitude and derivative (radians, radians/sec).
            DL=(DL0+DLD0*T+WLBR(1))*D2R
            DLD=(DLD0+WLBRD(1))*D2R/SPC

*        Heliocentric latitude and derivative (radians, radians/sec).
            DB=(DB0+WLBR(2))*D2R
            DBD=WLBRD(2)*D2R/SPC

*        Heliocentric radius vector and derivative (AU, AU/sec).
            DR=DR0+WLBR(3)
            DRD=WLBRD(3)/SPC

*        Functions of latitude, longitude, radius vector.
            SL=SIN(DL)
            CL=COS(DL)
            SB=SIN(DB)
            CB=COS(DB)
            SLCB=SL*CB
            CLCB=CL*CB

*        Heliocentric vector and derivative, J2000 ecliptic and equinox.
            X=DR*CLCB
            Y=DR*SLCB
            Z=DR*SB
            XD=DRD*CLCB-DR*(CL*SB*DBD+SLCB*DLD)
            YD=DRD*SLCB+DR*(-SL*SB*DBD+CLCB*DLD)
            ZD=DRD*SB+DR*CB*DBD

*        Transform to J2000 equator and equinox.
            PV(1)=X
            PV(2)=Y*CE-Z*SE
            PV(3)=Y*SE+Z*CE
            PV(4)=XD
            PV(5)=YD*CE-ZD*SE
            PV(6)=YD*SE+ZD*CE
         END IF
      END IF

      END
