Approximate heliocentric position and velocity of a nominated major planet: Mercury, Venus, EMB, Mars, Jupiter, Saturn, Uranus or Neptune (but not the Earth itself).
Status: support routine.
The TDB date DATE1+DATE2 is a Julian Date, apportioned in any convenient way between the two arguments. For example, JD(TDB)=2450123.7 could be expressed in any of these ways, among others:
DATE1 DATE2
2450123.7D0 0D0 (JD method)
2451545D0 -1421.3D0 (J2000 method)
2400000.5D0 50123.2D0 (MJD method)
2450123.5D0 0.2D0 (date & time method)
The JD method is the most natural and convenient to use in cases where the loss of several decimal digits of resolution is acceptable. The J2000 method is best matched to the way the argument is handled internally and will deliver the optimum resolution. The MJD method and the date & time methods are both good compromises between resolution and convenience. The limited accuracy of the present algorithm is such that any of the methods is satisfactory.
If an NP value outside the range 1-8 is supplied, an error status (J = -1) is returned and the PV vector set to zeroes.
For NP=3 the result is for the Earth-Moon Barycenter. To obtain the heliocentric position and velocity of the Earth, use instead the SOFA routine EPV00.
On successful return, the array PV contains the following:
PV(1,1) x }
PV(2,1) y } heliocentric position, au
PV(3,1) z }
PV(1,2) xdot }
PV(2,2) ydot } heliocentric velocity, au/d
PV(3,2) zdot }
The reference frame is equatorial and is with respect to the mean equator and equinox of epoch J2000.0.
The algorithm is due to J.L. Simon, P. Bretagnon, J. Chapront, M. Chapront-Touze, G. Francou and J. Laskar (Bureau des Longitudes, Paris, France). From comparisons with JPL ephemeris DE102, they quote the following maximum errors over the interval 1800-2050:
L (arcsec) B (arcsec) R (km)
Mercury 4 1 300
Venus 5 1 800
EMB 6 1 1000
Mars 17 1 7700
Jupiter 71 5 76000
Saturn 81 13 267000
Uranus 86 7 712000
Neptune 11 1 253000
Over the interval 1000-3000, they report that the accuracy is no worse than 1.5 times that over 1800-2050. Outside 1000-3000 the accuracy declines.
Comparisons of the present routine with the JPL DE200 ephemeris give the following RMS errors over the interval 1960-2025:
position (km) velocity (m/s)
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
Comparisons against DE200 over the interval 1800-2100 gave the following maximum absolute differences. (The results using DE406 were essentially the same.)
L (arcsec) B (arcsec) R (km) Rdot (m/s)
Mercury 7 1 500 0.7
Venus 7 1 1100 0.9
EMB 9 1 1300 1.0
Mars 26 1 9000 2.5
Jupiter 78 6 82000 8.2
Saturn 87 14 263000 24.6
Uranus 86 7 661000 27.4
Neptune 11 2 248000 21.4
The present SOFA re-implementation of the original Simon et al. Fortran code differs from the original in the following respects:
The date is supplied in two parts.
The result is returned only in equatorial Cartesian form; the ecliptic longitude, latitude and radius vector are not returned.
The result is in the J2000.0 equatorial frame, not ecliptic.
More is done in-line: there are fewer calls to other routines.
Different error/warning status values are used.
A different Kepler's-equation-solver is used (avoiding use of COMPLEX*16).
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 SOFA style conventions.
None of the above changes affects the result significantly.
The returned status, J, indicates the most serious condition encountered during execution of the routine. Illegal NP is considered the most serious, overriding failure to converge, which in turn takes precedence over the remote epoch warning.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=wp), | intent(in) | :: | date1 | TDB date part A (Note 1) |
||
| real(kind=wp), | intent(in) | :: | date2 | TDB date part B (Note 1) |
||
| integer, | intent(in) | :: | np | planet (1=Mercury, 2=Venus, 3=EMB ... 8=Neptune) |
||
| real(kind=wp), | intent(out), | dimension(3,2) | :: | pv | planet pos,vel (heliocentric, J2000.0, au, au/d) |
|
| integer, | intent(out) | :: | j |
|
subroutine PLAN94 ( date1, date2, np, pv, j )
implicit none
real(wp),intent(in) :: date1 !! TDB date part A (Note 1)
real(wp),intent(in) :: date2 !! TDB date part B (Note 1)
integer,intent(in) :: np !! planet (1=Mercury, 2=Venus, 3=EMB ... 8=Neptune)
real(wp),dimension(3,2),intent(out) :: pv !! planet pos,vel (heliocentric, J2000.0, au, au/d)
integer,intent(out) :: j !! status:
!! * -1 = illegal NP (outside 1-8)
!! * 0 = OK
!! * +1 = warning: date outside 1000-3000 AD
!! * +2 = warning: solution failed to converge
! Maximum number of iterations allowed to solve Kepler's equation
integer,parameter :: kmax = 10
! Days per Julian millennium
real(wp),parameter :: djm = 365250.0_wp
! Sin and cos of J2000.0 mean obliquity (IAU 1976)
real(wp),parameter :: sineps = 0.3977771559319137_wp
real(wp),parameter :: coseps = 0.9174820620691818_wp
! Gaussian constant
real(wp),parameter :: gk = 0.017202098950_wp
integer :: jstat, i, k
real(wp) :: t, da, dl, de, dp, di, dom, dmu, arga, argl, am, &
ae, dae, ae2, at, r, v, si2, xq, xp, tl, xsw, &
xcw, xm2, xf, ci2, xms, xmc, xpxq2, x, y, z
! Planetary inverse masses
real(wp),dimension(8),parameter :: amas = [ &
6023600.0_wp, 408523.5_wp, 328900.5_wp, 3098710.0_wp, &
1047.355_wp, 3498.5_wp, 22869.0_wp, 19314.0_wp ]
!
! Tables giving the mean Keplerian elements, limited to T**2 terms:
!
real(wp),dimension(3,8),parameter :: a = reshape([&
0.3870983098_wp, 0.0_wp, 0.0_wp, &
0.7233298200_wp, 0.0_wp, 0.0_wp, &
1.0000010178_wp, 0.0_wp, 0.0_wp, &
1.5236793419_wp, 3.0e-10_wp, 0.0_wp, &
5.2026032092_wp, 19132.0e-10_wp, -39.0e-10_wp, &
9.5549091915_wp, -0.0000213896_wp, 444.0e-10_wp, &
19.2184460618_wp, -3716e-10_wp, 979.0e-10_wp, &
30.1103868694_wp, -16635e-10_wp, 686.0e-10_wp ], [3,8]) !! semi-major axis (au)
!
real(wp),dimension(3,8),parameter :: dlm = reshape([&
252.25090552_wp, 5381016286.88982_wp, -1.92789_wp, &
181.97980085_wp, 2106641364.33548_wp, 0.59381_wp, &
100.46645683_wp, 1295977422.83429_wp, -2.04411_wp, &
355.43299958_wp, 689050774.93988_wp, 0.94264_wp, &
34.35151874_wp, 109256603.77991_wp, -30.60378_wp, &
50.07744430_wp, 43996098.55732_wp, 75.61614_wp, &
314.05500511_wp, 15424811.93933_wp, -1.75083_wp, &
304.34866548_wp, 7865503.20744_wp, 0.21103_wp ], [3,8]) !! mean longitude (degree and arcsecond)
!
real(wp),dimension(3,8),parameter :: e = reshape([&
0.2056317526_wp, 0.0002040653_wp, -28349e-10_wp, &
0.0067719164_wp, -0.0004776521_wp, 98127e-10_wp, &
0.0167086342_wp, -0.0004203654_wp, -0.0000126734_wp, &
0.0934006477_wp, 0.0009048438_wp, -80641e-10_wp, &
0.0484979255_wp, 0.0016322542_wp, -0.0000471366_wp, &
0.0555481426_wp, -0.0034664062_wp, -0.0000643639_wp, &
0.0463812221_wp, -0.0002729293_wp, 0.0000078913_wp, &
0.0094557470_wp, 0.0000603263_wp, 0.0_wp ], [3,8]) !! eccentricity
!
real(wp),dimension(3,8),parameter :: pi = reshape([&
77.45611904_wp, 5719.11590_wp, -4.83016_wp, &
131.56370300_wp, 175.48640_wp, -498.48184_wp, &
102.93734808_wp, 11612.35290_wp, 53.27577_wp, &
336.06023395_wp, 15980.45908_wp, -62.32800_wp, &
14.33120687_wp, 7758.75163_wp, 259.95938_wp, &
93.05723748_wp, 20395.49439_wp, 190.25952_wp, &
173.00529106_wp, 3215.56238_wp, -34.09288_wp, &
48.12027554_wp, 1050.71912_wp, 27.39717_wp ], [3,8]) !! longitude of the perihelion (degree and arcsecond)
!
real(wp),dimension(3,8),parameter :: dinc = reshape([&
7.00498625_wp, -214.25629_wp, 0.28977_wp, &
3.39466189_wp, -30.84437_wp, -11.67836_wp, &
0.0_wp, 469.97289_wp, -3.35053_wp, &
1.84972648_wp, -293.31722_wp, -8.11830_wp, &
1.30326698_wp, -71.55890_wp, 11.95297_wp, &
2.48887878_wp, 91.85195_wp, -17.66225_wp, &
0.77319689_wp, -60.72723_wp, 1.25759_wp, &
1.76995259_wp, 8.12333_wp, 0.08135_wp ], [3,8]) !! inclination (degree and arcsecond)
!
real(wp),dimension(3,8),parameter :: omega = reshape([&
48.33089304_wp, -4515.21727_wp, -31.79892_wp, &
76.67992019_wp, -10008.48154_wp, -51.32614_wp, &
174.87317577_wp, -8679.27034_wp, 15.34191_wp, &
49.55809321_wp, -10620.90088_wp, -230.57416_wp, &
100.46440702_wp, 6362.03561_wp, 326.52178_wp, &
113.66550252_wp, -9240.19942_wp, -66.23743_wp, &
74.00595701_wp, 2669.15033_wp, 145.93964_wp, &
131.78405702_wp, -221.94322_wp, -0.78728_wp ], [3,8]) !! longitude of the ascending node (degree and arcsecond)
!
! Tables for trigonometric terms to be added to the mean elements
! of the semi-major axes.
!
real(wp),dimension(9,8),parameter :: kp = reshape([&
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 ], [9,8])
!
real(wp),dimension(9,8),parameter :: ca = reshape([&
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 ], [9,8])
!
real(wp),dimension(9,8),parameter :: sa = reshape([&
-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 ], [9,8])
!
! Tables giving the trigonometric terms to be added to the mean
! elements of the mean longitudes.
!
real(wp),dimension(10,8),parameter :: kq = reshape([&
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 ], [10,8])
!
real(wp),dimension(10,8),parameter :: cl = reshape([&
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 ], [10,8])
!
real(wp),dimension(10,8),parameter :: sl = reshape([&
-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 ], [10,8])
! Validate the planet number.
if ( np<1 .or. np>8 ) then
jstat = -1
! Reset the result in case of failure.
do k=1,2
do i=1,3
pv(i,k) = 0.0_wp
end do
end do
else
! Time: Julian millennia since J2000.0.
t = ( ( date1-dj00 ) + date2 ) / djm
! OK status unless remote epoch.
if ( abs(t) <= 1.0_wp ) 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 = ( 3600.0_wp * dlm(1,np) + &
( dlm(2,np) + &
dlm(3,np) * t ) * t ) * das2r
de = e(1,np) + &
( e(2,np) + &
e(3,np) * t ) * t
dp = ANPM ( ( 3600.0_wp * pi(1,np) + &
( pi(2,np) + &
pi(3,np) * t ) * t ) * das2r )
di = ( 3600.0_wp * dinc(1,np) + &
( dinc(2,np) + &
dinc(3,np) * t ) * t ) * das2r
dom = ANPM ( ( 3600.0_wp * omega(1,np) &
+ ( omega(2,np) &
+ omega(3,np) * t ) * t ) * das2r )
! Apply the trigonometric terms.
dmu = 0.35953620_wp * t
do k=1,8
arga = kp(k,np) * dmu
argl = kq(k,np) * dmu
da = da + ( ca(k,np) * cos(arga) + &
sa(k,np) * sin(arga) ) * 1.0e-7_wp
dl = dl + ( cl(k,np) * cos(argl) + &
sl(k,np) * sin(argl) ) * 1.0e-7_wp
end do
arga = kp(9,np) * dmu
da = da + t * ( ca(9,np) * cos(arga) + &
sa(9,np) * sin(arga) ) * 1.0e-7_wp
do k=9,10
argl = kq(k,np) * dmu
dl = dl + t * ( cl(k,np) * cos(argl) + &
sl(k,np) * sin(argl) ) * 1.0e-7_wp
end do
dl = mod(dl, d2pi)
! Iterative solution of Kepler's equation to get eccentric anomaly.
am = dl - dp
ae = am + de*sin(am)
k = 0
do
dae = ( am - ae + de*sin(ae) ) / ( 1.0_wp - de*cos(ae) )
ae = ae + dae
k = k + 1
if ( k>=kmax ) jstat = 2
if ( k==kmax .or. abs(dae) <= 1.0e-12_wp ) exit
end do
! True anomaly.
ae2 = ae / 2.0_wp
at = 2.0_wp * atan2(sqrt((1.0_wp+de)/(1.0_wp-de)) * sin(ae2), &
cos(ae2))
! Distance (au) and speed (radians per day).
r = da * ( 1.0_wp - de*cos(ae) )
v = gk * sqrt( ( 1.0_wp + 1.0_wp/amas(np) ) / (da*da*da))
si2 = sin(di/2.0_wp)
xq = si2 * cos(dom)
xp = si2 * sin(dom)
tl = at + dp
xsw = sin(tl)
xcw = cos(tl)
xm2 = 2.0_wp * ( xp*xcw - xq*xsw )
xf = da / sqrt(1.0_wp - de*de)
ci2 = cos(di/2.0_wp)
xms = ( de * sin(dp) + xsw ) * xf
xmc = ( de * cos(dp) + xcw ) * xf
xpxq2 = 2.0_wp * xp * xq
! Position (J2000.0 ecliptic x,y,z in au).
x = r * ( xcw - xm2*xp )
y = r * ( xsw + xm2*xq )
z = r * ( -xm2 * ci2 )
! Rotate to equatorial.
pv(1,1) = x
pv(2,1) = y*coseps - z*sineps
pv(3,1) = y*sineps + z*coseps
! Velocity (J2000.0 ecliptic xdot,ydot,zdot in au/d).
x = v * ( ( -1.0_wp + 2.0_wp*xp*xp ) * xms + xpxq2 * xmc )
y = v * ( ( 1.0_wp - 2.0_wp*xq*xq ) * xmc - xpxq2 * xms )
z = v * ( 2.0_wp * ci2 * ( xp*xms + xq*xmc ) )
! Rotate to equatorial.
pv(1,2) = x
pv(2,2) = y*coseps - z*sineps
pv(3,2) = y*sineps + z*coseps
end if
! Return the status.
j = jstat
end subroutine PLAN94