\( \newcommand{\Matrix}[1]{\left( \begin{matrix} #1 \end{matrix} \right)} \)
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}) \label{eq:M} \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.
The angles \( i \), \( ω \), and \( Ω \) depend on the orientation of the chosen coordinate system. In this table they are given relative to the ecliptic coordinate system of the Earth.
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 Table 2.
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 \( ν_\text{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 Table 2.
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(273.867° + 144.637°) \\ & − \sin(100.464°) × \cos(1.303°) × \sin(273.867° + 144.637°)) = −5.04289 \\ y_\text{Jupiter} & = 5.40406 × (\sin(100.464°) × \cos(273.867° + 144.637°) \\ & + \cos(100.464°) × \cos(1.303°) × \sin(273.867° + 144.637°)) = 1.93965 \\ z_\text{Jupiter} & = 5.40406 × \sin(1.303°) × \sin(273.867° + 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 at 0 UTC
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{eqnarray} θ & = & Θ − l \pmod{360°} \label{eq:sterrentijd} \\ Θ & = & M_\text{Earth} + Π_\text{Earth} + 15° (t + t_\text{z}) \pmod{360°} \label{eq:Θ} \end{eqnarray}
where \( Θ \) is the sidereal time on the prime meridian, and
\begin{equation} Π_\text{Earth} = Ω_\text{Earth} + ω_\text{Earth} \end{equation}
and \( Π_\text{Earth} \) is shown in Table 2. 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 \( \pmod{360°} \) means that you can add or subtract multiples of 360°.
The sidereal time is often written as a time rather than an angle, just like for the right ascension.
Eq. \eqref{eq:Θ} is in terms of time \( t \) in units of 1/24th of a mean solar day. On the Earth that unit is equal to 1 hour. To use that formula on other planets, you'd have to measure \( t \) in units of 1/24th of the mean solar day of that planet.
An alternative for Eq. \eqref{eq:Θ} is
\begin{equation} Θ = θ_0 + θ_1 (J − J_{2000}) \pmod{360°} \label{eq:ΘJ} \end{equation}
where \( J_{2000} = 2451545 \), and for the Earth \( θ_0 = 280.1470° \) and \( θ_1 = 360.9856235°/\text{d} \).
With the value of \( Π_\text{Earth} \) from Table 2, formula \ref{eq:Θ} becomes
\[ Θ = M_\text{Earth} + 102.937° + 15° (t + t_\text{z}) \pmod{360°} \]
For 01:00 hours Central European Standard Time we have \( t = 1 \) and \( t_\text{z} = −1 \). In step 1 we found that \( M_\text{Earth} = 357.009° \). All in all, we find
\begin{eqnarray*} Θ & = & 357.009° + 102.937° + 15° (1 + (−1)) = 99.946° = \text{6h39m47s} \\ θ & = & 99.946° − (−5°) = 104.946° \end{eqnarray*}
(after subtracting multiples of 360° or 24h).
For the chosen date and time we had \( J − J_{2000} = 1460.5 \), so Eq. \eqref{eq:ΘJ} yields
\begin{eqnarray} Θ & = & 280.1470° + 360.9856235° × 1460.5 = 527494.650 = 99.650° \pmod{360°} \\ θ & = & 99.650° − (−5°) = 104.650° \end{eqnarray}
This result differs from that of Eq. \eqref{eq:Θ}. The difference is likely due to the different sources of the orbital elements.
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 Table 1, 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.
The subsolar point is the point on Earth from where the Sun appears to stand straight overhead, in the zenith. What are the geographical coordinates of that point?
We assume that the equatorial coordinates \( α, δ \) of the Sun as seen from Earth are already available. See above for how to calculate them.
The reverse of Equations \eqref{eq:H} and \eqref{eq:sina}ff is
\begin{align} \sin H \cos δ & = \sin A \cos h \\ \cos H \cos δ & = \cos A \sin φ \cos h + \sin h \cos φ \\ \sin δ & = \sin φ \sin h − \cos φ \cos h \cos A \\ H & = \arctan(\sin A, \cos A \sin φ + \tan h \cos φ) \\ θ & = α + H \end{align}
In the zenith \( h = 90° \) so \( \sin h = 1 \) and \( \cos h = 0 \) and
\begin{eqnarray} \sin δ & = & \sin φ \\ \sin H \cos δ & = & 0 \\ \cos H \cos δ & = & \cos φ \end{eqnarray}
from which follow \( φ = δ \) and \( θ = α \).
Rearranging Eq. \eqref{eq:sterrentijd} leads to
\begin{equation} l = Θ − θ \pmod{360°} \end{equation}
For the Sun we had \( α = 280.710° \) and \( δ = −23.074° \), so for the sub-solar point \( φ = δ = −23.074° \) and \( θ = α = 280.710° \).
For the Earth we had \( Θ = 99.650° \), so
\[ l = 99.650° − 280.710° = 178.940° \pmod{360°} \]
so the subsolar point on Earth at that time is at latitude South 23.074° and longitude West 178.940°.
At what point on planet Y does object X appear straight overhead, in the zenith? This is the same question as in Sec. 2.3, but with object X instead of the Sun and planet Y instead of the Earth. We calculate this in the same way as for the Sun and Earth, except that we need to know the equatorial coordinates of X relative to the equatorial coordinate system tied to the equator of planet Y instead of to the Earth, and likewise for the sidereal time.
For example, what is the sub-Earth point on Mars on 0 UTC at the beginning of January 1st, 2004? In other words, from what location on Mars does the Earth appear to stand in the zenith then?
We adopt the notation from the page about the position of the Sun. A subscript \( _q \) indicates the Earth-based equatorial coordinate system. A subscript \( _Q \) indicates the planet-Y-based equatorial coordinate system. Likewise, subscripts \( _c \) and \( _C \) indicate the Earth-based and planet-Y-based ecliptic coordinate systems.
Note that "Earth-based" and "planet-Y-based" do not imply that the origin of the coordinate system must be at the center of the Earth or planet Y, but only that the primary plane of the coordinate system is parallel to the relevant plane that is associated with the Earth or planet Y.
If all we have are the equatorial coordinates \( α_q, δ_q \) of object X as seen from planet Y, measured in the Earth-based equatorial coordinate system, then we proceed as follows:
Find or calculate \( θ_0 \), \( θ_1 \), and \( Π \) for the planet, and also \( C_{Qq} \) for transforming from Earth-based equatorial coordinates to planet-Y-based equatorial coordinates, as described on the aforementioned page.
On the other page we find \( C_{cq} \) in its Eq. 61. Using \( ε_0 = 23.4392911° \) for the obliquity of the ecliptic for Earth,
\[ C_{cq} = \Matrix{1 & 0 & 0 \\ 0 & +0.9174821 & +0.3977772 \\ 0 & −0.3977772 & +0.9174821} \]
Just after Eq. 79 on that page we find that for Mars
\[ C_{Qc} = \Matrix{−0.0864092 & −0.9960834 & −0.01874317 \\ +0.8907726 & −0.0688209 & −0.4492080 \\ +0.4461587 & −0.0555116 & +0.8932306} \]
Combined, these yield
\[ C_{Qq} = C_{Qc} C_{cq} = \Matrix{−0.0864092 & −0.9064330 & −0.4134157 \\ +0.8907726 & +0.1155427 & −0.4395157 \\ +0.4461587 & −0.4062376 & +0.7974418} \]
We find \( Π \) for Mars just after Eq. 76 on that page:
\[ Π = 71.0041° \]
We find \( θ_0, θ_1 \) for Mars in Table 5 on that page:
\begin{eqnarray*} θ_0 & = & 313.3827° \\ θ_1 & = & 350.89198226°/\text{d} \end{eqnarray*}
Calculate the rectangular coordinates \( \vec{r}_q \) of a unit vector in the direction indicated by \( α_q, δ_q \), relative to the Earth-based equatorial coordinate system \( q \):
\begin{equation} \vec{r}_q = \Matrix{r_{q1} \\ r_{q2} \\ r_{q3}} = \Matrix{ \cos α_q \cos δ_q \\ \sin α_q \cos δ_q \\ \sin δ_q } \end{equation}
From Table 3 we see that the Earth-based equatorial coordinates of Mars at that time are \( α_q = 8.335° \), \( δ_q = +3.660° \). Someone looking at Earth from Mars is looking in exactly the opposite direction than someone looking at Mars from the Earth, so the Earth-based equatorial coordinates of the Earth as seen from Mars are \( α_q = 188.335° \), \( δ_q = −3.660° \).
Then
\[ \vec{r}_q = \Matrix{ \cos(188.335°) × \cos(−3.660°) \\ \sin(188.335°) × \cos(−3.660°) \\ \sin(−3.660°)} = \Matrix{−0.9874194 \\ −0.1446650 \\ −0.0638356} \]
Transform into the planet-Y-based equatorial coordinate system \( Q \):
\begin{equation} \vec{r}_Q = C_{Qq} \vec{r}_q \end{equation}
We find
\[ \vec{r}_Q = \Matrix{−0.0864092 & −0.9064330 & −0.4134157 \\ +0.8907726 & +0.1155427 & −0.4395157 \\ +0.4461587 & −0.4062376 & +0.7974418} \Matrix{−0.9874194 \\ −0.1446650 \\ −0.0638356} = \Matrix{+0.2428419 \\ −0.8682244 \\ −0.4326826} \]
Calculate the corresponding spherical coordinates relative to the planet-Y-based equatorial coordinate system:
\begin{eqnarray} δ_Q & = & \arcsin(r_{Q3}) \\ α_Q & = & \arctan(r_{Q2}, r_{Q1}) \end{eqnarray}
We find
\begin{eqnarray*} δ_Q & = & \arcsin(−0.4326826) = −25.638° \\ α_Q & = & \arctan(−0.8682244, +0.2428419) = 285.626° \end{eqnarray*}
Now use the equations from Sec. 2.3 to calculate the sub-X point on planet Y, reading "planet Y" where those equations refer to the Earth, and "object X" where they refer to the Sun.
We find geographical latitude \( φ = δ_Q = −25.638° \) and sidereal time \( θ = α_Q = 285.626° \). To calculate the geographical longitude \( l \) we need to know \( Θ \).
\begin{eqnarray*} Θ & = & θ_0 + θ_1 (J − J_{2000}) = 313.3827° + 350.89198226° × 1460.5 = 512791.122791° = 151.123° \pmod{360°} \\ l & = & Θ − θ = 151.123° − 285.626° = −134.503 = 225.497° \pmod{360°} \end{eqnarray*}
so the sub-Earth point on Mars at that date and time is at Martian geographical latitude South 25.638°, and (West) longitude 225.497°.
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 table 2.
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).
//aa.quae.nl/en/reken/hemelpositie.html;
Last updated: 2017-10-21