Permalink
Switch branches/tags
Nothing to show
Find file
Fetching contributors…
Cannot retrieve contributors at this time
725 lines (681 sloc) 29 KB
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 sla_EPV) 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).
*
* This revision: 19 June 2004
*
* Copyright (C) 2004 P.T.Wallace. All rights reserved.
*
* 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