The position of the Sun in the sky as seen from a planet (such as the Earth) is determined by four things:
Below, I provide formulas that take all of these things into account, based on the rotation axes of the planets as defined by the IAU in "Report of the IAU/IAG Working Group on Cartographic Coordinates and Rotational Elements of the Planets and Satellites: 2000" by P.K. Seidelmann et al. (http://www.hnsky.org/iau-iag.htm). I neglected small effects here and there, to keep the formulas simple. Remember that you can add multiples of 360 degrees to any angle without changing its direction.
When in this page an ecliptic or pole or coordinates such as ecliptic longitude or declination are mentioned, then these are appropriate for the planet where the observer is. Such coordinates are based on the orbit and equator of that planet and are not identical to the coordinates of the same names that we use on Earth. For example, the ecliptic of Mars, the ecliptic that is appropriate for Mars, is not the same as the ecliptic of Earth, and an earth-based atlas of the stars is of no use if you have a right ascension and declination calculated for the Sun as seen from Mars.
For examples, we'll calculate the position of the Sun for 1 April 2004 at 12:00 UTC, as seen from 52° north latitude and 5° east longitude (Netherlands) on Earth, and as seen from 14°36' south latitude and 184°36' west longitude (Gusev crater) on Mars.
For this kind of astronomical calculations, it is convenient to express the date and time using an unending day numbering scheme. Such a scheme is provided by the Julian Day Number \( J \). The Julian Day Number Calculation Page explains how you can calculate the Julian Day Number for a date in the Gregorian calendar. For the calculations of the position of the Sun you should express time as Universal Time (UTC), and this includes the Julian Date. For example, JD 2453144.5 corresponds to 0:00 hours UTC on 19 May 2004, which is (for example) equivalent to 1:00 hours Central European Time on 19 May 2004, or 20:00 hours (8 pm) Eastern Standard Time on 18 May 2004.
Only on Earth do the seasons repeat themselves after about one calendar year (in the western ― Gregorian ― calendar), so, only for the Earth, the day number \( d \) in the calendar year can be used instead of the Julian Day Number: \( d \) = 1 corresponds to 0:00 UTC on January 1st, 2 to 0:00 UTC on January 2nd, 32 to 0:00 UTC on February 1st, and so on.
Because we see the Sun from the planet, we see the motion of the planet around the Sun reflected in the apparent motion of the Sun along the ecliptic, relative to the stars.
If the orbit of the planet were a perfect circle, then the planet as seen from the Sun would move along its orbit at a fixed speed, and then it would be simple to calculate its position (and also the position of the Sun as seen from the planet). The position that the planet would have relative to its perihelion if the orbit of the planet were a circle is called the mean anomaly, indicated in the formulas as \( M \).
You can calculate the mean anomaly of the planets (measured in degrees) to reasonable accuracy using the following formula, for a date given as a Julian Day Number (JD) \( J \):
\begin{align} M & = (M_0 + M_1×(J − J_{2000})) \bmod 360° \label{eq:m} \\ J_{2000} & = 2451545 \end{align}
You should take \( M_0 \) (in degrees) and \( M_1 \) (in degrees per day) from the following table:
Table 1: Planets: Mean Anomaly
\(M_0\) | \(M_1\) | |
---|---|---|
Mercury | 174.7948 | 4.09233445 |
Venus | 50.4161 | 1.60213034 |
Earth | 357.5291 | 0.98560028 |
Mars | 19.3730 | 0.52402068 |
Jupiter | 20.0202 | 0.08308529 |
Saturn | 317.0207 | 0.03344414 |
Uranus | 141.0498 | 0.01172834 |
Neptune | 256.2250 | 0.00598103 |
Pluto | 14.882 | 0.00396 |
For the Earth, you can also use the following formula:
\begin{equation} M = (−3.59° + 0.98560° × d) \bmod 360° \label{eq:m-aarde} \end{equation}
where \( d \) is the time since 00:00 UTC at the beginning of the most recent January 1st, measured in (whole and fractional) days.
The chosen date and time correspond to Julian Day Number 2453097, so \( J = 2453097 \), and also to \( d = 92.5 \). For the Earth, using equation \ref{eq:m} we find \( M_\text{earth} = 1887.1807° \bmod 360° = 87.1807° \) and with equation \ref{eq:m-aarde} we find \( M_\text{earth} = 87.58° \). Because equation \ref{eq:m} is a bit more accurate that equation \ref{eq:m-aarde}, we'll use equation \ref{eq:m}. For Mars, we find \( M_\text{mars} = 832.6531° \bmod 360° = 112.6531° \).
The orbits of the planets are not perfect circles but rather ellipses, so the speed of the planet in its orbit varies, and therefore the apparent speed of the Sun along the ecliptic also varies throughout the planet's year.
The true anomaly (symbol \( ν \), nu) is the angular distance of the planet from the perihelion of the planet, as seen from the Sun. For a circular orbit, the mean anomaly and the true anomaly are the same. The difference between the true anomaly and the mean anomaly is called the Equation of Center, written here as \( C \):
\begin{equation} ν = M + C \end{equation}
To calculate the Equation of Center or the true anomaly from the mean anomaly, you must first solve the Equation of Kepler. This equation cannot be solved in general, but you can find approximations to the solution that are more and more accurate. See the Calculation Page for the Equation of Kepler for more information about this. If the orbit looks more like a circle than like a parabola, then the following approximation is sufficiently accurate:
\begin{equation} C ≈ C_1 \sin M + C_2 \sin(2 M) + C_3 \sin(3 M) + C_4 \sin(4 M) + C_5 \sin(5 M) + C_6 \sin(6 M) \label{eq:c} \end{equation}
You can find the coefficients \( C_1 \) through \( C_6 \) in the following table. They depend on the eccentricity \( e \) of the planet's orbit. If a particular coefficient is not mentioned, then it is equal to zero (to four decimal places). The orbits of Mercury and Pluto deviate the most from circularity, so they need the most coefficients. The column \( E_C \) shows the maximum error that you get if you use the approximation with the coefficients from the table.
Table 2: Planets: Equation of Center
\(C_1\) | \(C_2\) | \(C_3\) | \(C_4\) | \(C_5\) | \(C_6\) | \(E_C\) | |
---|---|---|---|---|---|---|---|
Mercury | 23.4400 | 2.9818 | 0.5255 | 0.1058 | 0.0241 | 0.0055 | 0.0026 |
Venus | 0.7758 | 0.0033 | 0.0000 | ||||
Earth | 1.9148 | 0.0200 | 0.0003 | 0.0000 | |||
Mars | 10.6912 | 0.6228 | 0.0503 | 0.0046 | 0.0005 | 0.0001 | |
Jupiter | 5.5549 | 0.1683 | 0.0071 | 0.0003 | 0.0001 | ||
Saturn | 6.3585 | 0.2204 | 0.0106 | 0.0006 | 0.0001 | ||
Uranus | 5.3042 | 0.1534 | 0.0062 | 0.0003 | 0.0001 | ||
Neptune | 1.0302 | 0.0058 | 0.0001 | ||||
Pluto | 28.3150 | 4.3408 | 0.9214 | 0.2235 | 0.0627 | 0.0174 | 0.0096 |
For the Earth, we find
\begin{align*} C_\text{earth} & = 1.9148° × \sin(87.1807°) \\ & + 0.0200° × \sin(2×87.1807°) \\ & + 0.0003° × \sin(3×87.1807°) = 1.9142° \end{align*}
and thence \( ν_\text{earth} = 89.0949° \). For Mars we find
\begin{align*} C_\text{mars} & = 10.6912° × \sin(112.6531°) \\ & + 0.6228° × \sin(2×112.6531°) \\ & + 0.0503° × \sin(3 × 112.6531°) \\ & + 0.0046° × \sin(4 × 112.6531°) \\ & + 0.0005° × \sin(5 × 112.6531°) = 9.4092° \end{align*}
and \( ν_\text{mars} = 122.0623° \).
To find the position of the Sun in the sky we need to know what the ecliptic longitude \( Π \) is of the perihelion of the planet, relative to the ecliptic and vernal equinox (the ascending equinox) of the planet. The ecliptic of the planet is the plane of the orbit of the planet, which makes an angle with the orbit (ecliptic) of the Earth (and that angle is called the inclination of the orbit). The vernal equinox of the planet is the point where the Sun passes from south to north through the plane of the equator of the planet. We also need to know the obliquity \( ε \) of the equator of the planet compared to the orbit of the planet. These two values are listed in the following table for each planet, measured in degrees.
Table 3: Planets: Perihelion and Obliquity of the Ecliptic (Planetocentric)
\(Π\) | \(ε\) | |
---|---|---|
Mercury | 111.5943 | 0.02 |
Venus | 73.9519 | 2.64 |
Earth | 102.9372 | 23.45 |
Mars | 70.9812 | 25.19 |
Jupiter | 237.2074 | 3.12 |
Saturn | 99.4571 | 26.74 |
Uranus | 5.4639 | 82.22 |
Neptune | 182.1957 | 27.84 |
Pluto | 4.5433 | 57.46 |
These values are different than those given on the NASA Fact Sheet (http://nssdc.gsfc.nasa.gov/planetary/planetfact.html), because there they are measured relative to the ecliptic of the Earth (the plane of the orbit of the Earth), while here they are measured relative to the ecliptic of the planet (the plane of the orbit of the planet).
If you use the orbital elements of a planet that are measured relative to the ecliptic of the Earth, then the coordinates that you calculate from those elements are relative to the Earth's ecliptic, too. You can then meaningfully compare those coordinates to coordinates of stars and other objects in star atlases and catalogs made for use on Earth, to see where the planet is compared to those stars and other objects, but you cannot use them to calculate where the Sun is relative to the horizon on that planet.
On the other hand, if you use orbital elements that are measured relative to the ecliptic of the planet (as on this page), then you cannot meaningfully compare the resulting coordinates with star atlases and catalogs that were made for use on Earth, but you can use them to calculate where the Sun is relative to the horizon of that planet.
The ecliptical longitude \( λ \) (lambda) is the position along the ecliptic, relative to the vernal equinox (so relative to the stars). The mean longitude \( L \) is the ecliptical longitude that the planet would have if the orbit were a perfect circle. That is
\begin{equation} L = M + Π \label{L} \end{equation}
The ecliptic longitude of the planet as seen from the Sun is equal to
\begin{equation} λ = ν + Π = M + Π + C = L + C \end{equation}
If you look at the Sun from the planet, then you're looking in exactly the opposite direction than if you look at the planet from the Sun, so those directions are 180° apart. So, the ecliptic longitude of the Sun, as seen from the planet, is equal to
\begin{eqnarray} L_\text{sun} & = L + 180° = M + Π + 180° \\ λ_\text{sun} & = λ + 180° = ν + Π + 180° = M + Π + C + 180° = L_\text{sun} + C \label{eq:m2lambda} \end{eqnarray}
The value of \( λ_\text{sun} \) determines when the (astronomical) seasons begin: when \( λ_\text{sun} = 0° \), then spring begins in the northern hemisphere, and autumn in the southern hemisphere. Each next multiple of 90° brings the start of the next season.
The ecliptic latitude \( β_\text{sun} \) (beta) of the Sun, the perpendicular distance of the Sun from the ecliptic, is always so small that we can ignore it here. With this, we now have the ecliptic coordinates of the Sun.
As seen from Earth we find \( λ_\text{sun} = 372.0319° = 12.0321° \pmod{360°} \) and as seen from Mars \( λ_\text{sun} = 373.0435° = 13.0435° \pmod{360°} \).
The equatorial coordinate system in the sky is tied to the rotation axis of the planet. The equatorial coordinates are the right ascension \( α \) (alpha) and the declination \( δ \) (delta). The declination determines from which parts of the planet the object can be visible, and the right ascension determines (together with other things) when the object is visible.
With these formulas you can calculate the equatorial coordinates from the ecliptic coordinates:
\begin{align} \sin α \cos δ & = \sin λ \cos ε \cos β − \sin β \sin ε \label{eq:sinalpha} \\ \cos α \cos δ & = \cos λ \cos β \label{eq:cosalpha} \\ \sin δ & = \sin β \cos ε + \cos β \sin ε \sin λ \end{align}
For the Sun we have \( β_\text{sun} = 0 \), so then
\begin{align} α_\text{sun} & = \arctan(\sin λ_\text{sun} \cos ε, \cos λ_\text{sun}) \label{eq:alpha-zon} \\ δ_\text{sun} & = \arcsin(\sin λ_\text{sun} \sin ε) \label{eq:delta-zon} \end{align}
For future reference we define
\begin{equation} α_\text{sun} = λ_\text{sun} + S \label{eq:a} \end{equation}
If \( ε \) is close enough to 0° or 180°, and if we neglect small terms, then we can approximate the relationship between the right ascension \( α_\text{sun} \) and the ecliptic longitude \( λ_\text{sun} \) of the Sun as seen from the planet by
\begin{align} \arctan(\tan(λ) \cos(ε)) & = λ − \left( \frac{1}{4} ε^2 + \frac{1}{24} ε^4 + \frac{17}{2880} ε^6 \right) \sin(2 λ) \notag \\ & + \left( \frac{1}{32} ε^4 + \frac{1}{96} ε^6 \right) \sin(4 λ) \notag \\ & − \frac{1}{192} ε^6 \sin(6 λ) + O(ε^8) \label{eq:l2aeq} \\ α_\text{sun} & = λ_\text{sun} + S ≈ λ_\text{sun} + A_2 \sin(2 λ_\text{sun}) + A_4 \sin(4 λ_\text{sun}) + A_6 \sin(6 λ_\text{sun}) \label{eq:lambda2alpha} \end{align}
and the relationship between the declination \( δ_\text{sun} \) and the ecliptic longitude by
\begin{align} \arcsin(\sin(λ) \sin(ε)) & = (ε − \frac{1}{6} ε^3 + \frac{1}{120} ε^5) \sin(λ) \notag \\ & + \left( \frac{1}{6} ε^3 − \frac{1}{12} ε^5 \right) \sin(3 λ) \notag \\ & + \frac{3}{40} ε^5 \sin(5 λ) + O(ε^7) \\ δ_\text{sun} & ≈ D_1 \sin(λ_\text{sun}) + D_3 \sin^3(λ_\text{sun}) + D_5 \sin^5(λ_\text{sun}) \label{eq:lambda2delta} \end{align}
with \( A_2 \), \( A_4 \), \( A_6 \), \( D_1 \), \( D_3 \), and \( D_5 \) (measured in degrees) from the following table. The columns \( E_A \) and \( E_D \) show the greatest errors you make for \( α_\text{sun} \) and \( δ_\text{sun} \) when you use these approximations.
Table 4: Planets: Approximations for Right Ascension and Declination
\(A_2\) | \(A_4\) | \(A_6\) | \(E_A\) | \(D_1\) | \(D_3\) | \(D_5\) | \(E_D\) | |
---|---|---|---|---|---|---|---|---|
Mercury | −0.0000 | 0.0000 | 0.0000 | 0.0000 | ||||
Venus | −0.0305 | 0.0001 | 2.6427 | 0.0009 | 0.0036 | |||
Earth | −2.4680 | 0.0530 | −0.0014 | 0.0003 | 22.8008 | 0.5999 | 0.0493 | 0.0003 |
Mars | −2.8605 | 0.0712 | −0.0022 | 0.0004 | 24.3870 | 0.7331 | 0.0706 | 0.0011 |
Jupiter | −0.0424 | 0.0001 | 3.1151 | 0.0015 | 0.0034 | |||
Saturn | −3.2364 | 0.0911 | −0.0031 | 0.0009 | 25.7790 | 0.8649 | 0.0951 | 0.0010 |
Uranus | −42.5725 | 12.8039 | −2.6057 | 17.6902 | 56.9067 | −0.8355 | 26.1482 | 3.34 |
Neptune | −3.5195 | 0.1077 | −0.0039 | 0.0163 | 26.7577 | 0.9662 | 0.1164 | 0.060 |
Pluto | −17.1633 | 2.4178 | −0.3035 | 0.5052 | 48.3114 | 4.7880 | 4.3582 | 0.19 |
The approximations for Uranus aren't very good, because Uranus lies almost on its side: You'll do better to use the full formulas for Uranus.
As seen from Earth, and using equation \ref{eq:alpha-zon}, we find:
\[ α_\text{sun} = \arctan(\sin(12.0321°) × \cos(23.45°), \cos(12.0321°)) = 11.0639° \]
and with equation \ref{eq:delta-zon} we find
\[ δ_\text{sun} = \arcsin(\sin(12.0321°) × \sin(23.45°)) = 4.7585° \]
As seen from Mars, we find
\[ α_\text{sun} = \arctan(\sin(13.0435°) × \cos(25.19°), \cos(13.0435°)) = 11.8398° \]
and
\[ δ_\text{sun} = \arcsin(\sin(13.0435°) × \sin(25.19°)) = 5.5123° \]
Using equations \ref{eq:lambda2alpha} and \ref{eq:lambda2delta} we find for the Earth
\begin{align*} α_\text{sun} & = 12.0321° − 2.468° × \sin(2 × 12.0321°) \\ & + 0.053° × \sin(4 × 12.0321°) \\ & − 0.0014° × \sin(6 × 12.0321°) = 11.0639° \\ δ_\text{sun} & = 22.8008° × \sin(12.0321°) \\ & + 0.5999° × \sin(12.0321°)^3 \\ & + 0.0493° × \sin(12.0321°)^5 = 4.7585° \end{align*}
and for Mars
\begin{align*} α_\text{mars} & = 13.0435° − 2.8605° × \sin(2 × 13.0435°) \\ & + 0.0712° × \sin(4 × 13.0435°) \\ & − 0.0022° × \sin(6 × 13.0435°) = 11.8397° \\ δ_\text{mars} & = 24.387° × \sin(13.0435°) \\ & + 0.7331° × \sin(13.0435°)^3 \\ & + 0.0706° × \sin(13.0435°)^5 = 5.5124° \end{align*}
so the approximations yield practically the same results in these cases as the full equations \ref{eq:alpha-zon} and \ref{eq:delta-zon}.
Where a celestial body is in your sky depends on your geographical coordinates (latitude \( φ \) [phi] north, longitude \( l_\text{w} \) west), on the position of the body between the stars (its equatorial coordinates \( α \) and \( δ \)), and on the rotation angle of the planet at your location, relative to the stars. That latter angle is expressed in the sidereal time \( θ \) (theta). The sidereal time is the right ascension that is on the celestial meridian at that moment. The sidereal time is equal to
\begin{equation} \label{eq:sterrentijd} θ = (θ_0 + θ_1 × (J − J_{2000}) − l_\text{w}) \bmod{360°} \end{equation}
with \( θ_0 \) and \( θ_1 \) from the next table.
Table 5: Planets: Sidereal Time
\(θ_0\) | \(θ_1\) | |
---|---|---|
Mercury | 13.5964 | 6.1385025 |
Venus | 215.2995 | −1.4813688 |
Earth | 280.1600 | 360.9856235 |
Mars | 313.4803 | 350.89198226 |
Jupiter | 146.0727 | 870.5366420 |
Saturn | 174.3479 | 810.7939024 |
Uranus | 17.9705 | −501.1600928 |
Neptune | 52.3996 | 536.3128492 |
Pluto | 56.3183 | −56.3623195 |
For the Netherlands on Earth we find \( θ_\text{earth} = (560529.8477° − (−5°)) \bmod 360° = 14.8477° \) and for Gusev on Mars \( θ_\text{mars} = (544897.8367° − 184.6°) \bmod 360° = 33.2367° \).
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. To calculate the horizontal coordinates from the equatorial coordinates, you can use the following formulas:
\begin{align} \sin A \cos h & = \sin H \cos δ \\ \cos A \cos h & = \cos H \sin φ \cos δ − \sin δ \cos φ \\ \sin h & = \sin φ \sin δ + \cos φ \cos δ \cos H \label{eq:h} \\ H & = θ − α \label{eq:H} \\ A & = \arctan(\sin H, \cos H \sin φ − \tan δ \cos φ) \label{eq:A} \end{align}
The \( H \) is the hour angle, which indicates how long ago (measured in sidereal time) the celestial body passed through the celestial meridian.
For the Netherlands on Earth we find
\begin{align*} H & = 3.7838° \\ A & = \arctan(\sin(3.7838°), \cos(3.7838°) × \sin(52°) \\ & − \tan(4.7585°) × \cos(52°)) = 5.1302° \\ h & = \arcsin(\sin(52°) × \sin(4.7585°) \\ & + \cos(52°) × \cos(4.7585°) × \cos(3.7838°)) = 42.6542° \end{align*}
For Gusev on Mars we find
\begin{align*} H & = 21.3969° \\ A & = \arctan(\sin(21.3969°), \cos(21.3969°) × \sin(−14.6°) \\ & − \tan(5.5123°) × \cos(−14.6°)) = 131.9648° \\ h & = \arcsin(\sin(−14.6°) × \sin(5.5123°) \\ & + \cos(−14.6°) × \cos(5.5123°) × \cos(21.3969°)) = 60.7657° \end{align*}
At 12:00 UTC on 1 April 2004, the Sun as seen from the Netherlands stands about 5° west of south at 43° above the horizon, and as seen from the Gusev crater on Mars the Sun then stands about 3° south by northwest at 61° above the horizon.
The transit of a celestial body is the moment at which the body passes through the celestial meridian. The transit of the Sun is noon, the middle of the day, at 12 hours solar time. The hour angle \( H \) of the Sun is then equal to 0. We have
\begin{equation} θ = α_\text{sun} \pmod{360°} \label{eq:middag} \end{equation}
Using this and previous equations, you can find the time of transit by making a guess for a value of \( J \) for which equation \ref{eq:middag} holds, calculating \( θ \) and \( α_\text{sun} \) for it, and then checking whether they satisfy equation \ref{eq:middag}. If they do not, then you have to adjust \( J \). All in all, it is a big search operation for the correct value of \( J \). Moreover, equation \ref{eq:middag} does not help you to understand which things are most important to the time of solar transit. You do get such insight if you approximate the solution by neglecting smaller terms. The answer may then not be as accurate as the correct solution, but does indicate clearly what the solution looks like, is usually easier to calculate, and provides an excellent starting guess in the search for the real value of \( J \).
With equations \ref{eq:m}, \ref{eq:m2lambda}, \ref{eq:c}, \ref{eq:lambda2alpha}, and \ref{eq:sterrentijd}, and after omitting smaller terms, we find
\begin{equation} J_\text{transit} = J_{2000} + J_0 + l_\text{w} \frac{J_3}{360°} + J_1 \sin M + J_2 \sin(2 L_\text{sun}) \pmod{J_3} \label{eq:Jdoor} \end{equation}
where
\begin{align} L_\text{sun} & = M + Π + 180° \\ J_0 & = (M₀ + Π + 180° − θ_0) \frac{J_3}{360°} \pmod{J_3} \\ J_1 & = C_1 \frac{J_3}{360°} \\ J_2 & = A_2 \frac{J_3}{360°} \\ J_3 & = \frac{360°}{θ_1 − M_1} \end{align}
If you already have \( λ_\text{Sun} \), then you can use it instead of \( L_\text{Sun} \) in equation \ref{eq:Jdoor}: that is a tiny bit more accurate. \( J_0 \) provides the date and time of a transit of the Sun. \( J_1 \) shows by how much the time of transit can vary because of the eccentricity \( e \) of the orbit. \( J_2 \) indicates how much the time of transit can change because of the obliquity \( ε \) of the ecliptic. \( J_3 \) is the average length of the solar day (from one transit to the next). All values are measured in Earth days of 24 hours.
Table 6: Planets: Noon
\(J_0\) | \(J_1\) | \(J_2\) | \(J_3\) | |
---|---|---|---|---|
Mercury | 45.3495 | 11.4556 | 175.9386 | |
Venus | 87.8650 | −0.2516 | 0.0099 | −116.7505 |
Earth | 0.0009 | 0.0053 | −0.0069 | 1.0000000 |
Mars | 0.9044 | 0.0305 | −0.0082 | 1.027491 |
Jupiter | 0.3345 | 0.0064 | 0.4135775 | |
Saturn | 0.0766 | 0.0078 | −0.0040 | 0.4440276 |
Uranus | 0.1027 | −0.0106 | 0.0849 | −0.7183165 |
Neptune | 0.3841 | 0.0019 | −0.0066 | 0.6712575 |
Pluto | 3.8479 | −0.5023 | 0.3045 | −6.386797 |
To find the date and time of a solar transit near Julian Date \( J \), you now proceed as follows:
\begin{equation} n_{×} = \frac{J − J_{2000} − J_0}{J_3} − \frac{l_\text{w}}{360°} \end{equation}
and then take for \( n \) the whole number nearest to \( n_{×} \).
\begin{equation} J_{×} = J + J_3 × (n - n_{x}) \label{eq:schat-doorgang} \end{equation}
This \( J_{×} \) is a reasonable estimate for the date and time of the transit near \( J \), except that the \( J_1 \) and \( J_2 \) corrections are not in it yet.
\begin{equation} J_\text{transit} ≈ J_{×} + J_1 \sin M + J_2 \sin(2 L_\text{sun}) \end{equation}
If you want greater precision, then calculate the hour angle \( H \) of the Sun for \( J_\text{transit} \) and take
\begin{equation} J_\text{transit} − \frac{H}{360°} × J_3 \end{equation}
as improved value of \( J_\text{transit} \). You can repeat this until \( J_\text{transit} \) no longer changes.
For our example, we looked near \( J = 2453097 \). Which solar transit is closest to that in the Netherlands and in Gusev crater on Mars? For the Netherlands (\( l_\text{w} = −5° \)) we find
\[ n_{×} = \frac{2453097 − 2451545 − 0.0009}{1} − \frac{−5°}{360°} = 1552.0130 \]
so \( n = 1552 \), so
\[ J_{×} = 2453097 + 1 × (1552 − 1552.0130) = 2453096.9870 \]
For \( J_{x} \) we find \( M = 87.1679° \) and then
\[ L_\text{sun} = 87.1679° + 102.9372° + 180° = 370.1051° = 10.1051° \]
With that, we find
\begin{align*} J_\text{transit} & = 2453096.9869 + 0.0053 × \sin(87.1679°) \\ & − 0.0069 × \sin(2 × 10.1051°) = 2453096.9899 \end{align*}
The repetition method provides increased accuracy. For JD 2453096.9899 we find \( H = 0.1546° \), so a better estimate is \( J_\text{transit} = 2453096.9899 − \frac{0.1546°}{360°} × 1 = 2453096.9895 \). For that moment we find \( H \) equal to 0 to four digits after the decimal point, which is sufficiently precise. The solar transit at 5° east longitude happens on 1 April 2004 around 11:45 UTC.
For the Gusev crater on Mars (\( l_\text{w} = 184.6° \)) we find
\[ n_{×} = \frac{2453097 − 2451545 − 0.9044}{1.02749} − \frac{184.6°}{360°} = 1509.0840 \]
so \( n = 1509 \), so
\[ J_{×} = 2451545 + 1.02749 × (1509 − 1509.0840) = 2453096.9137 \]
For \( J_{x} \) we find \( M = 112.6079° \) and
\[ L_\text{sun} = 112.6079° + 70.9812° + 180° = 363.5891° = 3.5891° \]
With that, we find
\begin{align*} J_\text{transit} & = 2453096.9137 + 0.0305 × \sin(112.6079°) \\ & − 0.0082 × \sin(2 × 3.5891°) = 2453096.9408 \end{align*}
Now we use the repetition method. For JD 2453096.9408 we find \( H = 0.6556° \), so a better estimate is \( J_\text{transit} = 2453096.9408 − \frac{0.6556°}{360°} × 1.02749 = 2453096.9389 \). For that moment we find \( H = −0.0002° \) which yields no further change to \( J_\text{transit} \) (to four digits after the decimal point) so we are done. The solar transit in Gusev crater happens on 1 April 2004 around 10:32 UTC.
Above, we calculated when the Sun is highest in the sky, expressed in the Earthly time scale of the Julian Date. With this, you still don't know at what time that is on your clock, and that time scale is definitely not very handy if you are on another planet where the day has a different length from the day on Earth that the Julian Date is based on.
If you are not interested in transforming to times and dates for other planets, then you just want to know at what time the Sun is highest in the sky on your planet according to your clock, which is tied to your solar time and the season on your planet. Solar time is the time determined by the Sun. Three different kinds of solar time are relevant:
The answer to the question at what time the Sun is highest in the sky depends a lot on the time scale that you use. In true solar time, the answer is "always at exactly 12 o'clock", but true solar time is otherwise not very convenient at all. You can't use it if the Sun does not shine or if it is hidden behind clouds, and a mechanical clock or watch would have to be a lot more complicated if it had to show true solar time. In mean solar time, the answer is "on average at 12 o'clock", but it may differ from 12 o'clock from day to day. In official clock time, the answer is "on average at a fixed time", but that fixed time depends on your location, though it is usually near 12 o'clock.
So how large is the difference between the mean solar time and the true solar time? This is called the Equation of Time. The Equation of Time is how much you have to add to true solar time ("sundial time") to get mean solar time ("local time").
We found earlier (equation \ref{eq:middag}) for the moment at which the Sun is highest in the sky (transit) that \( θ = α_\text{Sun} \pmod{360°} \). For the right ascension \( α_\text{Sun} \) of the Sun we found \( α_\text{Sun} = L + C + S \). For the sidereal time \( θ \) measured in degrees we construct a different formula:
\begin{equation} θ = (L + t + 180°) \bmod 360° \end{equation}
Here \( t \) is the mean solar time measured in degrees (0° = midnight, 180° = noon). The explanation for this equation is as follows: The difference between the sidereal time and the solar time reflects the motion of the planet around the Sun so it increases by 360° in a planet year. The sidereal time is tied to the regular rotation of the planet around its axis, which does not notice the varying speed at which the planet orbits around the Sun (associated with eccentricity \( e \)) or the position of the Sun above or below the celestial equatior (associated with the obliquity \( ε \) of the ecliptic). Therefore, the sidereal time is tied to the mean longitude \( L \) of the Sun, which increases at a fixed rate by 360° during a planet year, and which is independent of \( e \) and \( ε \). During the descending equinox (when \( λ_\text{Sun} = L = 180° \)) the sidereal time is equal to the mean solar time, so we need the extra 180° in the formula. During a planet day, \( θ \) increases by 360° (the planet turns once around its axis) plus or minus a small amount because the sidereal time runs a little faster or slower than the mean solar time, but that difference is already caught in \( L \).
If we now set \( θ \) equal to \( α_\text{Sun} \) then we find
\begin{equation} L + t + 180° = θ = α_\text{Sun} = L + C + S \bmod 360° \end{equation}
or
\begin{equation} t = 180° + C + S \bmod 24 \end{equation}
The Equation of Time \( ∆t \) is equal to the difference between \( t \) and 180° (because 180° corresponds to 12:00:00 local mean solar time):
\begin{equation} ∆t = C + S \end{equation}
If you divide the Equation of Time measured in degrees by 15 then you get the Equation of Time measured in hours.
\( C \) depends on eccentricity \( e \) of the planet's orbit and on the mean anomaly \( M \) of the planet, i.e., on where the planet is compared to its perihelion. \( S \) depends on the obliquity \( ε \) of the ecliptic of the planet and on the longitude \( λ_\text{Sun} \) of the Sun, i.e., on the season. A first-order approximation is:
\begin{equation} ∆t ≈ C_1 \sin M + A_2 \sin(2 λ_\text{Sun}) \end{equation}
where we've omitted many smaller terms. The following table shows the amplitudes, i.e., the largest contribution that the \( C \) and \( S \) terms give to the Equation of Time for each planet, measured in planet minutes.
Table 7: Planets: Contributions to the Equation of Time
\(C\) | \(S\) | |
---|---|---|
Mercury | 94.5 | 0 |
Venus | 3.1 | 0.1 |
Earth | 7.7 | 9.9 |
Mars | 42.8 | 11.4 |
Jupiter | 22.2 | 0.2 |
Saturn | 25.4 | 13.0 |
Uranus | 21.2 | 178.1 |
Neptune | 4.1 | 14.1 |
Pluto | 114.6 | 69.3 |
For the Earth, the contributions of the orbit (\( C \)) and the season (\( S \)) are about equally large; for Uranus and Neptune the influence of the season is much greater than that of the orbit; and for the other planets the influence of the orbit is much greater than the influence of the seasons.
For the hour angle that corresponds to \( h = 0 \) we find
\begin{equation} H = \arccos(−\tan δ \tan φ) \label{eq:H0} \end{equation}
This equation has two solutions: if \( H \) is a solution, then \( −H \) is a solution, too. The solution with \( H \gt 0 \) is associated with sunset, and the solution with \( H \lt 0 \) is associated with sunrise. We write \( H_+ \) for the solution that has \( H ≥ 0 \). By using formulas found earlier and by omitting small terms, we find
\begin{align} H_+ & ≈ 90° + H_1 \sin λ_\text{sun} \tan φ \notag \\ & + H_3 \sin^3 λ_\text{sun} \tan φ × (3 + \tan^2 φ) \notag \\ & + H_5 \sin^5 λ_\text{sun} \tan φ × (15 + 10 \tan^2 φ + 3 \tan^4 φ) \label{eq:Happ} \end{align}
with
\begin{align} H_1 & = ε − \frac{1}{6} ε^3 + \frac{1}{120} ε^5 \\ H_3 & = \frac{1}{6} ε^3 − \frac{1}{12} ε^5 \\ H_5 & = \frac{1}{40} ε^5 \end{align}
where you must measure \( ε \) in radians instead of in degrees. You transform from degrees to radians by multiplying the number of degrees by \( π/180 ≈ 0.017453292 \).
The values of \( H_1 \), \( H_3 \), and \( H_5 \) are listed for all planets in the following table, measured in degrees. For Uranus and Pluto it is doubtful that the approximations yield reasonable results.
Table 8: Planets: Daytime Length Coefficients
\(H_1\) | \(H_3\) | \(H_5\) | |
---|---|---|---|
Mercury | 0.018 | ||
Venus | 2.643 | 0.001 | |
Earth | 22.801 | 0.600 | 0.016 |
Mars | 24.387 | 0.733 | 0.024 |
Jupiter | 3.115 | 0.002 | |
Saturn | 25.779 | 0.865 | 0.032 |
Uranus | 56.907 | −0.836 | 8.716 |
Neptune | 26.758 | 0.966 | 0.039 |
Pluto | 48.311 | 4.788 | 1.453 |
Sunrise is the moment at which the top of the solar disk touches the horizon in the morning. Sunset is the same, but at night. To calculate the times of sunrise and sunset, you must take into account (besides the effects mentioned earlier) also the following things:
To compensate for these two effects, you can set \( h \) equal to the \( h_0 \) from the following table (measured in degrees), instead of 0°. The refraction has been added in to the value for the Earth, but not for the other planets, because they have either no atmosphere of any consequence, or else they have one that is so dense that you cannot see the Sun at all from the surface. The table also gives the average diameter \( d_\text{Sun} \) of the solar disk (in degrees).
Table 9: Planets: Solar Disk Diameter and Refraction
\(h_0\) | \(d_\text{Sun}\) | \(\sin(h_0)\) | |
---|---|---|---|
Mercury | −0.69 | 1.38 | −0.0120 |
Venus | −0.37 | 0.74 | −0.0064 |
Earth | −0.83 | 0.53 | −0.0146 |
Mars | −0.17 | 0.35 | −0.0031 |
Jupiter | −0.05 | 0.10 | −0.0009 |
Saturn | −0.03 | 0.06 | −0.0005 |
Uranus | −0.01 | 0.03 | −0.0002 |
Neptune | −0.01 | 0.02 | −0.0002 |
Pluto | −0.01 | 0.01 | −0.0001 |
With equation \ref{eq:h} you can then find the hour angle at which the Sun with that \( α_\text{sun}, δ_\text{sun} \) rises or sets:
\begin{equation} H_\text{t} = \arccos\left( \frac{\sin h_0 − \sin φ \sin δ}{\cos φ \cos δ} \right) \label{eq:hbig} \end{equation}
Here again, \( H = H_\text{t} \gt 0 \) holds for sunset and \( H = −H_\text{t} \lt 0 \) for sunrise. The estimated JD of sunrise and sunset are then
\begin{eqnarray} J_\text{rise} & = J_\text{transit} − \frac{H_\text{t}}{360°} J_3 \\ J_\text{set} & = J_\text{transit} + \frac{H_\text{t}}{360°} J_3 \end{eqnarray}
You can improve these estimates by repetition similar to how we did it for the transit. For \( J_\text{rise} \) or \( J_\text{set} \) calculate the hour angle \( H \) and declination \( δ_\text{sun} \) of the Sun, and then use equation \ref{eq:hbig} to calculate the hour angle \( ±H_\text{t} \) that the Sun at those equatorial coordinates should have to be rising or setting. Then replace the estimates by
\begin{eqnarray} J_\text{rise} − \frac{H + H_t}{360°} J_3 \\ J_\text{set} − \frac{H − H_t}{360°} J_3 \end{eqnarray}
until they don't change anymore, to the desired precision.
We found earlier that the transit of the Sun as seen from the Earth occurs at \( J_\text{transit} = 2453096.9895 \). The declination of the Sun is then \( δ_\text{Sun} = 4.7545° \). According to the above table, we have for the Earth \( \sin h_0 = −0.0146 \). With \( φ = 52° \) we then find from equation \ref{eq:hbig} that
\[ H_\text{t} = \arccos\left( \frac{−0.0146 − \sin(52°) \sin(4.7545°)}{\cos(52°) \cos(4.7545°)} \right) = 97.4841° \]
The first estimate for the preceding sunrise is then
\[ J_\text{rise} = 2453096.9895 − \frac{97.4841}{360} × 1 = 2453096.7187 \]
For that instant of time we calculate the position of the Sun and find \( H = −97.5043° \) and \( δ_\text{sun} = 4.6501° \), and hence \( H_\text{t} = 97.3483° \). Then the JD correction is \( −\frac{−97.5043 + 97.3483}{360} × 1 = 0.0004 \), so a better estimate is \( J_\text{rise} = 2453096.7187 + 0.0004 = 2453096.7191 \).
If we repeat this improvement procedure once more then the correction is equal to 0 (to 4 digits after the decimal point), so we're done.
Sunrise occurs around 05:15 UTC.
The first estimate for the following sunset is
\[ J_\text{set} = 2453096.9895 + \frac{97.4841}{360} × 1 = 2453097.2603 \]
For that instant of time we calculate the position of the Sun and find \( H = 97.5042° \) and \( δ_\text{sun} = 4.8587° \), and hence \( H_\text{t} = 97.6199° \). Then the JD correction is \( −\frac{97.5042 − 97.6199}{360} × 1 = 0.0003 \), so a better estimate is \( J_\text{set} = 2453097.2606 + 0.0003 = 2453097.2606 \).
If we repeat this improvement procedure once more then the correction is equal to 0 (to 4 digits after the decimal point), so we're done.
Sunset occurs on 1 April 2004 around 18:15 UTC.
We found earlier that the transit of the Sun as seen from Gusev crater on Mars occurs at \( J_\text{transit} = 2453096.9389 \). The declination of the Sun is then \( δ_\text{sun} = 5.5000° \). According to the above table, we have for Mars \( \sin h_0 = −0.0031 \). With \( φ = −14.6° \) we then find from equation \ref{eq:hbig} that
\[ H_\text{t} = \arccos\left( \frac{−0.0031 − \sin(−14.6°) \sin(5.5000°)}{\cos(−14.6°) \cos(5.5000°)} \right) = 88.7472° \]
The first estimate for the preceding sunrise is then
\[ J_\text{rise} = 2453096.9389 − \frac{88.7472}{360} × 1.02749 = 2453096.6856 \]
For that instant of time we calculate the position of the Sun and find \( H = −88.7690° \) and \( δ_\text{sun} = 5.4494° \), and hence \( H_\text{t} = 88.7605° \). Then the JD correction is \( −\frac{−88.7690 + 88.7605}{360} × 1.02749 = 0.0000 \), so that first estimate was already accurate enough.
Sunrise in Gusev crater happens on 1 April 2004 around 04:27 UTC.
The first estimate for the following sunset is
\[ J_\text{set} = 2453096.9389 + \frac{88.7472}{360} × 1.02749 = 2453097.1922 \]
For that instant of time we calculate the position of the Sun and find \( H = 88.7690° \) and \( δ_\text{sun} = 5.5507° \), and hence \( H_\text{t} = 88.7339° \). Then the JD correction is \( −\frac{88.7690 − 88.7339}{360} × 1.02749 = −0.0001 \), so a better estimate is \( J_\text{set} = 2453097.1922 − 0.0001 = 2453097.1921 \). Another repeat does not improve this result anymore.
Sunset in Gusev crater happens on 1 April 2004 around 16:37 UTC.
If you watch sunrise or sunset from great altitude, then there are two more corrections to be made:
The correction \( ∆H \) to equation \ref{eq:H0} that is necessary because of the large Sun and the atmosphere is approximately equal to
\begin{equation} ∆H ≈ −\frac{h_0}{\sqrt{\cos^2(φ) − \sin^2(δ_\text{sun})}} \end{equation}
For the Netherlands and Belgium (with \( φ \) about 50°) this is about 5 to 6 minutes.
To calculate the duration of sunrise or sunset you must calculate the moment that the top of the solar disk touches the horizon, as described above (with \( h_0 \)), and also the moment that the bottom of the solar disk touches the horizon. To calculate that last one, you should use \( h_0 + d_0 \) instead of \( h_0 \).
For places on Earth, we can combine various approximations into the following (with results measured in hours UTC):
\begin{align} t_\text{transit} & ≈ \text{12h00m} + \frac{l_w}{15°} + 24 × (J_0 + J_1 \sin M + J_2 \sin 2 L_\text{Sun}) \notag \\ & = \text{12h01m} + \frac{l_w}{15°} + \text{7.6m} \sin M − \text{9.9m} \sin 2 L_\text{Sun} \\ t_\text{rise} & ≈ t_\text{transit} − \frac{H}{15°} \notag \\ & ≈ \text{6h00m} + \frac{l_w}{15°} + 24 × (J_0 + J_1 \sin M + J_2 \sin 2 L_\text{Sun}) \notag \\ & − \frac{H_1 \tan φ \sin L_\text{Sun} + H_3 \tan φ (3 + \tan^2 φ) \sin^3 L_\text{Sun}}{15°} \notag \\ & = \text{6h01m} + \frac{l_w}{15°} + \text{7.6m} \sin M − \text{9.9m} \sin 2 L_\text{Sun} \notag \\ & − \left( \text{1h31m} \tan φ \sin L_\text{Sun} + \text{2.2m} \tan φ (3 + \tan^2 φ) \sin^3 L_\text{Sun} + \frac{∆H}{15°} \right) \\ t_\text{set} & ≈ t_\text{transit} + \frac{H}{15°} \notag \\ & ≈ \text{18h00m} + \frac{l_w}{15°} + 24 × (J_0 + J_1 \sin M + J_2 \sin 2 L_\text{Sun}) \notag \\ & + \frac{H_1 \tan φ \sin L_\text{Sun} + H_3 \tan φ (3 + \tan^2 φ) \sin^3 L_\text{Sun}}{15°} \notag \\ & = \text{18h01m} + \frac{l_w}{15°} + \text{7.6m} \sin M − \text{9.9m} \sin 2 L_\text{Sun} \\ & + \left( \text{1h31m} \tan φ \sin L_\text{Sun} + \text{2.2m} \tan φ (3 + \tan^2 φ) \sin^3 L_\text{Sun} + \frac{∆H}{15°} \right) \end{align}
For a spot at \( φ = 50° \) and \( l_w = −5° \) (in the middle of the Netherlands) we find, approximately, measured in Central European Time (CET, equal to UTC plus one hour):
\begin{align} t_\text{rise} & ≈ \text{6h34m} − \text{1h48.6m} \sin L − \text{13.2m} \sin^3(L) + \text{7.6m} \sin M − \text{9.9m} \sin(2 L) \\ t_\text{transit} & ≈ \text{12h40m} + \text{7.6m} \sin M − \text{9.9m} \sin(2 L) \\ t_\text{set} & ≈ \text{18h46m} + \text{1h48.6m} \sin L + \text{13.2m} \sin^3(L) + \text{7.6m} \sin M − \text{9.9m} \sin(2 L) \end{align}
These formulas, though approximations, are reasonably accurate. A comparison, for all days of the year 2000, between the results of the full formulas and the results of the last approximate formulas shows differences of at most 3.5 minutes for the times of sunrise and sunset, and less than 1 minute for the time of transit.
From the preceding formulas we can deduce how the times of sunrise, transit, and sunset depend on the location (near any chosen spot). The times of sunrise, transit, and sunset (measured in a fixed time zone) get earlier by 4 minutes for each degree that you go to the east. This corresponds to \( \frac{2.16}{\cos φ} \) seconds per kilometer. The time of transit (i.e., high noon) doed not depend on the geographical latitude. In addition, sunrise gets earlier and sunset later by a number of minutes per degree that you go north that is approximately equal to \( \frac{1.59 \sin(L)}{\cos^2(φ)} \) and that is equivalent to about \( \frac{0.86 \sin(L)}{\cos^2(φ)} \) seconds per kilometer. The shift of the times of sunrise and sunset with latitude is greatest at the beginning of summer and winter. At the beginning of spring and autumn it is 0, so then the times of sunrise and sunset do not depend on the latitude (just like the time of transit).
The planetarium programs Redshift 3 and Redshift 5 and the current procedure (AA) yield the following values for the azimuth and altitude of the Sun in the sky at JD 2451545.0 (1−1-2000, 12:00 UTC) as seen from the location with latitude and longitude equal to zero on each of the planets:
Table 10: Solar Position From All Planets On 1−1-2000
Redshift 5 | Redshift 3 | AA | ||||
---|---|---|---|---|---|---|
\(A\) | \(h\) | \(A\) | \(h\) | \(A\) | \(h\) | |
Mercury | 89:59:36 | −4:28:58 | 90:00 | −4:21 | 270:00:00 | −4:29:32 |
Venus | 263:39:25 | −70:00:13 | 263:39 | −70:00 | 83:40:01 | −69:59:25 |
Earth | 178:03:33 | +66:57:05 | 177:37 | +66:57 | 357:20:14 | +66:55:48 |
Mars | 233:17:03 | +44:46:23 | 233:21 | +44:41 | 53:08:10 | +44:58:23 |
Jupiter | 275:36:38 | −56:47:58 | 273:19 | +22:27 | 93:19:59 | +23:07:15 |
Saturn | 115:08:60 | +33:21:18 | 115:08 | +33:18 | 294:58:39 | +32:59:52 |
Uranus | 223:07:12 | +45:26:37 | 223:09 | +45:25 | 44:13:28 | +45:46:26 |
Neptune | 217:28:20 | −54:09:28 | 218:03 | −54:32 | 37:50:23 | −54:44:57 |
Pluto | 125:32:00 | −43:59:37 | 125:32 | −44:00 | 305:34:32 | −43:59:11 |
And here are similar results for JD 2453097.0 (1−4-2004, 12:00 UTC):
Table 11: Solar Position From All Planets On 1−4-2004
Redshift 5 | Redshift 3 | AA | ||||
---|---|---|---|---|---|---|
\(A\) | \(h\) | \(A\) | \(h\) | \(A\) | \(h\) | |
Mercury | 89:52:29 | −87:19:22 | 89:01 | −87:06 | 270:00:00 | −87:19:36 |
Venus | 266:46:35 | +35:02:25 | 266:47 | +35:03 | 86:46:32 | +35:02:39 |
Earth | 11:08:57 | +85:07:29 | 11:10 | +85:07 | 194:18:16 | +85:05:20 |
Mars | 77:36:54 | −63:14:20 | 77:46 | −62:54 | 257:35:00 | −63:27:54 |
Jupiter | 92:10:30 | +46:10:15 | 91:36 | +20:10 | 271:36:17 | +20:04:19 |
Saturn | 231:05:23 | +47:32:41 | 231:24 | +47:11 | 50:41:50 | +47:56:45 |
Uranus | 143:35:09 | −72:11:32 | 144:19 | −72:22 | 321:08:05 | −72:43:02 |
Neptune | 173:45:46 | −61:31:02 | 172:54 | −61:57 | 353:01:20 | −61:59:07 |
Pluto | 135:36:03 | −41:02:35 | 135:37 | −41:04 | 315:36:48 | −41:01:38 |
If you take into account that the Redshift programs measure azimuth from the north point but I measure it from the south point (corresponding to a difference of 180°), then the results fit to within about a degree, except that the position from Jupiter as calculated by Redshift 5 is very different from that calculated by either Redshift 3 or myself, whereas the results from Redshift 3 and myself fit nicely. I therefore suspect that Redshift 5 is in error in that case.
http://aa.quae.nl/en/reken/zonpositie.html;
Last updated: 2016−04−25