How to compute planetary positions

0. Foreword

1. Introduction

2. A few words about accuracy

3. The time scale

d = 367*y - 7 * ( y + (m+9)/12 ) / 4 + 275*m/9 + D - 730530

d = 367*y - 7 * ( y + (m+9)/12 ) / 4 - 3 * ( ( y + (m-9)/7 ) / 100 + 1 ) / 4 + 275*m/9 + D - 730515

d = d + UT/24.0 (this is a floating-point division)

4. The orbital elements

N = longitude of the ascending node i = inclination to the ecliptic (plane of the Earth's orbit) w = argument of perihelion a = semi-major axis, or mean distance from Sun e = eccentricity (0=circle, 0-1=ellipse, 1=parabola) M = mean anomaly (0 at perihelion; increases uniformly with time)

w1 = N + w = longitude of perihelion L = M + w1 = mean longitude q = a*(1-e) = perihelion distance Q = a*(1+e) = aphelion distance P = a ^ 1.5 = orbital period (years if a is in AU, astronomical units) T = Epoch_of_M - (M(deg)/360_deg) / P = time of perihelion v = true anomaly (angle between position and perihelion) E = eccentric anomaly

ecl = 23.4393 - 3.563E-7 * d

N = 0.0 i = 0.0 w = 282.9404 + 4.70935E-5 * d a = 1.000000 (AU) e = 0.016709 - 1.151E-9 * d M = 356.0470 + 0.9856002585 * d

N = 125.1228 - 0.0529538083 * d i = 5.1454 w = 318.0634 + 0.1643573223 * d a = 60.2666 (Earth radii) e = 0.054900 M = 115.3654 + 13.0649929509 * d

N = 48.3313 + 3.24587E-5 * d i = 7.0047 + 5.00E-8 * d w = 29.1241 + 1.01444E-5 * d a = 0.387098 (AU) e = 0.205635 + 5.59E-10 * d M = 168.6562 + 4.0923344368 * d

N = 76.6799 + 2.46590E-5 * d i = 3.3946 + 2.75E-8 * d w = 54.8910 + 1.38374E-5 * d a = 0.723330 (AU) e = 0.006773 - 1.302E-9 * d M = 48.0052 + 1.6021302244 * d

N = 49.5574 + 2.11081E-5 * d i = 1.8497 - 1.78E-8 * d w = 286.5016 + 2.92961E-5 * d a = 1.523688 (AU) e = 0.093405 + 2.516E-9 * d M = 18.6021 + 0.5240207766 * d

N = 100.4542 + 2.76854E-5 * d i = 1.3030 - 1.557E-7 * d w = 273.8777 + 1.64505E-5 * d a = 5.20256 (AU) e = 0.048498 + 4.469E-9 * d M = 19.8950 + 0.0830853001 * d

N = 113.6634 + 2.38980E-5 * d i = 2.4886 - 1.081E-7 * d w = 339.3939 + 2.97661E-5 * d a = 9.55475 (AU) e = 0.055546 - 9.499E-9 * d M = 316.9670 + 0.0334442282 * d

N = 74.0005 + 1.3978E-5 * d i = 0.7733 + 1.9E-8 * d w = 96.6612 + 3.0565E-5 * d a = 19.18171 - 1.55E-8 * d (AU) e = 0.047318 + 7.45E-9 * d M = 142.5905 + 0.011725806 * d

N = 131.7806 + 3.0173E-5 * d i = 1.7700 - 2.55E-7 * d w = 272.8461 - 6.027E-6 * d a = 30.05826 + 3.313E-8 * d (AU) e = 0.008606 + 2.15E-9 * d M = 260.2471 + 0.005995147 * d

5. The position of the Sun

E = M + e*(180/pi) * sin(M) * ( 1.0 + e * cos(M) )

E = M + e * sin(M) * ( 1.0 + e * cos(M) )

xv = r * cos(v) = cos(E) - e yv = r * sin(v) = sqrt(1.0 - e*e) * sin(E) v = atan2( yv, xv ) r = sqrt( xv*xv + yv*yv )

atan2( y, x ) = atan(y/x) if x positive atan2( y, x ) = atan(y/x) +- 180 degrees if x negative atan2( y, x ) = sign(y) * 90 degrees if x zero

lonsun = v + w

xs = r * cos(lonsun) ys = r * sin(lonsun)

xe = xs ye = ys * cos(ecl) ze = ys * sin(ecl)

RA = atan2( ye, xe ) Dec = atan2( ze, sqrt(xe*xe+ye*ye) )

5b. The Sidereal Time

GMST = GMST0 + UT

Ls = M + w

GMST0 = Ls + 180_degrees GMST = GMST0 + UT LST = GMST + local_longitude

GMST0 = (Ls + 180_degrees)/15 = Ls/15 + 12_hours GMST = GMST0 + UT LST = GMST + local_longitude/15

6. The position of the Moon and of the planets

M = e * sin(E) - E

E = M + e * sin(M) * ( 1.0 + e * cos(M) )

E = M + e*(180/pi) * sin(M) * ( 1.0 + e * cos(M) )

E1 = E0 - ( E0 - e*(180/pi) * sin(E0) - M ) / ( 1 - e * cos(E0) )

E1 = E0 - ( E0 - e * sin(E0) - M ) / ( 1 - e * cos(E0) )

xv = r * cos(v) = a * ( cos(E) - e ) yv = r * sin(v) = a * ( sqrt(1.0 - e*e) * sin(E) ) v = atan2( yv, xv ) r = sqrt( xv*xv + yv*yv )

7. The position in space

xh = r * ( cos(N) * cos(v+w) - sin(N) * sin(v+w) * cos(i) ) yh = r * ( sin(N) * cos(v+w) + cos(N) * sin(v+w) * cos(i) ) zh = r * ( sin(v+w) * sin(i) )

lonecl = atan2( yh, xh ) latecl = atan2( zh, sqrt(xh*xh+yh*yh) )

8. Precession

lon_corr = 3.82394E-5 * ( 365.2422 * ( Epoch - 2000.0 ) - d )

9. Perturbations of the Moon

Ms, Mm Mean Anomaly of the Sun and the Moon Nm Longitude of the Moon's node ws, wm Argument of perihelion for the Sun and the Moon Ls = Ms + ws Mean Longitude of the Sun (Ns=0) Lm = Mm + wm + Nm Mean longitude of the Moon D = Lm - Ls Mean elongation of the Moon F = Lm - Nm Argument of latitude for the Moon

-1.274 * sin(Mm - 2*D) (the Evection) +0.658 * sin(2*D) (the Variation) -0.186 * sin(Ms) (the Yearly Equation) -0.059 * sin(2*Mm - 2*D) -0.057 * sin(Mm - 2*D + Ms) +0.053 * sin(Mm + 2*D) +0.046 * sin(2*D - Ms) +0.041 * sin(Mm - Ms) -0.035 * sin(D) (the Parallactic Equation) -0.031 * sin(Mm + Ms) -0.015 * sin(2*F - 2*D) +0.011 * sin(Mm - 4*D)

-0.173 * sin(F - 2*D) -0.055 * sin(Mm - F - 2*D) -0.046 * sin(Mm + F - 2*D) +0.033 * sin(F + 2*D) +0.017 * sin(2*Mm + F)

-0.58 * cos(Mm - 2*D) -0.46 * cos(2*D)

10. Perturbations of Jupiter, Saturn and Uranus

Mj Mean anomaly of Jupiter Ms Mean anomaly of Saturn Mu Mean anomaly of Uranus (needed for Uranus only)

-0.332 * sin(2*Mj - 5*Ms - 67.6 degrees) -0.056 * sin(2*Mj - 2*Ms + 21 degrees) +0.042 * sin(3*Mj - 5*Ms + 21 degrees) -0.036 * sin(Mj - 2*Ms) +0.022 * cos(Mj - Ms) +0.023 * sin(2*Mj - 3*Ms + 52 degrees) -0.016 * sin(Mj - 5*Ms - 69 degrees)

+0.812 * sin(2*Mj - 5*Ms - 67.6 degrees) -0.229 * cos(2*Mj - 4*Ms - 2 degrees) +0.119 * sin(Mj - 2*Ms - 3 degrees) +0.046 * sin(2*Mj - 6*Ms - 69 degrees) +0.014 * sin(Mj - 3*Ms + 32 degrees)

-0.020 * cos(2*Mj - 4*Ms - 2 degrees) +0.018 * sin(2*Mj - 6*Ms - 49 degrees)

+0.040 * sin(Ms - 2*Mu + 6 degrees) +0.035 * sin(Ms - 3*Mu + 33 degrees) -0.015 * sin(Mj - Mu + 20 degrees)

11. Geocentric (Earth-centered) coordinates

xh = r * cos(lonecl) * cos(latecl) yh = r * sin(lonecl) * cos(latecl) zh = r * sin(latecl)

xs = rs * cos(lonsun) ys = rs * sin(lonsun)

xg = xh + xs yg = yh + ys zg = zh

12. Equatorial coordinates

xe = xg ye = yg * cos(ecl) - zg * sin(ecl) ze = yg * sin(ecl) + zg * cos(ecl)

RA = atan2( ye, xe ) Dec = atan2( ze, sqrt(xe*xe+ye*ye) )

rg = sqrt(xg*xg+yg*yg+zg*zg) = sqrt(xe*xe+ye*ye+ze*ze)

12b. Azimuthal coordinates

HA = LST - RA

x = cos(HA) * cos(Decl) y = sin(HA) * cos(Decl) z = sin(Decl) xhor = x * sin(lat) - z * cos(lat) yhor = y zhor = x * cos(lat) + z * sin(lat) az = atan2( yhor, xhor ) + 180_degrees alt = asin( zhor ) = atan2( zhor, sqrt(xhor*xhor+yhor*yhor) )

13. The Moon's topocentric position

mpar = asin( 1/r )

alt_topoc = alt_geoc - mpar * cos(alt_geoc)

gclat = lat, rho = 1.0

gclat = lat - 0.1924_deg * sin(2*lat) rho = 0.99833 + 0.00167 * cos(2*lat)

HA = LST - RA

g = atan( tan(gclat) / cos(HA) )

topRA = RA - mpar * rho * cos(gclat) * sin(HA) / cos(Decl) topDecl = Decl - mpar * rho * sin(gclat) * sin(g - Decl) / sin(g)

topDecl = Decl - mpar * rho * sin(-Decl) * cos(HA)

ppar = (8.794/3600)_deg / r

14. The position of Pluto

S = 50.03 + 0.033459652 * d P = 238.95 + 0.003968789 * d

lonecl = 238.9508 + 0.00400703 * d - 19.799 * sin(P) + 19.848 * cos(P) + 0.897 * sin(2*P) - 4.956 * cos(2*P) + 0.610 * sin(3*P) + 1.211 * cos(3*P) - 0.341 * sin(4*P) - 0.190 * cos(4*P) + 0.128 * sin(5*P) - 0.034 * cos(5*P) - 0.038 * sin(6*P) + 0.031 * cos(6*P) + 0.020 * sin(S-P) - 0.010 * cos(S-P) latecl = -3.9082 - 5.453 * sin(P) - 14.975 * cos(P) + 3.527 * sin(2*P) + 1.673 * cos(2*P) - 1.051 * sin(3*P) + 0.328 * cos(3*P) + 0.179 * sin(4*P) - 0.292 * cos(4*P) + 0.019 * sin(5*P) + 0.100 * cos(5*P) - 0.031 * sin(6*P) - 0.026 * cos(6*P) + 0.011 * cos(S-P) r = 40.72 + 6.68 * sin(P) + 6.90 * cos(P) - 1.18 * sin(2*P) - 0.03 * cos(2*P) + 0.15 * sin(3*P) - 0.14 * cos(3*P)

15. The elongation and physical ephemerides of the planets

d = d0 / R

Mercury 6.74" Venus 16.92" Earth 17.59" equ 17.53" pol Mars 9.36" equ 9.28" pol Jupiter 196.94" equ 185.08" pol Saturn 165.6" equ 150.8" pol Uranus 65.8" equ 62.1" pol Neptune 62.2" equ 60.9" pol

d = 1873.7" * 60 / r

elong = acos( ( s*s + R*R - r*r ) / (2*s*R) ) FV = acos( ( r*r + R*R - s*s ) / (2*r*R) )

phase = ( 1 + cos(FV) ) / 2 = hav(180_deg - FV)

hav(x) = ( 1 - cos(x) ) / 2 = sin^2 (x/2)

elong = acos( cos(slon - mlon) * cos(mlat) ) FV = 180_deg - elong

Mercury: -0.36 + 5*log10(r*R) + 0.027 * FV + 2.2E-13 * FV**6 Venus: -4.34 + 5*log10(r*R) + 0.013 * FV + 4.2E-7 * FV**3 Mars: -1.51 + 5*log10(r*R) + 0.016 * FV Jupiter: -9.25 + 5*log10(r*R) + 0.014 * FV Saturn: -9.0 + 5*log10(r*R) + 0.044 * FV + ring_magn Uranus: -7.15 + 5*log10(r*R) + 0.001 * FV Neptune: -6.90 + 5*log10(r*R) + 0.001 * FV Moon: +0.23 + 5*log10(r*R) + 0.026 * FV + 4.0E-9 * FV**4

Moon: -21.62 + 5*log10(r*R) + 0.026 * FV + 4.0E-9 * FV**4

ring_magn = -2.6 * sin(abs(B)) + 1.2 * (sin(B))**2

los = atan2( yg, xg ) las = atan2( zg, sqrt( xg*xg + yg*yg ) )

ir = 28.06_deg Nr = 169.51_deg + 3.82E-5_deg * d

B = asin( sin(las) * cos(ir) - cos(las) * sin(ir) * sin(los-Nr) )

16. Positions of asteroids

N = N_Epoch + 0.013967 * ( 2000.0 - Epoch ) + 3.82394E-5 * d

P = 365.2568984 * a**1.5 (days) = 1.00004024 * a**1.5 (years)

17. Position of comets.

M = 360.0 * (d-dT)/P (degrees)

a = q / (1.0 - e)

18. Parabolic orbits

H = (d-dT) * (k/sqrt(2)) / q**1.5

h = 1.5 * H g = sqrt( 1.0 + h*h ) s = cbrt( g + h ) - cbrt( g - h )

v = 2.0 * atan(s) r = q * ( 1.0 + s*s )

19. Near-parabolic orbits.

a = 0.75 * (d-dT) * k * sqrt( (1 + e) / (q*q*q) ) b = sqrt( 1 + a*a ) W = cbrt(b + a) - cbrt(b - a) f = (1 - e) / (1 + e) a1 = (2/3) + (2/5) * W*W a2 = (7/5) + (33/35) * W*W + (37/175) * W**4 a3 = W*W * ( (432/175) + (956/1125) * W*W + (84/1575) * W**4 ) C = W*W / (1 + W*W) g = f * C*C w = W * ( 1 + f * C * ( a1 + a2*g + a3*g*g ) ) v = 2 * atan(w) r = q * ( 1 + w*w ) / ( 1 + w*w * f )

20. Hyperbolic orbits

a = q / (1 - e)

M = (t-T) / (-a)**1.5

sinh(x) = ( exp(x) - exp(-x) ) / 2 cosh(x) = ( exp(x) + exp(-x) ) / 2

M = e * sinh(F) - F

F = M

F1 = ( M + e * ( F0 * cosh(F0) - sinh(F0) ) ) / ( e * cosh(F0) - 1 )

v = 2 * arctan( sqrt((e+1)/(e-1)) ) * tanh(F/2) r = a * ( 1 - e*e ) / ( 1 + e * cos(v) )

21. Rise and set times.

22. Validity of orbital elements.

23. Links to other sites.