aboutsummaryrefslogtreecommitdiff
path: root/math/slalib/planet.f
diff options
context:
space:
mode:
authorJoseph Hunkeler <jhunkeler@gmail.com>2015-07-08 20:46:52 -0400
committerJoseph Hunkeler <jhunkeler@gmail.com>2015-07-08 20:46:52 -0400
commitfa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 (patch)
treebdda434976bc09c864f2e4fa6f16ba1952b1e555 /math/slalib/planet.f
downloadiraf-linux-fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4.tar.gz
Initial commit
Diffstat (limited to 'math/slalib/planet.f')
-rw-r--r--math/slalib/planet.f725
1 files changed, 725 insertions, 0 deletions
diff --git a/math/slalib/planet.f b/math/slalib/planet.f
new file mode 100644
index 00000000..02590c65
--- /dev/null
+++ b/math/slalib/planet.f
@@ -0,0 +1,725 @@
+ SUBROUTINE slPLNT (DATE, NP, PV, JSTAT)
+*+
+* - - - - - - -
+* P L N 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: slPLNE
+*
+* 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 slEVP (or slEPV) or call
+* slDMON 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., 51 Franklin Street, Fifth Floor,
+* Boston, MA 02110-1301 USA
+*
+* Copyright (C) 1995 Association of Universities for Research in Astronomy Inc.
+*-
+
+ 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 slPLNE(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
generated by cgit (git 2.34.1) at 2025-09-20 10:18:15 -0400