Astronomy Answers
Positions in the Sky

Here it is explained how you can calculate the approximate position of the Sun or a planet as seen from another planet. Further down the page it is explained how you can calculate the approximate position of the Moon as seen from the Earth.

1. Positions of the Sun and the Planets

If you want to calculate the position of a celestial body in the sky or as seen from the Sun (heliocentric), then you need information about that body and perhaps about the Earth. Which information you need depends on what exactly you want to calculate and how accurate the results should be. This page explains how you can calculate the position in the sky as seen from a particular place on Earth or from the Sun for a celestial body that orbits the Sun in a circle or ellipse. For this, you need the following information:

• the length a of the semimajor axis of the orbit, which indicates the size of the orbit. a should be measured in Astronomical Units (AU). Sometimes the perihelion distance q is given instead of a. Then you can calculate a as follows:

  (Eq. 1) a = q/(1 − e)

• the eccentricity e of the orbit, which indicates the shape of the orbit.
• the argument ω [omega] of the perihelion, which indicates how far the perihelion lies beyond the ascending node of the orbit. Warning: The symbol ω is also sometimes used for the sum of "our" ω and Ω. The symbol ϖ is also sometimes used for that sum, and that symbol resembles ω. We use Π = ω + Ω for that sum.
• the ecliptic longitude Ω [large omega] of the ascending node of the orbit.
• the inclination i of the orbit, which indicates how far the orbit is tilted compared to the ecliptic.

• the mean anomaly M0 at a specific date d0, which indicates the position of the planet in its orbit at that date. Sometimes L0 = M0 + Ω + ω is given instead of M0, and sometimes a date dperihelion is given instead when the object was in its perihelion (where M is equal to 0). Here are the transformation formulas:

  (Eq. 2) M0 = L0 − Ω − ω

  (Eq. 3) M0 = n (d0 − dperihelion)

In addition, you need the following information:

• the obliquity ε [epsilon] of the ecliptic
• your geographical longitude l (measured westward from Greenwich) and latitude φ [phi] (measured northward from the equator)
• your time zone tz, measured as the number of hours that you have to add to your local clock time t to get Universal Time (UTC)

All in all, you need (besides the desired date and time) 16 separate numbers: 6 to describe the orbit and position of the planet in space, 6 to describe the orbit and position of the Earth in space, and 4 that describe the orientation of the Earth and your position on Earth. If you only want to calculate the heliocentric position, then you need fewer numbers.

The orbital elements of the planets (at noon UTC on 1 January 2000, at Julian Day 2451545) are listed in the following table (with a measured in AU, and i, ω, Ω, and M0 in degrees relative to the equinox of that date).

For the Sun, a is equal to 0, so also xSun = ySun = zSun = 0 so you can skip steps 1 through 4 below.

The inclination of the orbit of the Earth relative to the equinox of the date is zero by definition, so that orbit has no nodes and the value of the ecliptic longitude Ω of the ascending node is immaterial. The value of Ω that I provide for the Earth in the table is the average of the values that you find a few years before and after 1 January 2000, when the inclination (compared to the equinox of J2000.0) is different from zero.

If you multiply the inclination i by −1, add 180 degrees to Ω, and suitably adjust the argument ω of the perihelion, then you still have the same orbit (but with the roles of ascending and descending nodes exchanged). For the other planets than the Earth, the inclination is very different from zero, so you can always choose it to be positive (by applying the factor of −1 when needed). Only for the orbit of the Earth can be inclination be positive or negative, and one source can choose to take the future inclination to be positive while another source can take the past inclination to be positive.

When the inclination is very small, then the length of the ascending node is not determined very accurately, so the values that different sources give for the length of the ascending node of the Earth can differ by many degrees (even if you take a possible sign difference for the inclination into account) but this does not make much difference to the position of the Earth that you calculate using that value, because the ecliptic latitude cannot be greater (in absolute sense) than the inclination.

Here is a table with some more handy numbers for each planet that depend only on the orbital elements. n is measured in degrees per day, a(1−e2) in AU, and Π in degrees. These numbers are explained later.

The steps you must follow to calculate the position of a planet as seen from Earth at date d are then:

