Computing planetary positionsHow to compute planetary positionsBy Paul Schlyter, Stockholm, Swedenemail: pausch@stjarnhimlen.se or WWW: http://stjarnhimlen.se/Break out of a frame0. Foreword1. Introduction2. A few words about accuracy3. The time scale4. The orbital elements5. The position of the Sun6. The position of the Moon and of the planets7. The position in space8. Precession9. Perturbations of the Moon10. Perturbations of Jupiter, Saturn and Uranus11. Geocentric (Earth-centered) coordinates12. Equatorial coordinates13. The Moon's topocentric position14. The position of Pluto15. The elongation and physical ephemerides of the planets16. Positions of asteroids17. Position of comets18. Parabolic orbits19. Near-parabolic orbits20. Rise and set times21. Validity of orbital elements22. Links to other sitesTutorial with numerical test cases0. ForewordBelow is a description of how to compute the positions for the Sun andMoon and the major planets, as well as for comets and minor planets,from a set of orbital elements.The algorithms have been simplified as much as possible while stillkeeping a fairly good accuracy. The accuracy of the computedpositions is a fraction of an arc minute for the Sun and the innerplanets, about one arc minute for the outer planets, and 1-2 arcminutes for the Moon. If we limit our accuracy demands to thislevel, one can simplify further by e.g. ignoring the differencebetween mean, true and apparent positions.The positions computed below are for the 'equinox of the day', whichis suitable for computing rise/set times, but not for plotting theposition on a star map drawn for a fixed epoch. In the latter case,correction for precession must be applied, which is most simplyperformed as a rotation along the ecliptic.These algortihms were developed by myself back in 1979, based on apreprint from T. van Flandern's and K. Pulkkinen's paper "Lowprecision formulae for planetary positions", published in theAstrophysical Journal Supplement Series, 1980. It's basically asimplification of these algorithms, while keeping a reasonableaccuracy. They were first implemented on a HP-41C programmablepocket calculator, in 1979, and ran in less than 2 KBytes of RAM!Nowadays considerable more accurate algorithms are available ofcourse, as well as more powerful computers. Nevertheless I'veretained these algorithms as what I believe is the simplest way tocompute solar/lunar positions with an accuracy of 1-2 arc minutes.1. IntroductionThe text below describes how to compute the positions in the sky ofthe Sun, Moon and the major planets out to Neptune. The algorithmfor Pluto is taken from a fourier fit to Pluto's position as computedby numerical integration at JPL. Positions of other celestial bodiesas well (i.e. comets and asteroids) can also be computed, if theirorbital elements are available.These formulae may seem complicated, but I believe this is thesimplest method to compute planetary positions with the fairly goodaccuracy of about one arc minute (=1/60 degree). Any furthersimplifications will yield lower accuracy, but of course that may beok, depending on the application.2. A few words about accuracyThe accuracy requirements are modest: a final position with an errorof no more than 1-2 arc minutes (one arc minute = 1/60 degree). Thisaccuracy is in one respect quite optimal: it is the highest accuracyone can strive for, while still being able to do manysimplifications. 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.3. The time scaleThe time scale in these formulae are counted in days. Hours,minutes, seconds are expressed as fractions of a day. Day 0.0 occursat 2000 Jan 0.0 UT (or 1999 Dec 31, 0:00 UT). This "day number" d iscomputed as follows (y=year, m=month, D=date, UT=UT inhours+decimals): 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 "/". InFortran, 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)4. The orbital elementsThe 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 anomalyOne Astronomical Unit (AU) is the Earth's mean distance to theSun, or 149.6 million km. When closest to the Sun, a planet is inperihelion, and when most distant from the Sun it's inaphelion. For the Moon, an artificial satellite, or any otherbody orbiting the Earth, one says perigee and apogeeinstead, for the points in orbit least and most distant from Earth.To describe the position in the orbit, we use three angles: MeanAnomaly, True Anomaly, and Eccentric Anomaly. They are all zero when theplanet is in perihelion:Mean Anomaly (M): This angle increases uniformly over time, by360 degrees per orbital period. It's zero at perihelion. It'seasily computed from the orbital period and the time since lastperihelion.True Anomaly (v): This is the actual angle between the planetand the perihelion, as seen from the central body (in this case theSun). It increases non-uniformly with time, changing most rapidly atperihelion.Eccentric Anomaly (E): This is an auxiliary angle used inKepler's Equation, when computing the True Anomaly from the MeanAnomaly and the orbital eccentricity.Note that for a circular orbit (eccentricity=0), these three anglesare all equal to each other.Another quantity we will need is ecl, the obliquity of theecliptic, i.e. the "tilt" of the Earth's axis of rotation(currently ca 23.4 degrees and slowly decreasing). First, compute the"d" of the moment of interest (section 3). Then,compute the obliquity of the ecliptic: ecl = 23.4393 - 3.563E-7 * dNow compute the orbital elements of the planet of interest. If youwant the position of the Sun or the Moon, you only need to computethe orbital elements for the Sun or the Moon. If you want theposition of any other planet, you must compute the orbital elementsfor that planet and for the Sun (of course the orbitalelements for the Sun are really the orbital elements for the Earth;however it's customary to here pretend that the Sun orbits theEarth). This is necessary to be able to compute the geocentricposition of the planet.Please note that a, the semi-major axis, is given in Earth radii forthe Moon, but in Astronomical Units for the Sun and all the planets.When computing M (and, for the Moon, when computing N and w as well),one will quite often get a result that is larger than 360 degrees, ornegative (all angles are here computed in degrees). If negative,add 360 degrees until positive. If larger than 360 degrees, subtract360 degrees until the value is less than 360 degrees. Note that, inmost programming languages, one must then multiply these angles withpi/180 to convert them to radians, before taking the sine or cosineof them.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 * 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 givenhere are somewhat less accurate. They include a long periodperturbation between Uranus and Neptune. The period of theperturbation is about 4200 years. Therefore, these elements shouldnot be expected to give results within the stated accuracy for morethan a few centuries in the past and into the future.5. The position of the SunThe position of the Sun is computed just like the position of anyother planet, but since the Sun always is moving in the ecliptic, andsince the eccentricity of the orbit is quite small, a fewsimplifications can be made. Therefore, a separate presentation forthe Sun is given.Of course, we're here really computing the position of the Earth inits orbit around the Sun, but since we're viewing the sky from anEarth-centered perspective, we'll pretend that the Sun is in orbitaround the Earth instead.First, compute the eccentric anomaly E from the mean anomaly Mand 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'reaccurate 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 thecorrect angle in all four quadrants. It is available as a libraryfunction in Fortran, C and C++. In other languages, one has to writeone'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 zeroSee these links for some code in Basic orPascal. Fortran and C/C++ already hasatan2() as a standard library function. Now, compute the Sun's true longitude: lonsun = v + wConvert lonsun,r to ecliptic rectangular geocentric coordinatesxs,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 ofthe ecliptic. To convert this to equatorial, rectangular, geocentriccoordinates, 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) )6. The position of the Moon and of the planetsFirst, compute the eccentric anomaly, E, from M, the mean anomaly,and e, the eccentricity. As a first approximation, do (E and M indegrees): 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, thisapproximation is sufficiently accurate. If the eccentricity islarger, set E0=E and then use this iteration formula (E and M indegrees): 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 E1are sufficiently close together (about 0.001 degrees). For cometorbits with eccentricites close to one, a difference of less than1E-4 or 1E-5 degrees should be required.If this iteration formula won't converge, the eccentricity isprobably too close to one. Then you should instead use the formulaefor near-parabolic or parabolicorbits.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 )7. The position in spaceCompute the planet's position in 3-dimensional 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) )For the Moon, this is the geocentric (Earth-centered) position in theecliptic coordinate system. For the planets, this is theheliocentric (Sun-centered) position, also in the ecliptic coordinatesystem. If one wishes, one can compute the ecliptic longitude andlatitude (this must be done if one wishes to correct forperturbations, or if one wants to precess the position to a standardepoch): 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 courseshould equal r (except for small round-off errors).8. PrecessionIf one wishes to compute the planet's position for some standardepoch, such as 1950.0 or 2000.0 (e.g. to be able to plot the positionon a star atlas), one must add the correction below to lonecl. If aplanet's and not the Moon's position is computed, one must also addthe same correction to lonsun, the Sun's longitude. The desiredEpoch is expressed as the year, possibly with a fraction. lon_corr = 3.82394E-5 * ( 365.2422 * ( Epoch - 2000.0 ) - d )If one wishes the position for today's epoch (useful when computingrising/setting times and the like), no corrections need to be done.9. Perturbations of the MoonIf the position of the Moon is computed, and one wishes a betteraccuracy than about 2 degrees, the most important perturbations hasto be taken into account. If one wishes 2 arc minute accuracy, allthe following terms should be accounted for. If less accuracy isneeded, some of the smaller terms can be omitted.First compute: 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 Add 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 inlongitude or latitude and smaller than 0.1 Earth radii in distancehave been omitted here. A few of the largest perturbation terms evenhave their own names! The Evection (the largest perturbation) wasdiscovered already by Ptolemy a few thousand years ago (the Evectionwas one of Ptolemy's epicycles). The Variation and the YearlyEquation were both discovered by Tycho Brahe in the 16'thcentury.The computations can be simplified by omitting the smallerperturbation terms. The error introduced by this seldom exceeds thesum of the amplitudes of the 4-5 largest omitted terms. If one onlycomputes the three largest perturbation terms in longitude and thelargest term in latitude, the error in longitude will rarley exceed0.25 degrees, and in latitude 0.15 degrees.10. Perturbations of Jupiter, Saturn and UranusThe only planets having perturbations larger than 0.01 degrees areJupiter, Saturn and Uranus. First compute: 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 bothJupiter and Saturn. Its period is 918 years, and its amplitude is0.332 degrees for Jupiter and 0.812 degrees for Saturn. These isalso a "great Saturn-Uranus term", period 560 years, amplitude 0.035degrees for Uranus, less than 0.01 degrees for Saturn (and thereforeomitted). The other perturbations have periods between 14 and 100years. One should also mention the "great Uranus-Neptune term",which has a period of 4220 years and an amplitude of about onedegree. It is not included here, instead it is included in theorbital elements of Uranus and Neptune.For Mercury, Venus and Mars we can ignore all perturbations. ForNeptune the only significant perturbation is already included in theorbital elements, as mentioned above, and therefore no furtherperturbation terms need to be accounted for.11. Geocentric (Earth-centered) coordinatesNow we have computed the heliocentric (Sun-centered) coordinate ofthe planet, and we have included the most important perturbations.We want to compute the geocentric (Earth-centerd) position. Weshould convert the perturbed lonecl, latecl, r to (perturbed) xh, yh,zh: 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 geocentricposition, and thus we simply set xg=xh, yg=yh, zg=zh. Otherwise wemust also compute the Sun's position: convert lonsun, rs (where rs isthe r computed here) to xs, ys: xs = rs * cos(lonsun) ys = rs * sin(lonsun)(Of course, any correction for precession should be added to lonecland lonsun before converting to xh,yh,zh and xs,ys).Now convert from heliocentric to geocentric position: xg = xh + xs yg = yh + ys zg = zhWe now have the planet's geocentric (Earth centered) position inrectangular, ecliptic coordinates.12. Equatorial coordinatesLet's convert our rectangular, ecliptic coordinates to rectangular,equatorial coordinates: simply rotate the y-z-plane by ecl, the angleof the obliquity of the ecliptic: 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)13. The Moon's topocentric positionThe Moon's position, as computed earlier, is geocentric, i.e. as seenby an imaginary observer at the center of the Earth. Real observersdwell on the surface of the Earth, though, and they will see adifferent position - the topocentric position. This position candiffer by more than one degree from the geocentric position. Tocompute the topocentric positions, we must add a correction to thegeocentric position.Let's start by computing the Moon's parallax, i.e. the apparentsize of the (equatorial) radius of the Earth, as seen from the Moon: mpar = asin( 1/r )where r is the Moon's distance in Earth radii. It's simplest to applythe correction in horizontal coordinates (azimuth and altitude):within our accuracy aim of 1-2 arc minutes, no correction need to beapplied to the azimuth. One need only apply a correction to thealtitude above the horizon: alt_topoc = alt_geoc - mpar * cos(alt_geoc)Sometimes one need to correct for topocentric position directly inequatorial coordinates though, e.g. if one wants to draw on a starmap how the Moon passes in front of the Pleiades, as seen from somespecific location. Then we need to know the Moon's geocentricRight Ascension and Declination (RA, Decl), the Local SiderealTime (LST), and our latitude (lat).Our astronomical latitude (lat) must first be converted to ageocentric latitude (gclat), and distance from the center of the Earth(rho) in Earth equatorial radii. If we only want an approximatetopocentric position, it's simplest to pretend that the Earth isa perfect sphere, and simply set: 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.99833 + 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'sgeographical longitude (east longitude positive, west negative).GMST0 is the Greenwich Mean Sidereal Time at 0h UT, and is mosteasily computed by adding, or subtracting, 180 degrees to Ls, theSun's mean longitude, and then converting from degrees to hours byduvuding 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 andDeclination (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 adivision by zero when computing topRA, but that formula breaks downvery close to the celestial poles anyway. Also if gclat is preciselyzero, g becomes zero too, and we get a division by zero when computingtopDecl. 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 gclatextremely close to zero).This correction to topocentric position can also be applied to theSun and the planets. But since they're much farther away, thecorrection becomes much smaller. It's largest for Venus at inferiorconjunction, when Venus' parallax is somewhat larger than 32 arcseconds. Within our aim of obtaining a final accuracy of 1-2 arcminutes, it might barely be justified to correct to topocentricposition when Venus is close to inferior conjunction, and perhapsalso when Mars is at a favourable opposition. But in all other casesthis correction can safely be ignored within our accuracy aim. Weonly need to worry about the Moon in this case.If you want to compute topocentric coordinates for the planets anyway,you do it the same way as for the Moon, with one exception: theparallax of the planet (ppar) is computed from this formula: ppar = (8.794/3600)_deg / rwhere r is the distance of the planet from the Earth, in astronomicalunits.14. The position of PlutoNo analytical theory has ever been constructed for the planet Pluto.Our most accurate representation of the motion of this planet is fromnumerical integrations. Yet, a "curve fit" may be performed to thesenumerical integrations, and the result will be the formulae below,valid from ca 1800 to ca 2100.Compute d, our day number, as usual (section 3). Thencompute these angles: 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/latitudefor Pluto. To convert to geocentric coordinates, do as for the otherplanets.15. The elongation and physical ephemerides of the planetsWhen we finally have completed our computation of the heliocentricand geocentric coordinates of the planets, it could also beinteresting to know what the planet will look like. How large willit appear? What's its phase and magnitude (brightness)? Thesecomputations are much simpler than the computations of thepositions.Let's start by computing the apparent diameter of the planet: d = d0 / RR is the planet's geocentric distance in astronomical units, andd is the planet's apparent diameter at a distance of 1 astronomicalunit. d0 is of course different for each planet. The valuesbelow are given in seconds of arc. Some planets have differentequatorial 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". TheMoon's apparent diameter is: d = 1873.7" * 60 / rwhere r is the Moon's distance in Earth radii.Two other quantities we'd like to know are the phase angle and theelongation.The phase angle tells us the phase: if it's zero the planet appears"full", if it's 90 degrees it appears "half", and if it's 180 degreesit appears "new". Only the Moon and the inferior planets (i.e.Mercury and Venus) can have phase angles exceeding about 50 degrees.The elongation is the apparent angular distance of the planet fromthe Sun. If the elongation is smaller than ca 20 degrees, the planetis hard to observe, and if it's smaller than ca 10 degrees it'susually not possible to observe the planet.To compute phase angle and elongation we need to know the planet'sheliocentric distance, r, its geocentric distance, R, and thedistance to the Sun, s. Now we can compute the phase angle, FV,and the elongation, elong: 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 theMoon is so close to the Earth, the procedure above would introducetoo big errors. Instead we use the Moon's ecliptic longitude andlatitude, mlon and mlat, and the Sun's ecliptic longitude, mlon, tocompute 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. FVis the phase angle (in degrees), r is the heliocentric and R thegeocentric 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 becomesca 1.14E+13, which is a quite large number.For the Moon, we also need the heliocentric distance, r, andgeocentric distance, R, in AU (astronomical units). Here r can be setequal to the Sun's geocentric distance in AU. The Moon's geocentricdistance, R, previously computed i Earth radii, must be convertedto AU's - we do this by multiplying by sin(17.59"/2) = 1/23450. Orwe could modify the magnitude formula for the Moon so it usesr in AU's and R in Earth radii: Moon: -21.62 + 5*log10(r*R) + 0.026 * FV + 4.0E-9 * FV**4Saturn needs special treatment due to its rings: when Saturn'srings are "open" then Saturn will appear much brighter than whenwe view Saturn's rings edgewise. We'll compute ring_mang likethis: 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 andlatitude (los, las) which we've already computed. We also need thetilt of the rings to the ecliptic, ir, and the "ascending node" ofthe 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.16. Positions of asteroidsFor asteroids, the orbital elements are often given as: N,i,w,a,e,M,where N,i,w are valid for a specific epoch (nowadays usually 2000.0).In our simplified computational scheme, the only significant changeswith the epoch occurs in N. To convert N_Epoch to the N (today'sepoch) we want to use, simply add a correction for precession: 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.0Most often M, the mean anomaly, is given for another day than the daywe want to compute the asteroid's position for. If the daily motion,n, is given, simply add n * (time difference in days) to M. If n isnot given, but the period P (in days) is given, then n = 360.0/P. IfP is not given, it can be computed from: 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).When all orbital elements has been computed, proceed as with theother planets (section 6).17. Position of comets.For comets having elliptical orbits, M is usually not given. InsteadT, the time of perihelion, is given. At perihelion M is zero. Tocompute M for any other moment, first compute the "day number" d of T(section 3), let's call this dT. Then compute the"day number" d of the moment for which you want to compute aposition, let's call this d. Then M, the mean anomaly, is computedlike: M = 360.0 * (d-dT)/P (degrees)where P is given in days, and d-dT of course is the time since lastperihelion, also in days.Also, a, the semi-major axis, is usually not given. Instead q, theperihelion distance, is given. a can easily be computed from q ande: a = q / (1.0 - e)Then proceed as with an asteroid (section 16).18. Parabolic orbitsIf the comet has a parabolic orbit, a different method has to beused. Then the orbital period of the comet is infinite, and M (themean anomaly) is always zero. The eccentricity, e, is always exactly1. Since the semi-major axis, a, is infinite, we must insteaddirectly use the perihelion distance, q. To compute a parabolicorbit, we proceed like this:Compute the "day number", d, for T, the moment of perihelion, callthis dT. Compute d for the moment we want to compute a position,call it d (section 3). The constant k is the Gaussiangravitational constant: k = 0.01720209895 exactly!Then compute: H = (d-dT) * (k/sqrt(2)) / q**1.5where 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). Theformulae has been devised so that both g+h and g-h always arepositive. Therefore one can here safely compute cbrt(x) asexp(log(x)/3.0) . In general, cbrt(-x) = -cbrt(x) and of coursecbrt(0) = 0.Instead of trying to compute some eccentric anomaly, we compute thetrue anomaly and the heliocentric distance directly: v = 2.0 * atan(s) r = q * ( 1.0 + s*s )When we know the true anomaly and the heliocentric distance, we canproceed by computing the position in space (section 7).19. Near-parabolic orbits. The most common case for a newly discovered comet is that the orbitisn't an exact parabola, but very nearly so. It's eccentricity isslightly below, or slightly above, one. The algorithm presented herecan be used for eccentricities between about 0.98 and 1.02. If theeccentricity is smaller than 0.98 the elliptic algorithm (Kepler'sequation/etc) should be used instead. No known comet has aneccentricity exceeding 1.02.As for the purely parabolic orbit, we start by computing the timesince perihelion in days, d-dT, and the perihelion distance, q. Wealso need to know the eccentricity, e. The constant k is theGaussian gravitational constant: k = 0.01720209895 exactly!Then we can proceed as: 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 )This algorithm yields the true anomaly, v, and the heliocentricdistance, r, for a nearly-parabolic orbit. Note that this algorithmwill fail very far from the perihelion; however the accuracy issufficient for all comets closer than Pluto.20. Rise and set times.(this subject has received a document of itsown)21. Validity of orbital elements.Due to perturbations from mainly the giant planets, like Jupiter andSaturn, the orbital elements of celestial bodies are constantlychanging. The orbital elements for the Sun, Moon and the majorplanets, as given here, are valid for a long time period. However,orbital elemets given for a comet or an asteroid are valid only for alimited time. How long they are valid is hard to say generally. Itdepends on many factors, such as the accuracy you need, and themagnitude of the perturbations the comet or asteroid is subjected tofrom, say, Jupiter. A comet might travel in roughly the same orbitseveral orbital periods, experiencing only slight perturbations, butsuddenly it might pass very close to Jupiter and get its orbitchanged drastically. To compute this in a reliable way is quitecomplicated and completely out of scope for this description. As arule of thumb, though, one can assume that an asteriod, if one usesthe orbital elements for a specific epoch, one or a few revolutionsaway from that moment will have an error in its computed position ofat least one or a few arc minutes, and possibly more. The errorswill accumulate with time.22. Links to other sites.Astronomical Calculations by Keith Burnett: http://www.xylem.f2s.com/kepler/Free BASIC programs can be found at ftp://seds.lpl.arizona.edu/pub/software/pc/general/ in: ast.exe (needs GWBASIC!) and duff2ed.exe (Pete Duffett-Smiths programs)Books from Willmann-Bell about Math and Celestial Mechanics: http://www.willbell.com/math/index.htmJohn Walker's freeware program Home Planet + other stuff:http://www.fourmilab.ch/Elwood Downey's Xephem and Ephem programs, with C source code:http://www.clearskyinstitute.com/xephem/.Ephem can also be found at ftp://seds.lpl.arizona.edu/pub/software/pc/general/ as ephem421.zipSteven Moshier: Astronomy and numerical software source codes: http://www.moshier.net/Dan Bruton's astronomical software links: http://www.physics.sfasu.edu/astro/software.htmlMel Bartel's software (much ATM stuff): http://www.efn.org/~mbartels/tm/software.htmlAlmanac data from USNO: http://aa.usno.navy.mil/data/Asteroid orbital elements from Lowell Observatory: http://asteroid.lowell.edu/SAC downloads: http://www.saguaroastro.org/content/downloads.htmEarth Satellite software from AMSAT: http://www.amsat.org/amsat/ftpsoft.htmlBureau des Longitudes: (site is reconstructured so these links do not work right now) http://www.bdl.fr/ VSOP87: ftp://www.bdl.fr/pub/ephem/planets/vsop87DE200, DE450 and DE403/404 appear to no longer be available:DE200 at NRAO: ftp://sadira.gb.nrao.edu/incoming/SSEphem/DE403/404 at JPL: ftp://navigator.jpl.nasa.gov:/ephem/export/Some catalogues at CDS, Strasbourg, France - high accuracy orbital theories: Overview: http://cdsweb.u-strasbg.fr/htbin/myqcat3?VI/ Precession & mean orbital elements: http://cdsweb.u-strasbg.fr/htbin/myqcat3?VI/66/ ELP2000-82 (orbital theory of Moon): http://cdsweb.u-strasbg.fr/htbin/myqcat3?VI/79/ VSOP87 (orbital theories of planets): http://cdsweb.u-strasbg.fr/htbin/myqcat3?VI/81/ Astronomical Data Center http://adc.gsfc.nasa.gov/adc.html had formerly lots of catalogs, but that service has now been discontinued. Similar services are available at: CDS (France) - ADAC (Japan) - CAD (Russia) Formerly, these catalogs were available at the ADC - maybe I'll find equivalents elsewhere later: Asteroid orbital elements 1998: ftp://adc.gsfc.nasa.gov/pub/adc/archives/catalogs/1/1245/ JPL ephemeris DE118/LE62: ftp://adc.gsfc.nasa.gov/pub/adc/archives/catalogs/1/1093A/ JPL ephemeris DE200/LE200: ftp://adc.gsfc.nasa.gov/pub/adc/archives/catalogs/1/1094A/ USNO ZZCAT: ftp://adc.gsfc.nasa.gov/pub/adc/archives/catalogs/1/1157/ XZ catalog of zodiacal stars: ftp://adc.gsfc.nasa.gov/pub/adc/archives/catalogs/1/1201/ Tycho Reference Catalog: ftp://adc.gsfc.nasa.gov/pub/adc/archives/catalogs/1/1250/ Tycho 2 catalog: ftp://adc.gsfc.nasa.gov/pub/adc/archives/catalogs/1/1259/ USNO A2.0 catalog (very large): ftp://adc.gsfc.nasa.gov/pub/adc/archives/catalogs/1/1252/ HST Guide Star catalog 1.2: ftp://adc.gsfc.nasa.gov/pub/adc/archives/catalogs/1/1254/ HST Guide Star catalog 1.3: ftp://adc.gsfc.nasa.gov/pub/adc/archives/catalogs/1/1255/ Tycho 2 catalog: ftp://adc.gsfc.nasa.gov/pub/adc/archives/catalogs/1/1259/ AC 2000.2 catalog: ftp://adc.gsfc.nasa.gov/pub/adc/archives/catalogs/1/1275/ The original ZC (Zodiacal Catalogue):http://stjarnhimlen.se/zc/http://sorry.vse.cz/~ludek/zakryty/pub.phtml#zc |
|
