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.
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:
\begin{equation} a = \frac{q}{1 - e} \end{equation}
the mean anomaly \( M_0 \) at a specific date \( d_0 \), which indicates the position of the planet in its orbit at that date. Sometimes \( L_0 = M_0 + Ω + ω \) is given instead of \( M_0 \), and sometimes a date \( d_\text{perihelion} \) is given instead when the object was in its perihelion (where \( M \) is equal to 0). Here are the transformation formulas:
\begin{align} M_0 & = L_0 - Ω - ω \\ M_0 & = n (d_0 - d_\text{perihelion}) \end{align}
In addition, you need the following information:
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 \( M_0 \) in degrees relative to the equinox of that date). The numbers in the table come from [Meeus]. For the major planets they are based on the VSOP87 model, and for Pluto they are based on calculations by Goffin, Meeus, and Steyaert from 1986.
Table 1: Planets: Orbital Elements
\({a}\) | \({e}\) | \({i}\) | \({ω}\) | \({Ω}\) | \({M_0}\) | |
---|---|---|---|---|---|---|
Mercury | 0.38710 | 0.20563 | 7.005 | 29.125 | 48.331 | 174.795 |
Venus | 0.72333 | 0.00677 | 3.395 | 54.884 | 76.680 | 50.416 |
Earth | 1.00000 | 0.01671 | 0.000 | 288.064 | 174.873 | 357.529 |
Mars | 1.52368 | 0.09340 | 1.850 | 286.502 | 49.558 | 19.373 |
Jupiter | 5.20260 | 0.04849 | 1.303 | 273.867 | 100.464 | 20.020 |
Saturn | 9.55491 | 0.05551 | 2.489 | 339.391 | 113.666 | 317.021 |
Uranus | 19.21845 | 0.04630 | 0.773 | 98.999 | 74.006 | 141.050 |
Neptune | 30.11039 | 0.00899 | 1.770 | 276.340 | 131.784 | 256.225 |
Pluto | 39.543 | 0.2490 | 17.140 | 113.768 | 110.307 | 14.882 |
For the Sun, \( a \) is equal to 0, so also \( x_\text{Sun} = y_\text{Sun} = z_\text{Sun} = 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-e^2) \) in AU, and \( Π \) in degrees. These numbers are explained later.
Table 2: Planets: Other Orbital Numbers
\({n}\) | \({a(1-e^2)}\) | \({Π}\) | |
---|---|---|---|
Mercury | 4.092317 | 0.37073 | 77.456 |
Venus | 1.602136 | 0.72330 | 131.564 |
Earth | 0.985608 | 0.99972 | 102.937 |
Mars | 0.524039 | 1.51039 | 336.060 |
Jupiter | 0.083056 | 5.19037 | 14.331 |
Saturn | 0.033371 | 9.52547 | 93.057 |
Uranus | 0.011698 | 19.17725 | 173.005 |
Neptune | 0.005965 | 30.10796 | 48.124 |
Pluto | 0.003964 | 37.09129 | 224.075 |
The steps you must follow to calculate the position of a planet as seen from Earth at date \( d \) are then:
If you want to calculate the heliocentric position of a planet, then you should take \( (x_\text{planet}, y_\text{planet}, z_\text{planet}) \) 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.
First, calculate which angle \( n \) the planet traverses on average per day, as seen from the Sun:
\begin{equation} n = \frac{0.9856076686°}{a^{3/2}} \label{eq:n} \end{equation}
and then the mean anomaly \( M \) for date \( d \) (measured in days) from the start date \( d_0 \), the mean anomaly \( M_0 \) on the start date, and \( n \):
\begin{equation} M = M_0 + n (d - d_0) \end{equation}
For \( d \) and \( d_0 \) 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.
0:00 UTC at 1 January 2004 is 1460.6 days after the time for which \( M_0 \) from the table is valid, so \( d - d_0 \) is equal to 1460.5. The mean anomaly of Jupiter at the desired time is equal to \( M_\text{Jupiter} = 20.020 + 0.083056 × 1460.5 = 141.324° \) and for the Earth it is \( M_\text{Earth} = 357.529 + 0.985608 × 1460.5 = 357.009° \) (after subtracting multiples of 360°).
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.
For Jupiter, the true anomaly that goes with \( M_\text{Jupiter} = 141.324° \) and \( e_\text{Jupiter} = 0.04849 \) is equal to \( ν_\text{Jupiter} = 144.637° \). For the Earth, the true anomaly that goes with \( M_\text{Earth} = 357.009° \) and \( e_\text{Earth} = 0.01671 \) is equal to \( ν_earth = 356.907° \).
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:
\begin{equation} r = \frac{a (1 - e^2)}{1 + e \cos ν} \end{equation}
The values of \( a(1-e^2) \) for the planets are shown in the second table at the top of this page.
The distance between Jupiter and the Sun, measured in AU, is then
\[ r_\text{Jupiter} = \frac{5.19037}{1 + 0.04849 × \cos(144.637°)} = 5.40406 \]
and for the Earth
\[ r_\text{Earth} = \frac{0.99972}{1 + 0.01671 × \cos(356.907°)} = 0.98331 \]
We can now calculate the rectangular heliocentric ecliptic coordinates, all measured in AU. Of those, \( x_\text{planet} \) is the distance from the Sun as measured along the line from the Sun to the vernal equinox, \( y_\text{planet} \) the distance to the line from the Sun to the vernal equinox, measured in the plane of the ecliptic, and \( z_\text{planet} \) the distance above the plane of the ecliptic.
\begin{align} x_\text{planet} & = r (\cos Ω \cos(ω + ν) - \sin Ω \cos i \sin(ω + ν)) \\ y_\text{planet} & = r (\sin Ω \cos(ω + ν) + \cos Ω \cos i \sin(ω + ν)) \\ z_\text{planet} & = r \sin i \sin(ω + ν) \end{align}
For Jupiter we find
\begin{align*} x_\text{Jupiter} & = 5.40406 × (\cos(100.464°) × \cos(20.020° + 144.637°) \\ & - \sin(100.464°) × \cos(1.303°) × \sin(20.020° + 144.637°)) = −5.04289 \\ y_\text{Jupiter} & = 5.40406 × (\sin(100.464°) × \cos(20.020° + 144.637°) \\ & + \cos(100.464°) × \cos(1.303°) × \sin(20.020° + 144.637°)) = 1.93965 \\ z_\text{Jupiter} & = 5.40406 × \sin(1.303°) × \sin(20.020° + 144.637°) = 0.10478 \end{align*}
and for the Earth
\begin{align*} x_\text{Earth} & = 0.98331 × (\cos(174.873°) × \cos(288.064° \\ & + 356.907°) - \sin(174.873°) × \cos(0°) × \sin(288.064° + 356.907°)) = −0.16811 \\ y_\text{Earth} & = 0.98331 × (\sin(174.873°) × \cos(288.064° \\ & + 356.907°) + \cos(174.873°) × \cos(0°) × \sin(288.064° + 356.907°)) = 0.96884 \\ z_\text{Earth} & = 0.98331 × \sin(0°) × \sin(288.064° + 356.907°) = 0.00000 \end{align*}
If you have done steps 1 through 4 for the planet and also for the Earth, then you have calculated \( x_\text{planet}, y_\text{planet}, z_\text{planet}, x_\text{Earth}, y_\text{Earth}, z_\text{Earth} \). Then we can calculate the coordinates of the planet relative to the Earth, measured in AU:
\begin{align} x & = x_\text{planet} - x_\text{Earth} \\ y & = y_\text{planet} - y_\text{Earth} \\ z & = z_\text{planet} - z_\text{Earth} \end{align}
For the rectangular geocentric coordinates of Jupiter, measured in AU, we find
\begin{align*} x & = −5.04289 - (−0.16811) = −4.87477 \\ y & = 1.93965 - 0.96884 = 0.97081 \\ z & = 0.10478 - 0.00000 = 0.10478 \end{align*}
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
\begin{equation} Δ = \sqrt{x^2 + y^2 + z^2} \end{equation}
and then
\begin{align} x & = Δ \cos λ \cos β \\ y & = Δ \sin λ \cos β \\ z & = Δ \sin β \end{align}
It follows that
\begin{align} λ & = \arctan(y, x) \\ β & = \arcsin\left( \frac{z}{Δ} \right) \end{align}
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
\begin{equation} λ_0 = \arctan\left( \frac{y}{x} \right) \end{equation}
If \( \sin λ_0 \) has the same sign (plus or minus) as \( y \) and \( \cos λ \) the same sign as \( x \), then \( λ = λ_0 \). Otherwise, \( λ = λ_0 + 180° \).
We find
\begin{align*} Δ & = \sqrt{(−4.87477)^2 + 0.97081^2 + 0.10478^2} = 4.97161 \\ λ & = \arctan(0.97081, −4.87477) = 168.737° \\ β & = \arcsin\left( \frac{0.10478}{4.97161} \right) = 1.208° \end{align*}
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
\begin{equation} ε = 23.4397° \end{equation}
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:
\begin{align} \sin α \cos δ & = \sin λ \cos β \cos ε - \sin β \sin ε \label{eq:sinalpha} \\ \cos α \cos δ & = \cos λ \cos β \label{eq:cosalpha} \\ δ & = \arcsin(\sin β \cos ε + \cos β \sin ε \sin λ) \end{align}
You can combine formulas \ref{eq:sinalpha} and \ref{eq:cosalpha} to:
\begin{equation} α = \arctan(\sin λ \cos ε - \tan β \sin ε, \cos λ) \label{eq:alpha} \end{equation}
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.
For Jupiter we find
\begin{align*} α & = \arctan(\sin(168.737°) × \cos(23.4397°) \\ & - \tan(1.208°) × \sin(23.4397°), \cos(168.737°)) = 170.120° = \text{11h20m29s} \\ δ & = \arcsin(\sin(1.208°) × \cos(23.4397°) \\ & + \cos(1.208°) × \sin(23.4397°) × \sin(168.737°)) = 5.567° \end{align*}
If you calculate the right ascension, declination, and distance for the Sun and all planets from the table, then you find
Table 3: Planet Positions on 1-1-2004
Name | right ascension | declination | distance | |
---|---|---|---|---|
degrees | time | degrees | AU | |
Sun | 280.710 | 18h42m50s | −23.074 | 0.98331 |
Mercury | 268.693 | 17h54m46s | −20.296 | 0.70403 |
Venus | 316.189 | 21h04m45s | −18.614 | 1.3061 |
Mars | 8.335 | 0h33m20s | +3.660 | 1.1115 |
Jupiter | 170.120 | 11h20m29s | +5.567 | 4.9716 |
Saturn | 100.256 | 6h41m01s | +22.420 | 8.0443 |
Uranus | 333.148 | 22h12m36s | −11.868 | 20.654 |
Neptune | 313.525 | 20h54m06s | −17.459 | 30.973 |
Pluto | 260.277 | 17h21m07s | −14.497 | 31.700 |
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 \( t_\text{z} \) 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
\begin{equation} \label{eq:sterrentijd} θ = M_\text{Earth} + Π_\text{Earth} + 15° (t + t_\text{z}) - l \end{equation}
where
\begin{equation} Π_\text{Earth} = Ω_\text{Earth} + ω_\text{Earth} \end{equation}
and \( Π_\text{Earth} \) is shown in the second table. The \( t_\text{z} \) 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.
With the value of \( Π_\text{Earth} \) from the second table, formula \ref{eq:sterrentijd} becomes
\[ θ = M_\text{Earth} + 102.937° + 15° (t + t_\text{z}) - l \]
For 01:00 hours Central European Standard Time at 5° east longitude we have \( t = 1 \), \( t_\text{z} = −1 \) and \( l = −5° \). In step 1 we found that \( M_\text{Earth} = 357.009° \). All in all, we find
\[ θ = 357.009° + 102.937° + 15° (1 + (−1)) - (−5°) = 104.946° = \text{6h59m47s} \]
(after subtracting multiples of 360° or 24h).
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.
\begin{equation} \label{eq:H} H = θ - α \end{equation}
For Jupiter, we find
\[ H = 104.946° - 170.120° = −65.174° = \text{6h59m47s} - \text{11h20m29s} = \text{−4h20m42s} \]
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:
\begin{align} \sin A \cos h & = \sin H \cos δ \label{eq:sina} \\ \cos A \cos h & = \cos H \cos δ \sin φ - \sin δ \cos φ \label{eq:cosa} \\ \label{eq:h} \sin h & = \sin φ \sin δ + \cos φ \cos δ \cos H \end{align}
It follows that
\begin{align} A & = \arctan(\sin H, \cos H \sin φ - \tan δ \cos φ) \label{eq:a} \\ h & = \arcsin(\sin φ \sin δ + \cos φ \cos δ \cos H) \end{align}
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):
\begin{equation} h_\text{apparent} = h + \frac{0.017}{\tan\left( h + \frac{10.26}{h + 5.10} \right)} \end{equation}
For Jupiter at 01:00 hours Central European Standard Time on 1 January 2004 as seen from 52° north latitude and 5° east longitude we find
\begin{align*} h & = \arcsin(\sin(52°) \sin(5.567°) + \cos(52°) \cos(5.567°) \cos(−65.174°)) = 19.495° \\ A & = \arctan(\sin(−65.174°), \cos(−65.174°) \sin(52°) - \tan(5.567°) \cos(52°)) = −73.383° \end{align*}
Jupiter is then about 19.5 degrees above the horizon, about 16 degrees south of east. The correction for refraction is about 0.05°.
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 \( λ_\text{Sun} \) of the Sun and also the ecliptic longitude \( λ \) and latitude \( β \) of the planet, then you can calculate elongation (1) with
\begin{equation} ψ_1 = \arccos(\cos β \cos(λ - λ_\text{Sun})) \end{equation}
Elongation (2) is then equal to
\begin{equation} ψ_2 = λ - λ_\text{Sun} \end{equation}
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 \( (x_\text{Earth}, y_\text{Earth}, z_\text{Earth}) \) of Earth and the cartesian geocentric ecliptic coordinates \( (x, y, z) \) of the planet, then you can calculate elongation (1) also like this:
\begin{align} Δ & = \sqrt{x^2 + y^2 + z^2} \\ r_\text{Earth} & = \sqrt{x_\text{Earth}^2 + y_\text{Earth}^2 + z_\text{Earth}^2} \\ ψ_1 & = \arccos\left( -\frac{x_\text{Earth} x + y_\text{Earth} y + z_\text{Earth} z}{r_\text{Earth} Δ} \right) \end{align}
We found earlier for Jupiter that \( λ = 168.737° \) and \( β = 1.208° \). Similarly, we find for the Sun that \( λ_\text{Sun} = 279.844° \). It follows that
\begin{align*} ψ_1 & = \arccos(\cos(1.208°) × \cos(168.737° - 279.844°)) = 111.102° \\ ψ_2 & = 168.737° - 279.844° = −111.063° \end{align*}
Earlier, we found \( (x_\text{Earth,} y_\text{Earth}, z_\text{Earth}) = (−0.16811, 0.96884, 0.00000) \) and \( (x, y, z) = (−4.87477, 0.97081, 0.10478) \) and, from that, \( Δ = 4.97161 \) and \( r_\text{Earth} = 0.98331 \), and with these we now find
\[ ψ_1 = \arccos\left( -\frac{−0.16811 × −4.87477 + 0.96884 × 0.97081 + 0 × 0.10478}{4.97161 × 0.98331} \right) = 111.102° \]
just like before.
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.
Table 4: Planets: Position Calculation Accuracy
Name | \({α}\) | \({δ}\) | \({Δ}\) |
---|---|---|---|
degrees | degrees | AU | |
Sun | 0.03 | 0.01 | 0.0000 |
Mercury | 0.09 | 0.04 | 0.0013 |
Venus | 0.17 | 0.05 | 0.0008 |
Mars | 0.26 | 0.07 | 0.0018 |
Jupiter | 0.32 | 0.12 | 0.0093 |
Saturn | 1.08 | 0.43 | 0.049 |
Uranus | 1.00 | 0.35 | 0.047 |
Neptune | 0.68 | 0.2 | 0.072 |
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
\begin{equation} θ_\text{transit} = α \end{equation}
so
\begin{equation} t_\text{transit} = \frac{α + l - M_\text{Earth} - Π_\text{Earth}}{15} - t_\text{z} \end{equation}
if you measure the angles in degrees and the times in hours. With the values from the first table, this becomes
\begin{equation} t_\text{transit} = \frac{α + l - M_\text{Earth}}{15} + \text{17h8m15s} \end{equation}
if you measure \( t_\text{transit} \) 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.
For Jupiter on 1 January 2004 we found in step 7 that \( α \) is equal to 170.120°, and for the Earth we found in step 1 that \( M_\text{Earth} \) is equal to 357.009°. For 5° east longitude (\( l = −5) \) we then find
\[ t_\text{transit} = \frac{170.120 + (−5) - 357.009}{15} + \text{17h8m15s} = \text{4h20m42s} \]
so Jupiter culminates (is highest in the sky) at 04:21 hours Universal Time, which is equivalent to 05:21 hours Central European Standard Time.
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 \( h_0 \). This yields:
\begin{equation} H_\text{horizon} = \arccos\left( \frac{\sin h_0 - \sin φ \sin δ}{\cos φ \cos δ} \right) \end{equation}
so
\begin{align} t_\text{rise} & = t_\text{transit} - \frac{H_\text{horizon}}{15} \\ t_\text{set} & = t_\text{transit} + \frac{H_\text{horizon}}{15} \end{align}
if you measure the times in hours and the hour angle in degrees.
For Jupiter we found, in step 7, that \( δ \) is equal to 5.567°. For 52° north latitude (\( φ = 52° \)) and with \( h_0 = 0 \) we then find
\begin{align*} H_\text{horizon} & = \arccos(-\tan(52°) × \tan(5.567°)) = 97.167° = \text{6h28m40s} \\ t_\text{rise} & = t_\text{transit} - \frac{H_\text{horizon}}{15} = \text{5h20m42s} - \text{6h28m40s} = \text{−1h7m58s} = \text{22h52m2s} \\ t_\text{set} & = t_\text{transit} + \frac{H_\text{horizon}}{15} = \text{5h20m42s} + \text{6h28m40s} = \text{11h49m22s} \end{align*}
I added 24 hours to get the final result for \( t_\text{rise} \), so that time is one on the previous day. All in all, Jupiter is above the horizon from 22h52m to 11h49m CET. There is no point in quoting these times more precise than to the nearest minute.
There are several approximations that one can make to simplify the calculations:
assume that the orbit is nearly a circle, so \( e \) is much smaller than 1. Then, approximately,
\begin{equation} ν ≈ M + 2 e \sin M + \frac{5}{4} e^2 \sin 2M \end{equation}
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
\begin{align} x_\text{planet} & = r \cos(Ω + ω + ν) \\ y_\text{planet} & = r \sin(Ω + ω + ν) \\ z_\text{planet} & = 0 \end{align}
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
\begin{align} α & = \arctan(\sin λ \cos ε, \cos λ) \\ δ & = \arcsin(\sin ε \sin λ) \end{align}
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,
\begin{align} α & ≈ λ - \frac{1}{4} ε^2 \sin 2λ \\ δ & ≈ ε \sin λ \end{align}
assume that \( h_0 \) is much smaller than 10°. Then, approximately,
\begin{equation} H_\text{horizon} ≈ \arccos(-\tan δ \tan φ) - \frac{h_0}{\sqrt{\cos^2 δ - \sin^2 φ}} \end{equation}
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.
Table 5: Planets: Position Approximation Accuracy
\(α\) (degrees) | \(δ\) (degrees) | \(Δ\) (AU) | ||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|
1+2+3 | 1+2 | 1 | none | 1+2+3 | 1+2 | 1 | none | 1+2+3 | 1+2 | 1 | none | |
Sun | 0.14 | 0.032 | 0.032 | 0.032 | 0.26 | 0.011 | 0.011 | 0.011 | 0.0001 | 0.0001 | 0.0001 | 0.0001 |
Mercury | 2.3 | 2.3 | 0.45 | 0.088 | 5.0 | 4.9 | 0.23 | 0.045 | 0.0084 | 0.0084 | 0.0056 | 0.0013 |
Venus | 3.4 | 3.4 | 0.17 | 0.17 | 8.2 | 7.9 | 0.05 | 0.05 | 0.0045 | 0.0045 | 0.00078 | 0.00077 |
Mars | 2.6 | 2.6 | 0.41 | 0.26 | 6.3 | 6.5 | 0.14 | 0.074 | 0.0028 | 0.0028 | 0.0026 | 0.0017 |
Jupiter | 0.87 | 0.86 | 0.32 | 0.32 | 1.8 | 1.6 | 0.12 | 0.12 | 0.0095 | 0.0095 | 0.0094 | 0.0093 |
Saturn | 1.2 | 1.2 | 1.1 | 1.1 | 3.2 | 3.0 | 0.44 | 0.43 | 0.049 | 0.049 | 0.049 | 0.049 |
Uranus | 1.1 | 0.99 | 1.0 | 1.0 | 1.1 | 1.1 | 0.35 | 0.35 | 0.047 | 0.047 | 0.047 | 0.047 |
Neptune | 1.1 | 1.1 | 0.68 | 0.68 | 1.4 | 1.4 | 0.27 | 0.27 | 0.072 | 0.072 | 0.072 | 0.072 |
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°.
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 \( c_0 + c_1 (d - d_0) \) where \( c_0 \) and \( c_1 \) come from the table and \( d - d_0 \) is the number of days since 1 January 2000, 12:00 UTC, just like in the calculations for the planets.
\({c_0}\) | \({c_1}\) | |
---|---|---|
\({L}\) | 218.316 | 13.176396 |
\({M}\) | 134.963 | 13.064993 |
\({F}\) | 93.272 | 13.229350 |
Now calculate the geocentric ecliptic coordinates \( λ \) (longitude), \( β \) (latitude), and \( Δ \) (distance), with the angles in degrees and the distance in kilometers:
\begin{align} λ & = L + 6.289 \sin M \\ β & = 5.128 \sin F \\ Δ & = 385001 - 20905 \cos M \end{align}
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.
We calculate the position of the Moon for 0:00 UTC on 1 January 2004, which is 1460.5 days after \( d_0 \). Then
\begin{align*} L & = 218.316 + 13.176396 × 1460.5 = 19462.44° = 22.44° \\ M & = 134.963 + 13.064993 × 1460.5 = 19216.39° = 136.39° \\ F & = 93.272 + 13.229350 × 1460.5 = 19414.74° = 334.74° \end{align*}
This yields
\begin{align*} λ & = 22.44 + 6.289 × \sin(136.39°) = 26.78° \\ β & = 5.128 × \sin(334.74°) = −2.19° \\ Δ & = 385001 - 20905 × \cos(136.39°) = 400136 \text{km} \end{align*}
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).
http://aa.quae.nl/en/reken/hemelpositie.html;
Last updated: 2016-10-08