1. Calculate the mean anomaly M from d, M0, d0 and a.
2. Calculate the true anomaly ν [nu] from M and e.
3. Calculate the distance r to the Sun, from a, e and ν.
4. Calculate the rectangular (= cartesian) heliocentric ecliptic coordinates (xplanet, yplanet, zplanet) from ν, r, ω, Ω and i.
5. Calculate the rectangular geocentric ecliptic coordinates (x, y, z) from (xplanet, yplanet, zplanet) and (xearth, yearth, zearth).
6. Calculate the geocentric ecliptic longitude λ [lambda] and latitude β [beta] from (x, y, z).
7. Calculate the right ascension α [alpha] and declination δ [delta] from λ, β and the obliquity ε [epsilon] of the ecliptic.
8. Calculate the sidereal time θ [theta] from t, your geographical longitude l and your time zone tz.
9. Calculate the hour angle H from α and θ.
10. Calculate the height h and the azimuth A from H, δ and your geographical latitude φ [phi].
11. You can calculate the elongation ψ [psi] of the planet from the geocentric ecliptic coordinates (cartesian or polar) of the Sun and the planet.

If you want to calculate the heliocentric position of a planet, then you should take (xplanet, yplanet, zplanet) for (x, y, z) in step 5.

These steps are explained in more detail below. As an example, we'll calculate the position of Jupiter at 01:00 hours Central European Standard Time (0:00 hours UTC) on 1 January 2004 from 52° north latitude and 5° east longitude (which is near Utrecht in the Netherlands).

The section about approximations details approximations that provide less accurate results with less calculations.

WARNING: These calculations are based on the assumption that the planets orbit around the Sun in elliptical orbits that never change. In reality these orbits do change, but slowly and not by much. These formulas are very suitable to determine whether a particular planet will be visible at a certain time, or what star it will be close to in the sky, but not (for example) to calculate whether a particular planet will eclipse a particular star. The accuracy of the results of these formulas is discussed later.

WARNING 2: Make sure that you use the right units for the trigonometric functions such as sin and arctan. In computer programs those functions commonly work with radians. Many calculators can be configured to work with degrees or radians. Calculate sin(1) with your calculator or computer program. If the result is approximately 0.8415, then the unit is the radian. If the result is approximately 0.01745, then the unit is the degree. To transform from degrees to radians, multiply the number of degrees by π/180 ≈ 0.01745329.

1.1. The Mean Anomaly

First, calculate which angle n the planet traverses on average per day, as seen from the Sun:

(Eq. 4) n = 0,9856076686°/(a√a)

and then the mean anomaly M for date d (measured in days) from the start date d0, the mean anomaly M0 on the start date, and n:

(Eq. 5) M = M0 + n (d − d0)

For d and d0 you could very well use Julian Day Numbers. The values of n for the planets are shown in the second table at the top of this page.

1.2. The True Anomaly

The true anomaly ν [nu] is the angle between the line from the focus of the orbit (the Sun) to the perihelion of the orbit and the line from the focus to the planet. To calculate the true anomaly, you need to solve the Equation of Kepler. See the Page about the Equation of Kepler about this.

1.3. The Distance to the Sun

The eccentricity e indicates the shape of the orbit and cannot be negative. If e is equal to 0, then the orbit is a circle. If e is greater than 0 but less than 1, then the orbit is an ellipse. If e is equal to 1, then the orbit is a parabola or a straight line. If e is greater than 1, then the orbit is a hyperbola.

The size a of the semimajor axis of the orbit indicates the size of the orbit. For a circular orbit, a is equal to the distance between the planet and the Sun. For an elliptic orbit, a is equal to half of the length of the longest line segment that fits within the ellipse.

The distance r of the planet from the Sun is:

(Eq. 6) r = a (1 − e2)/(1 + e cos ν)

The values of a(1−e2) for the planets are shown in the second table at the top of this page.

1.4. The Rectangular Heliocentric Ecliptical Coordinates

We can now calculate the rectangular heliocentric ecliptic coordinates, all measured in AU. Of those, xplanet is the distance from the Sun as measured along the line from the Sun to the vernal equinox, yplanet the distance to the line from the Sun to the vernal equinox, measured in the plane of the ecliptic, and zplanet the distance above the plane of the ecliptic.

(Eq. 7) xplanet = r (cos Ω cos(ω + ν) − sin Ω cos i sin(ω + ν))

(Eq. 8) yplanet = r (sin Ω cos(ω + ν) + cos Ω cos i sin(ω + ν))

(Eq. 9) zplanet = r sin i sin(ω + ν)

1.5. The Rectangular Geocentric Ecliptical Coordinates

If you have done steps 1 through 4 for the planet and also for the Earth, then you have calculated xplanet, yplanet, zplanet, xEarth, yEarth, zEarth. Then we can calculate the coordinates of the planet relative to the Earth, measured in AU:

