The position of the Moon and of the planets

First, compute the eccentric anomaly, E, from M, the mean anomaly, and e, the eccentricity. As a first approximation, do (E and M in degrees):

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

or, if E and M are in radians:

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

If e, the eccentricity, is less than about 0.05-0.06, this approximation is sufficiently accurate. If the eccentricity is larger, set E0=E and then use this iteration formula (E and M in degrees):

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

or (E and M in radians):

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

For each new iteration, replace E0 with E1. Iterate until E0 and E1 are sufficiently close together (about 0.001 degrees). For comet orbits with eccentricites close to one, a difference of less than 1E-4 or 1E-5 degrees should be required.

If this iteration formula won't converge, the eccentricity is probably too close to one. Then you should instead use the formulae for near-parabolic or parabolic orbits.

Now compute the planet's distance and true anomaly:

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 )

The position of the Sun

The position of the Sun is computed just like the position of any other planet, but since the Sun always is moving in the ecliptic, and since the eccentricity of the orbit is quite small, a few simplifications can be made. Therefore, a separate presentation for the Sun is given.

Of course, we're here really computing the position of the Earth in its orbit around the Sun, but since we're viewing the sky from an Earth-centered perspective, we'll pretend that the Sun is in orbit around the Earth instead.

First, compute the eccentric anomaly E from the mean anomaly M and from the eccentricity e (E and M in degrees):

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

or (if E and M are expressed in radians):

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

Note that the formulae for computing E are not exact; however they're accurate enough here.

Then compute the Sun's distance r and its true anomaly v from:

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 )

(note that the r computed here is later used as rs)

atan2() is a function that converts an x,y coordinate pair to the correct angle in all four quadrants. It is available as a library function in Fortran, C and C++. In other languages, one has to write one's own atan2() function. It's not that difficult:

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
Now, compute the Sun's true longitude:

lonsun = v + w

Convert lonsun,r to ecliptic rectangular geocentric coordinates xs,ys:

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

(since the Sun always is in the ecliptic plane, zs is of course zero). xs,ys is the Sun's position in a coordinate system in the plane of the ecliptic. To convert this to equatorial, rectangular, geocentric coordinates, compute:

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

Finally, compute the Sun's Right Ascension (RA) and Declination (Dec):

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

The orbital elements

The primary orbital elements are here denoted as:

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)

Related orbital elements are:

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

One Astronomical Unit (AU) is the Earth's mean distance to the Sun, or 149.6 million km. When closest to the Sun, a planet is in perihelion, and when most distant from the Sun it's in aphelion. For the Moon, an artificial satellite, or any other body orbiting the Earth, one says perigee and apogee instead, for the points in orbit least and most distant from Earth.

To describe the position in the orbit, we use three angles: Mean Anomaly, True Anomaly, and Eccentric Anomaly. They are all zero when the planet is in perihelion:
Mean Anomaly (M): This angle increases uniformly over time, by 360 degrees per orbital period. It's zero at perihelion. It's easily computed from the orbital period and the time since last perihelion.
True Anomaly (v): This is the actual angle between the planet and the perihelion, as seen from the central body (in this case the Sun). It increases non-uniformly with time, changing most rapidly at perihelion.
Eccentric Anomaly (E): This is an auxiliary angle used in Kepler's Equation, when computing the True Anomaly from the Mean Anomaly and the orbital eccentricity.
Note that for a circular orbit (eccentricity=0), these three angles are all equal to each other.


Orbital elements of the Sun:

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

Orbital elements of the Moon:

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

Orbital elements of Mercury:

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

Orbital elements of Venus:

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

Orbital elements of Mars:

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

Orbital elements of Jupiter:

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

Orbital elements of Saturn:

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

Orbital elements of Uranus:

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

Orbital elements of Neptune:

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

The time scale

The time scale in these formulae are counted in days. Hours, minutes, seconds are expressed as fractions of a day. Day 0.0 occurs at 2000 Jan 0.0 UT (or 1999 Dec 31, 0:00 UT). This "day number" d is computed as follows (y=year, m=month, D=date, UT=UT in hours+decimals):

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

Note that ALL divisions here should be INTEGER divisions. In Pascal, use "div" instead of "/", in MS-Basic, use "\" instead of "/". In Fortran, C and C++ "/" can be used if both y and m are integers. Finally, include the time of the day, by adding:

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

A few words about accuracy

The accuracy requirements are modest: a final position with an error of no more than 1-2 arc minutes (one arc minute = 1/60 degree). This accuracy is in one respect quite optimal: it is the highest accuracy one can strive for, while still being able to do many simplifications. The simplifications made here are:

1: Nutation and aberration are both ignored.
2: Planetary aberration (i.e. light travel time) is ignored.
3: The difference between Terrestial Time/Ephemeris Time (TT/ET), and Universal Time (UT) is ignored.
4: Precession is computed in a simplified way, by a simple addition to the ecliptic longitude.
5: Higher-order terms in the planetary orbital elements are ignored. This will give an additional error of up to 2 arc min in 1000 years from now. For the Moon, the error will be larger: 7 arc min 1000 years from now. This error will grow as the square of the time from the present.
6: Most planetary perturbations are ignored. Only the major perturbation terms for the Moon, Jupiter, Saturn, and Uranus, are included. If still lower accuracy is acceptable, these perturbations can be ignored as well.
7: The largest Uranus-Neptune perturbation is accounted for in the orbital elements of these planets. Therefore, the orbital elements of Uranus and Neptune are less accurace, especially in the distant past and future. The elements for these planets should therefore only be used for at most a few centuries into the past and the future.

Introduction

The text below describes how to compute the positions in the sky of the Sun, Moon and the major planets out to Neptune. The algorithm for Pluto is taken from a fourier fit to Pluto's position as computed by numerical integration at JPL. Positions of other celestial bodies as well (i.e. comets and asteroids) can also be computed, if their orbital elements are available.

These formulae may seem complicated, but I believe this is the simplest method to compute planetary positions with the fairly good accuracy of about one arc minute (=1/60 degree). Any further simplifications will yield lower accuracy, but of course that may be ok, depending on the application.