Computing planetary positions
d = 367*y - 7 * ( y + (m+9)/12 ) / 4 + 275*m/9 + D - 730530Note 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)
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 anomalyOne 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.
ecl = 23.4393 - 3.563E-7 * dNow compute the orbital elements of the planet of interest. If you want the position of the Sun or the Moon, you only need to compute the orbital elements for the Sun or the Moon. If you want the position of any other planet, you must compute the orbital elements for that planet and for the Sun (of course the orbital elements for the Sun are really the orbital elements for the Earth; however it's customary to here pretend that the Sun orbits the Earth). This is necessary to be able to compute the geocentric position of the planet.
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 * dOrbital 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 * dOrbital 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 * dOrbital 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 * dOrbital 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 * dOrbital 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 * dOrbital 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 * dOrbital 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 * dOrbital 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 * dPlease note than the orbital elements of Uranus and Neptune as given here are somewhat less accurate. They include a long period perturbation between Uranus and Neptune. The period of the perturbation is about 4200 years. Therefore, these elements should not be expected to give results within the stated accuracy for more than a few centuries in the past and into the future.
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.
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( 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 zeroSee these links for some code in Basic or Pascal. Fortran and C/C++ already has atan2() as a standard library function.
lonsun = v + wConvert 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) )
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.
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 )
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) )For the Moon, this is the geocentric (Earth-centered) position in the ecliptic coordinate system. For the planets, this is the heliocentric (Sun-centered) position, also in the ecliptic coordinate system. If one wishes, one can compute the ecliptic longitude and latitude (this must be done if one wishes to correct for perturbations, or if one wants to precess the position to a standard epoch):
lonecl = atan2( yh, xh ) latecl = atan2( zh, sqrt(xh*xh+yh*yh) )As a check one can compute sqrt(xh*xh+yh*yh+zh*zh), which of course should equal r (except for small round-off errors).
lon_corr = 3.82394E-5 * ( 365.2422 * ( Epoch - 2000.0 ) - d )If one wishes the position for today's epoch (useful when computing rising/setting times and the like), no corrections need to be done.
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 MoonAdd these terms to the Moon's longitude (degrees):
-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)Add these terms to the Moon's latitude (degrees):
-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)Add these terms to the Moon's distance (Earth radii):
-0.58 * cos(Mm - 2*D) -0.46 * cos(2*D)All perturbation terms that are smaller than 0.01 degrees in longitude or latitude and smaller than 0.1 Earth radii in distance have been omitted here. A few of the largest perturbation terms even have their own names! The Evection (the largest perturbation) was discovered already by Ptolemy a few thousand years ago (the Evection was one of Ptolemy's epicycles). The Variation and the Yearly Equation were both discovered by Tycho Brahe in the 16'th century.
Mj Mean anomaly of Jupiter Ms Mean anomaly of Saturn Mu Mean anomaly of Uranus (needed for Uranus only)Perturbations for Jupiter. Add these terms to the longitude:
-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)Perturbations for Saturn. Add these terms to the longitude:
+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)For Saturn: also add these terms to the latitude:
-0.020 * cos(2*Mj - 4*Ms - 2 degrees) +0.018 * sin(2*Mj - 6*Ms - 49 degrees)Perturbations for Uranus: Add these terms to the longitude:
+0.040 * sin(Ms - 2*Mu + 6 degrees) +0.035 * sin(Ms - 3*Mu + 33 degrees) -0.015 * sin(Mj - Mu + 20 degrees)The "great Jupiter-Saturn term" is the largest perturbation for both Jupiter and Saturn. Its period is 918 years, and its amplitude is 0.332 degrees for Jupiter and 0.812 degrees for Saturn. These is also a "great Saturn-Uranus term", period 560 years, amplitude 0.035 degrees for Uranus, less than 0.01 degrees for Saturn (and therefore omitted). The other perturbations have periods between 14 and 100 years. One should also mention the "great Uranus-Neptune term", which has a period of 4220 years and an amplitude of about one degree. It is not included here, instead it is included in the orbital elements of Uranus and Neptune.
xh = r * cos(lonecl) * cos(latecl) yh = r * sin(lonecl) * cos(latecl) zh = r * sin(latecl)If we are computing the Moon's position, this is already the geocentric position, and thus we simply set xg=xh, yg=yh, zg=zh. Otherwise we must also compute the Sun's position: convert lonsun, rs (where rs is the r computed here) to xs, ys:
xs = rs * cos(lonsun) ys = rs * sin(lonsun)(Of course, any correction for precession should be added to lonecl and lonsun before converting to xh,yh,zh and xs,ys).
xg = xh + xs yg = yh + ys zg = zhWe now have the planet's geocentric (Earth centered) position in rectangular, ecliptic coordinates.
xe = xg ye = yg * cos(ecl) - zg * sin(ecl) ze = yg * sin(ecl) + zg * cos(ecl)Finally, compute the planet's Right Ascension (RA) and Declination (Dec):
RA = atan2( ye, xe ) Dec = atan2( ze, sqrt(xe*xe+ye*ye) )Compute the geocentric distance:
rg = sqrt(xg*xg+yg*yg+zg*zg) = sqrt(xe*xe+ye*ye+ze*ze)
mpar = asin( 1/r )where r is the Moon's distance in Earth radii. It's simplest to apply the correction in horizontal coordinates (azimuth and altitude): within our accuracy aim of 1-2 arc minutes, no correction need to be applied to the azimuth. One need only apply a correction to the altitude above the horizon:
alt_topoc = alt_geoc - mpar * cos(alt_geoc)Sometimes one need to correct for topocentric position directly in equatorial coordinates though, e.g. if one wants to draw on a star map how the Moon passes in front of the Pleiades, as seen from some specific location. Then we need to know the Moon's geocentric Right Ascension and Declination (RA, Decl), the Local Sidereal Time (LST), and our latitude (lat).
gclat = lat, rho = 1.0However, if we do wish to account for the flattening of the Earth, we instead compute:
gclat = lat - 0.1924_deg * sin(2*lat) rho = 0.99883 + 0.00167 * cos(2*lat)Next we compute the Moon's geocentric Hour Angle (HA):
HA = LST - RAwhere LST is our Local Sidereal Time, computed as:
LST = GMST0 + UT + LON/15where UT is the Universal Time in hours, LON is the observer's geographical longitude (east longitude positive, west negative). GMST0 is the Greenwich Mean Sidereal Time at 0h UT, and is most easily computed by adding, or subtracting, 180 degrees to Ls, the Sun's mean longitude, and then converting from degrees to hours by duvuding by 15:
GMST0 = ( Ls + 180_deg ) / 15We also need an auxiliary angle, g:
g = atan( tan(gclat) / cos(HA) )Now we're ready to convert the geocentric Right Ascention and Declination (RA, Decl) to their topocentric values (topRA, topDecl):
topRA = RA - mpar * rho * cos(gclat) * sin(HA) / cos(Decl) topDecl = Decl - mpar * rho * sin(gclat) * sin(g - Decl) / sin(g)(Note that if decl is exactly 90 deg, cos(Decl) becomes zero and we get a division by zero when computing topRA, but that formula breaks down very close to the celestial poles anyway. Also if gclat is precisely zero, g becomes zero too, and we get a division by zero when computing topDecl. In that case, replace the formula for topDecl with
topDecl = Decl - mpar * rho * sin(-Decl) * cos(HA)which is valid for gclat equal to zero; it can also be used for gclat extremely close to zero).
ppar = (8.794/3600)_deg / rwhere r is the distance of the planet from the Earth, in astronomical units.
S = 50.03 + 0.033459652 * d P = 238.95 + 0.003968789 * dNext compute the heliocentric ecliptic longitude and latitude (degrees), and distance (a.u.):
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)Now we know the heliocentric distance and ecliptic longitude/latitude for Pluto. To convert to geocentric coordinates, do as for the other planets.
d = d0 / RR is the planet's geocentric distance in astronomical units, and d is the planet's apparent diameter at a distance of 1 astronomical unit. d0 is of course different for each planet. The values below are given in seconds of arc. Some planets have different equatorial and polar diameters:
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" polThe Sun's apparent diameter at 1 astronomical unit is 1919.26". The Moon's apparent diameter is:
d = 1873.7" * 60 / rwhere r is the Moon's distance in Earth radii.
elong = acos( ( s*s + R*R - r*r ) / (2*s*R) ) FV = acos( ( r*r + R*R - s*s ) / (2*r*R) )When we know the phase angle, we can easily compute the phase:
phase = ( 1 + cos(FV) ) / 2 = hav(180_deg - FV)hav(FV) is the "haversine" of the phase angle. The "haversine" (or "half versine") is an old and now obsolete trigonometric function; it's defined as:
hav(x) = ( 1 - cos(x) ) / 2 = sin^2 (x/2)As usual we must use a different procedure for the Moon. Since the Moon is so close to the Earth, the procedure above would introduce too big errors. Instead we use the Moon's ecliptic longitude and latitude, mlon and mlat, and the Sun's ecliptic longitude, mlon, to compute first the elongation, then the phase angle, of the Moon:
elong = acos( cos(slon - mlon) * cos(mlat) ) FV = 180_deg - elongFinally we'll compute the magnitude (or brightness) of the planets. Here we need to use a formula that's different for each planet. FV is the phase angle (in degrees), r is the heliocentric and R the geocentric distance (both in AU):
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** is the power operator, thus FV**6 is the phase angle (in degrees) raised to the sixth power. If FV is 150 degrees then FV**6 becomes ca 1.14E+13, which is a quite large number.
Moon: -21.62 + 5*log10(r*R) + 0.026 * FV + 4.0E-9 * FV**4Saturn needs special treatment due to its rings: when Saturn's rings are "open" then Saturn will appear much brighter than when we view Saturn's rings edgewise. We'll compute ring_mang like this:
ring_magn = -2.6 * sin(abs(B)) + 1.2 * (sin(B))**2Here B is the tilt of Saturn's rings which we also need to compute. Then we start with Saturn's geocentric ecliptic longitude and latitude (los, las) which we've already computed. We also need the tilt of the rings to the ecliptic, ir, and the "ascending node" of the plane of the rings, Nr:
ir = 28.06_deg Nr = 169.51_deg + 3.82E-5_deg * dHere d is our "day number" which we've used so many times before. Now let's compute the tilt of the rings:
B = asin( sin(las) * cos(ir) - cos(las) * sin(ir) * sin(los-Nr) )This concludes our computation of the magnitudes of the planets.
N = N_Epoch + 0.013967 * ( 2000.0 - Epoch ) + 3.82394E-5 * dwhere Epoch is expressed as a year with fractions, e.g. 1950.0 or 2000.0
P = 365.2568984 * a**1.5 (days) = 1.00004024 * a**1.5 (years)** is the power-of operator. a**1.5 is the same as sqrt(a*a*a).
M = 360.0 * (d-dT)/P (degrees)where P is given in days, and d-dT of course is the time since last perihelion, also in days.
a = q / (1.0 - e)Then proceed as with an asteroid (section 16).
H = (d-dT) / (82.21168627 * q**1.5)where q**1.5 is the same as sqrt(q*q*q). Also compute:
h = 1.5 * H g = sqrt( 1.0 + h*h ) s = cbrt( g + h ) - cbrt( g - h )cbrt() is the cube root function: cbrt(x) = x**(1.0/3.0). The formulae has been devised so that both g+h and g-h always are positive. Therefore one can here safely compute cbrt(x) as exp(log(x)/3.0) . In general, cbrt(-x) = -cbrt(x) and of course cbrt(0) = 0.
v = 2.0 * atan(s) r = q * ( 1.0 + s*s )When we know the true anomaly and the heliocentric distance, we can proceed by computing the position in space (section 7).
a = 0.75 * (d-dT) * k * sqrt( (1 + e) / (q*q*q) ) b = sqrt( 1 + a*a ) W = cbrt(b + a) - cbrt(b - a) c = 1 + 1/(W*W) f = (1 - e) / (1 + e) g = f / (c*c) 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 ) w = W * ( 1 + g * c * ( a1 + a2*g + a3*g*g ) ) v = 2 * atan(w) r = q * ( 1 + w*w ) / ( 1 + w*w * f )This algorithm yields the true anomaly, v, and the heliocentric distance, r, for a nearly-parabolic orbit. Note that this algorithm will fail very far from the perihelion; however the accuracy is sufficient for all comets closer than Pluto.