(Eq. 10) x = xplanet − xEarth

(Eq. 11) y = yplanet − yEarth

(Eq. 12) z = zplanet − zEarth

1.6. The Geocentric Ecliptical Longitude and Latitude

The ecliptic latitude β [beta] indicates how far the celestial body is from the ecliptic. The ecliptic longitude λ [lambda] shows how far the celestial body is from the vernal equinox, measured along the ecliptic. For the distance Δ [large delta] of the planet from the Earth we find

(Eq. 13) Δ = √(x² + y² + z²)

and then

(Eq. 14) x = Δ cos λ cos β

(Eq. 15) y = Δ sin λ cos β

(Eq. 16) z = Δ sin β

It follows that

(Eq. 17) λ = arctan(y, x)

(Eq. 18) β = arcsin(z/Δ)

Here, arctan(y, x) is the special arc tangent function with two arguments that calculates the angle between the x axis and the line that goes from the point with coordinates (0,0) to the point with coordinates (x,y). This special arc tangent function is usually not available under that name on electronic calculators, but is often available in programming languages (with a name like "atan2"). Some calculators have a function that transforms rectangular coordinates to polar ones (a distance and an angle), and you can use that to calculate the angle that goes with x and y, and that angle is λ. That function often has a name like "R→P" or "Pol". To check: if you enter x = 1 and y = 0, then it should yield 0 for the angle, if you enter x = 0 and y = 1 then the result should be 90°, and if you enter x = −1 and y = −1 then the result should be −135° or 225° (those differ by 360° so they indicate the same angle).

If that function is not on your calculator or in your programming language, then you can use the ordinary arc tangent function, but then you have to do some extra work. In that case, calculate

(Eq. 19) λ0 = arctan(y/x)

If sin λ0 has the same sign (plus or minus) as y and cos λ the same sign as x, then λ = λ0. Otherwise, λ = λ0 + 180°.

1.7. The Right Ascension and Declination

The axis of rotation of the Earth is not at right angles to the plane of the orbit of the Earth, so the equator of the Earth makes an angle with the plane of the orbit of the Earth. This angle ε [epsilon] is called the obliquity of the ecliptic and its value at the beginning of 2000 was

(Eq. 20) ε = 23.4397°

The equatorial coordinate system in the sky is tied to the rotation axis of the Earth. The equatorial coordinates are the right ascension α [alpha] and the declination δ [delta]. The declination shows how far the body is from the celestial equator and determines from which parts of the Earth the object can be visible. The right ascension shows how far the body is from the vernal equinox, as measured along the celestial equator, and determines (together with other things) when the object is visible. We have:

(Eq. 21) sin α cos δ = sin λ cos β cos ε − sin β sin ε

(Eq. 22) cos α cos δ = cos λ cos β

(Eq. 23) δ = arcsin(sin β cos ε + cos β sin ε sin λ)

You can combine formulas 21 and 22 to:

(Eq. 24) α = arctan(sin λ cos ε − tan β sin ε, cos λ)

where the special arc tangent function was used again that was described at step 6. The right ascension is usually written as a time, with hours, minutes, and seconds. To transform it from a time measured in hours to an angle measured in degrees (so you can use it in formulas with things such as sine functions), you should multiply it by 15.

1.8. The Sidereal Time

Where a celestial body appears to be in the sky for you depends on the orientation of the Earth at your location compared to the stars. This angle is captured in the sidereal time θ (theta). The sidereal time at a certain moment is equal to the right ascension that transits (culminates, passes through the celestial meridian) at that moment. If t is the local clock time and tz is how much you have to add to the local time to get Universal Time (UTC), both measured in hours, then the sidereal time measured in degrees is equal to

(Eq. 25) θ = MEarth + ΠEarth + 15° (t + tz) − l

where

(Eq. 26) ΠEarth = ΩEarth + ωEarth

and ΠEarth is shown in the second table. The tz is equal to −1 hour for Dutch and Belgian standard time (CET, Central European Time), and to −2 hours for Dutch and Belgian daylight savings time (CEDT, Central European Daylight Savings Time). It is usually close to l/15° hours (reduced by one extra hour for daylight savings time).

The sidereal time is often written as a time rather than an angle, just like for the right ascension.

1.9. The Hour Angle

The hour angle H indicates how far the celestial body has passed beyond the celestial meridian. If the hour angle is zero, then the body transits and is highest in the sky, on the celestial meridian. The hour angle is often written as a time rather than an angle, just like for the sidereal time and the right ascension. That time is approximately equal to how long ago the body was highest in the sky.

(Eq. 27) H = θ − α

1.10. The Height and Azimuth

The position of a celestial body in the sky is specified by its altitude h above the horizon, and its azimuth A. The altitude is 0° at the horizon, +90° in the zenith (straight over your head), and −90° in the nadir (straight down). The azimuth is the direction along the horizon, which we measure from south to west. South has azimuth 0°, west +90°, north +180°, and east +270° (or −90°, that's the same thing). The altitude and azimuth are the horizontal coordinates. If the Earth had no atmosphere, then we'd have, for the horizontal coordinates:

(Eq. 28) sin A cos h = sin H cos δ

(Eq. 29) cos A cos h = cos H cos δ sin φ − sin δ cos φ

(Eq. 30) sin h = sin φ sin δ + cos φ cos δ cos H

It follows that

(Eq. 31) A = arctan(sin H, cos H sin φ − tan δ cos φ)

(Eq. 32) h = arcsin(sin φ sin δ + cos φ cos δ cos H)

where the special arc tangent function is used which was described earlier.

The atmosphere of the Earth bends light upward a bit so that it seems as if celestial bodies are a bit higher in the sky than they would be without an atmosphere. This effect is on average about 0.57° at the horizon and decreases rapidly with height. Moreover, it varies significantly close to the horizon, depending on the circumstances in the atmosphere. If you want to apply a correction to h for this refraction, then you can use the following formula (if you measure h in degrees):

(Eq. 33) happarent = h + 0.017/tan(h + 10.26/(h + 5.10))

1.11. The Elongation

The elongation ψ of a planet is a measure for the distance of that planet from the Sun along the sky. You can measure that distance (1) directly between the planet and the Sun, or (2) only along the ecliptic. These two methods in general give slightly different values, except for planets that are exactly on the ecliptic.

If you know the ecliptic longitude λSun of the Sun and also the ecliptic longitude λ and latitude β of the planet, then you can calculate elongation (1) with

(Eq. 34) ψ1 = arccos(cos β cos(λ − λSun))

Elongation (2) is then equal to

(Eq. 35) ψ2 = λ − λSun

It is customary to add or subtract multiples of 360° to or from ψ2 until the value lies between −180° and 180°. Elongation (1) can only be positive or zero, but elongation (2) can be negative, too. Positive means "east of the Sun" and negative "west of the Sun".

If you know the cartesian heliocentric ecliptic coordinates (xEarth, yEarth, zEarth) of Earth and the cartesian geocentric ecliptic coordinates (x, y, z) of the planet, then you can calculate elongation (1) also like this:

(Eq. 36) Δ = √(x2 + y2 + z2)

(Eq. 37) rEarth = √(xEarth2 + yEarth2 + zEarth2)

(Eq. 38) ψ1 = arccos(−(xEarth x + yEarth y + zEarth z)/(rEarth Δ))

1.12. Accuracy

How accurate are the results that these formulas yield? I compared those results for every day from 1 January 1980 through 1 January 2020 with the results that you get if you do the calculations using a truncated version of the VSOP theory according to the book "Astronomical Algorithms" by J. Meeus, who includes the effects of the mutual gravitational attraction of the planets, of abberation, nutation, the travel time of light, and the difference between UTC and TDT. The greatest differences are shown in the next table.

1.13. Application

1.13.1. Culmination or Transit

The culmination or transit of a celestial body is the moment at which the body passes through the celestial meridian and is highest in the sky. The hour angle H of the body is then 0. From this it follows that

(Eq. 39) θtransit = α

so

(Eq. 40) ttransit = (α + l − MEarth − ΠEarth)/15 − tz

if you measure the angles in degrees and the times in hours. With the values from the first table, this becomes

(Eq. 41) ttransit = (α + l − MEarth)/15 + 17h8m15s

if you measure ttransit in Universal Time (UTC). For Central European Standard Time (which is valid in winter in most of western Europe, including the Netherlands and Belgium) you should add 1 hour, but for Daylight Savings Time 2 hours.

1.13.2. Rise and Set

When the celestial body rises or sets, then it has a particular height h as calculated for an observer on a smooth spherical Earth without an atmosphere and for a celestial body that looks like a point. The Earth does in fact have an atmosphere, and the other conditions need not be met, either, so h need not be equal to 0 when the celestial body rises or sets. We therefore define that the rise or set happens when h is equal to some fixed value h0. This yields:

(Eq. 42) Hhorizon = arccos((sin h0 − sin φ sin δ)/(cos φ cos δ))

so

(Eq. 43) trise = ttransit − Hhorizon/15

(Eq. 44) tset = ttransit + Hhorizon/15

if you measure the times in hours and the hour angle in degrees.

1.14. Approximations

There are several approximations that one can make to simplify the calculations:

• approximate the elliptical orbit by a circular one (i.e., set e equal to 0). Then ν = M and r = a so steps 2 and 3 have become very simple.

• assume that the orbit is nearly a circle, so e is much smaller than 1. Then, approximately,

  (Eq. 45) ν ≈ M + 2 e sin M + (5⁄4) e2 sin 2M

  and then step 2 has become simpler.

• ignore the angle between the plane of the orbit and the plane of the orbit of the Earth (i.e., set i equal to 0). Then step 4 becomes much simpler, namely

  (Eq. 46) xplanet = r cos(Ω + ω + ν)

  (Eq. 47) yplanet = r sin(Ω + ω + ν)

  (Eq. 48) zplanet = 0

  where Π = Ω + ω and the value of Π for the planets is displayed in the second table.

  This held for the Earth anyway (according to our model), so then we also have β equal to 0, and then step 7 becomes

  (Eq. 49) α = arctan(sin λ cos ε, cos λ)

  (Eq. 50) δ = arcsin(sin ε sin λ)

• approximate sines and cosines of ε with polynomials. That only yields short formulas if you also assume that i is equal to 0 (which holds, for example, for the Sun as seen from Earth). Then, approximately,

  (Eq. 51) α ≈ λ − (¼) ε2 sin 2λ

  (Eq. 52) δ ≈ ε sin λ

• assume that h0 is much smaller than 10°. Then, approximately,

  (Eq. 53) Hhorizon ≈ arccos(−tan δ tan φ) − h0/√((cos δ)2 − (sin φ)2)

The next table shows the largest differences between the results that you get using some of these approximations and the results if you use the complete formulas mentioned before the section about approximations. The approximations of which the errors are shown here are (1) assume that the orbit is almost a circle, (2) ignore the angle between the orbit of the planet and the orbit of the Earth, and (3) approximate the sines and cosines of ε with polynomials. In general, approximation 1 yields the greatest reduction in the calculation time.

For example: if you calculate the declination δ of Mecury without approximations 1 - 3 (but with the approximation that its orbit is a fixed ellipse, which is assumed everywhere on this page), then the greatest difference relative to the results of the VSOP theory is 0.045°. If you use only approximation 1, then it becomes 0.23°. If you use approximations 1 and 2, then it becomes 4.9°. If you use approximations 1 - 3 then the greatest difference is 5.0°.

2. Position of the Moon

The position of the Moon cannot be calculated in the same fashion as the positions of the planets, because the orbit of the Moon around the Earth itself spins around with a period of only about 18 years, so that the nodes of the orbit of the Earth shift once around the whole ecliptic in that period. If you want to calculate the position of the Moon very accurately, then you must take into account thousands of terms and add them up. Here, I ignore all of those terms except for the very largest ones.

Use the numbers from the next table to calculate the mean geocentric ecliptic longitude L compared to the equinox of the date, the mean anomaly M, and the mean distance F of the Moon from its ascending node, measured in degrees. The value of such a quantity is equal to c0 + c1 (d − d0) where c0 and c1 come from the table and d − d0 is the number of days since 1 January 2000, 12:00 UTC, just like in the calculations for the planets.

Now calculate the geocentric ecliptic coordinates λ (longitude), β (latitude), and Δ (distance), with the angles in degrees and the distance in kilometers:

(Eq. 54) λ = L + 6.289 sin M

(Eq. 55) β = 5.128 sin F

(Eq. 56) Δ = 385001 − 20905 cos M

You can now use step 7 and later ones from the method for the planets to calculate the right ascension and declination or the height and azimuth.

I calculated the position of the Moon (λ, β, Δ) for every day between 1950 and 2050 according to the above method and also according to a much more accurate method (of which the above method includes only the very largest terms). The largest error in λ was 2.57 degrees (standard deviation 1.04 degrees), that in β 0.81 degrees (standard deviation 0.31 degrees), and that in Δ 7645 km (standard deviation 3388 km).

languages: [en] [nl]

http://aa.quae.nl/en/reken/hemelpositie.html;
Last updated: 2012-01-13