Astronomy Answers: Position of the Sun

Astronomy Answers
Position of the Sun


[Astronomy Answers] [Dictionary] [AnswerBook] [Universe Family Tree] [Science] [Starry Sky] [Planet Positions] [Calculate] [Colophon]

1. Time ... 2. The Mean Anomaly ... 3. The Equation of Center ... 4. The Perihelion and the Obliquity of the Ecliptic ... 5. The Ecliptical Coordinates ... 6. The Equatorial coordinates ... 7. The Observer ... 8. Solar Transit ... 9. Equation of Time ... 10. Sunrise and Sunset ... 11. The Duration of Sunrise and Sunset ... 12. Earth ... 13. The Netherlands ... 14. Local dependence on location ... 15. Comparison with Planetarium Programs ... 16. Derivation of (Π) and (ε) ... 17. Derivation of ( θ_0 ) ... 18. Planet-based Topographic Coordinates ... 19. Planet-Bound Horizontal Coordinates ... 20. Coordinate Transformations Between All Systems ... 21. North Pole, South Pole, Primary Pole

\(\def\|{&}\DeclareMathOperator{\D}{\bigtriangleup\!} \DeclareMathOperator{\d}{\text{d}\!}\)

\( \DeclareMathOperator{\Rotx}{Rotx} \DeclareMathOperator{\Roty}{Roty} \DeclareMathOperator{\Rotz}{Rotz} \newcommand{\Matrix}[1]{\left( \begin{matrix} #1 \end{matrix} \right)} \)

The position of the Sun in the sky as seen from a planet (such as the Earth) is determined by four things:

  1. the time.
  2. the motion of the planet in its orbit around the Sun, which does not happen at constant speed, because of the eccentricity of the orbit of the planet.
  3. the angle between the axis of rotation of the planet and the plane of the orbit of the planet, which is not equal to 90 degrees. This causes the seasons.
  4. the location of the observer on the planet, which determines how high the Sun can get in the sky.

Below, I provide formulas that take all of these things into account. 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.

I updated the tables on this page in October 2016 so that they are based on the relevant planet characteristics published by the IAU in 2011, rather than those of 2000. The most important change is that Pluto is no longer regarded as a major planet but now as a dwarf planet, which means that now (since 2006) the other pole of Pluto is regarded as its north pole than before. See Chapter 21.

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.

1. Time

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.

2. The Mean Anomaly

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:

\({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° \).

3. The Equation of Center

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.

\({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° \).

4. The Perihelion and the Obliquity of the Ecliptic

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.

\({Π}\) \({ε}\)
Mercury 230.3265 0.0351
Venus 73.7576 2.6376
Earth 102.9373 23.4393
Mars 71.0041 25.1918
Jupiter 237.1015 3.1189
Saturn 99.4587 26.7285
Uranus 5.4634 82.2298
Neptune 182.2100 27.8477
Pluto 184.5484 119.6075

These values are different than those given on the NASA Fact Sheet (//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.

5. The Ecliptical Coordinates

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.0322° = 12.0322° \pmod{360°} \) and as seen from Mars \( λ_\text{sun} = 373.0664° = 13.0664° \pmod{360°} \).

6. The Equatorial coordinates

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.

\({A_2}\) \({A_4}\) \({A_6}\) \({E_A}\) \({D_1}\) \({D_3}\) \({D_5}\) \({E_D}\)
Mercury −0.0000 0.0000 0.0351 0.0000
Venus −0.0304 0.0001 2.6367 0.0009 0.0036
Earth −2.4657 0.0529 −0.0014 0.0003 22.7908 0.5991 0.0492 0.0003
Mars −2.8608 0.0713 −0.0022 0.0004 24.3880 0.7332 0.0706 0.0011
Jupiter −0.0425 0.0001 3.1173 0.0015 0.0034
Saturn −3.2338 0.0909 −0.0031 0.0009 25.7696 0.8640 0.0949 0.0010
Uranus −42.5874 12.8117 −2.6077 17.6902 56.9083 −0.8433 26.1648 3.34
Neptune −3.5214 0.1078 −0.0039 0.0163 26.7643 0.9669 0.1166 0.060
Pluto −19.3248 3.0286 −0.4092 0.5052 49.8309 4.9707 5.5910 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.0322°) × \cos(23.4393°), \cos(12.0322°)) = 11.0649° \]

and with equation \ref{eq:delta-zon} we find

\[ δ_\text{sun} = \arcsin(\sin(12.0322°) × \sin(23.4393°)) = 4.7565° \]

As seen from Mars, we find

\[ α_\text{sun} = \arctan(\sin(13.0664°) × \cos(25.1918°), \cos(13.0664°)) = 11.8605° \]

and

\[ δ_\text{sun} = \arcsin(\sin(13.0664°) × \sin(25.1918°)) = 5.5222° \]

Using equations \ref{eq:lambda2alpha} and \ref{eq:lambda2delta} we find for the Earth

\begin{align*} α_\text{sun} \| = 12.0322° − 2.4657° × \sin(2 × 12.0322°) \\ \| + 0.0529° × \sin(4 × 12.0322°) \\ \| − 0.0014° × \sin(6 × 12.0322°) = 11.0649° \\ δ_\text{sun} \| = 22.7908° × \sin(12.0322°) \\ \| + 0.5991° × \sin(12.0322°)^3 \\ \| + 0.0492° × \sin(12.0322°)^5 = 4.7565° \end{align*}

and for Mars

\begin{align*} α_\text{sun} \| = 13.0664° − 2.8608° × \sin(2 × 13.0664°) \\ \| + 0.0713° × \sin(4 × 13.0664°) \\ \| − 0.0022° × \sin(6 × 13.0664°) = 11.8605° \\ δ_\text{sun} \| = 24.3880° × \sin(13.0664°) \\ \| + 0.7332° × \sin(13.0664°)^3 \\ \| + 0.0706° × \sin(13.0664°)^5 = 5.5221° \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}.

7. The Observer

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} \pmod{360°} \end{equation}

with \( θ_0 \) and \( θ_1 \) from the next table. The \( \pmod{360°} \) means that on both sides of the \( = \) you can add or subtract any multiple of 360°. For \( θ \) the custom is to choose the value between 0° and 360°.

\({θ_0}\) \({θ_1}\)
Mercury 132.3282 6.1385025
Venus 104.9067 −1.4813688
Earth 280.1470 360.9856235
Mars 313.3827 350.89198226
Jupiter 145.9722 870.5360000
Saturn 174.3508 810.7939024
Uranus 29.6474 −501.1600928
Neptune 52.4160 536.3128662
Pluto 122.2370 56.3625225

\( θ_0 \) is the sidereal time (in degrees) at longitude 0° at the instant defined by \( J_{2000} \). \( θ_1 \) is the rate of change of the sidereal time, in degrees per day.

For the Netherlands on Earth we find

\[ θ_\text{aarde} = 560529.8347° − (−5°) = 14.8347° \pmod{360°} \]

and for Gusev on Mars

\[ θ_\text{mars} = 544897.7392° − 184.6° = 33.1392° \pmod{360°} \]

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 δ \label{eq:horizontal} \\ \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 \| = 14.8347° − 11.0649° = 3.7698° \\ A \| = \arctan(\sin(3.7698°), \cos(3.7698°) × \sin(52°) \\ \| − \tan(4.7565°) × \cos(52°)) = 5.1111° \\ h \| = \arcsin(\sin(52°) × \sin(4.7565°) \\ \| + \cos(52°) × \cos(4.7565°) × \cos(3.7698°)) = 42.6530° \end{align*}

For Gusev on Mars we find

\begin{align*} H \| = 21.2786° \\ A \| = \arctan(\sin(21.2786°), \cos(21.2786°) × \sin(−14.6°) \\ \| − \tan(5.5222°) × \cos(−14.6°)) = 132.1463° \\ h \| = \arcsin(\sin(−14.6°) × \sin(5.5222°) \\ \| + \cos(−14.6°) × \cos(5.5222°) × \cos(21.2786°)) = 60.8439° \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.

8. Solar Transit

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_0 + Π + 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.

\({J_0}\) \({J_1}\) \({J_2}\) \({J_3}\)
Mercury 45.3497 11.4556 175.9386
Venus 52.1268 −0.2516 0.0099 −116.7505
Earth 0.0009 0.0053 −0.0068 1.0000000
Mars 0.9047 0.0305 −0.0082 1.027491
Jupiter 0.3345 0.0064 0.4135778
Saturn 0.0766 0.0078 −0.0040 0.4440276
Uranus 0.1260 −0.0106 0.0850 −0.7183165
Neptune 0.3841 0.0019 −0.0066 0.6712575
Pluto 4.5635 −0.5024 0.3429 6.387672

To find the date and time of a solar transit near Julian Date \( J \), you now proceed as follows:

  1. Calculate

    \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_{×} \).

  2. Then calculate

    \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.

  3. Calculate \( M \) and \( L_\text{sun} \) for \( J_{×} \) and then get a better estimate for the date and time of the solar transit from

    \begin{equation} J_\text{transit} ≈ J_{×} + J_1 \sin M + J_2 \sin(2 L_\text{sun}) \end{equation}

  4. 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.

  5. \( J_\text{transit} \) is a Julian Date, which counts days and which is equal to a whole number at 12:00 UTC (i.e., noon, UTC). So, a Julian Date like 2453096.9898 that ends in .9898 is .0102 days before the next whole number, i.e., 0.0102 days = 0.0102×24 = 0.245 hours = 0.0102×24×60 = about 15 minutes before 12:00 UTC, i.e., about 11:45 UTC.

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.9870 + 0.0053 × \sin(87.1679°) \\ \| − 0.0068 × \sin(2 × 10.1051°) = 2453096.9766 \end{align*}

The repetition method provides increased accuracy. For JD 2453096.9766 we find \( H = −4.6663° \), so a better estimate is \( J_\text{transit} = 2453096.9766 − \frac{−4.6663°}{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.9047}{1.027491} − \frac{184.6°}{360°} = 1509.0822 \]

so \( n = 1509 \), so

\[ J_{×} = 2451545 + 1.027491 × (1509 − 1509.0822) = 2453096.9155 \]

For \( J_{x} \) we find \( M = 112.5978° \) and

\[ L_\text{sun} = 112.5978° + 70.9812° + 180° = 363.5790° = 3.5790° \]

With that, we find

\begin{align*} J_\text{transit} \| = 2453096.9155 + 0.0305 × \sin(112.5978°) \\ \| − 0.0082 × \sin(2 × 3.5891°) = 2453096.8945 \end{align*}

Now we use the repetition method. For JD 2453096.8945 we find \( H = −15.6787° \), so a better estimate is \( J_\text{transit} = 2453096.8945 − \frac{−15.6787°}{360°} × 1.027491 = 2453096.9393 \). That 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.

9. Equation of Time

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.

\({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.

10. Sunrise and Sunset

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.

\({H_1}\) \({H_3}\) \({H_5}\)
Mercury 0.035
Venus 2.636 0.001
Earth 22.137 0.599 0.016
Mars 23.576 0.733 0.024
Jupiter 3.116 0.002
Saturn 24.800 0.864 0.032
Uranus 28.680 −0.843 8.722
Neptune 26.668 0.967 0.039
Pluto 38.648 4.971 1.864

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:

  1. The Sun does not look like a point but like a disk, so when the center of the Sun sets (\( h = 0 \)), then half of the Sun is still above the horizon. The size of the Sun must be taken into account.
  2. The atmosphere of the planet (if it has one) bends light downward so that a celestial body appears slightly higher in the sky than it would have done without the atmosphere. This refraction of the light is strongest near the horizon, where the light has had to travel furthest through the air. On Earth, at the horizon as seen from sea level, the refraction amounts to some 0.57° on average. This means that a celestial body that with the atmosphere just set would have already been 0.57° below the horizon without the atmosphere.

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).

\({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.027491 = 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.027491 = 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.027491 = 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.027491 = −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:

  1. The horizon appears to be lower, by an amount in degrees that, on Earth, is approximately equal to the square root of the altitude in kilometers. For an altitude of 10 km (6 mi), the horizon is about 3.2° lower.
  2. There is more refraction of sunlight, because the sunlight almost touches the ground but then goes back up through the atmosphere to get to the observer. The total refraction is at most twice as great as it is at sea level, and I expect that the refraction from jet cruising altitude (10 km) above the Earth is almost at that maximum, so it is about 0.5° more than at sea level.

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.

11. The Duration of Sunrise and Sunset

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 \).

12. Earth

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} \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) \end{align}

13. The Netherlands

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.

14. Local dependence on location

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).

15. Comparison with Planetarium Programs

The planetarium programs Redshift 3, Redshift 5, Stellarium 0.15.0, the NASA web page JPL HORIZONS, 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:

Redshift 5 Redshift 3 AA Stellarium 0.15.0 JPL HORIZONS
\({A}\) \({h}\) \({A}\) \({h}\) \({A}\) \({h}\) \({A}\) \({h}\) \({A}\) \({h}\)
Mercury 89 59 36 −4 28 58 90 00 −4 21 270 01 32 −4 29 32 89 59 42.8 −4 22 03.8 90.0256 −4.4830
Venus 263 39 25 −70 00 13 263 39 −70 00 83 39 27 −69 59 26 92 14 33.2 −14 59 18.5 263.6545 −70.0006
Earth 178 03 33 +66 57 05 177 37 +66 57 357 18 25 +66 56 23 178 02 03.4 +66 57 37.5 178.0722 +66.9528
Mars 233 17 03 +44 46 23 233 21 +44 41 53 03 40 +45 04 15 233 20 09.0 +44 43 21.0 233.2109 +44.8716
Jupiter 275 36 38 −56 47 58 273 19 +22 27 93 19 50 +23 06 54 283 53 42.5 +77 07 12.8 273.3132 +22.3831
Saturn 115 08 60 +33 21 18 115 08 +33 18 294 58 03 +32 59 53 115 08 41.7 +33 21 41.8 115.1490 +33.3541
Uranus 223 07 12 +45 26 37 223 09 +45 25 44 13 02 +45 46 42 213 51 11.4 −51 47 54.9 223.1205 +45.4433
Neptune 217 28 20 −54 09 28 218 03 −54 32 37 49 37 −54 44 36 217 55 18.7 −53 54 11.7 217.4714 −54.1581
Pluto 125 32 00 −43 59 37 125 32 −44 00 125 35 22 −42 06 27 119 07 54.6 +30 48 30.8 305.5559 −42.1332

And here are similar results for JD 2453097.0 (1-4-2004, 12:00 UTC):

Redshift 5 Redshift 3 AA Stellarium 0.15.0 JPL HORIZONS
\({A}\) \({h}\) \({A}\) \({h}\) \({A}\) \({h}\) \({A}\) \({h}\) \({A}\) \({h}\)
Mercury 89 52 29 −87 19 22 89 01 −87 06 269 19 40 −87 19 34 89 56 12.4 −87 12 34.2 89.3290 −87.3182
Venus 266 46 35 +35 02 25 266 47 +35 03 86 46 41 +35 02 37 264 45 51.5 −59 43 08.3 266.7781 +35.0387
Earth 11 08 57 +85 07 29 11 10 +85 07 194 28 06 +85 05 14 11 13 57.3 +85 07 14.8 11.1353 +85.1259
Mars 77 36 54 −63 14 20 77 46 −62 54 257 30 37 −63 34 42 77 38 47.5 −63 07 28.1 77.5625 −63.3588
Jupiter 92 10 30 +46 10 15 91 36 +20 10 271 35 20 +19 04 53 96 19 20.0 +76 11 48.7 91.5977 +19.6703
Saturn 231 05 23 +47 32 41 231 24 +47 11 50 42 48 +47 56 56 231 16 56.8 +47 20 24.3 231.0880 +47.5457
Uranus 143 35 09 −72 11 32 144 19 −72 22 321 08 50 −72 43 10 109 18 45.9 +41 28 47.3 143.5871 −72.1924
Neptune 173 45 46 −61 31 02 172 54 −61 57 352 58 27 −61 58 26 174 30 03.9 −61 33 40.9 173.7614 −61.5171
Pluto 135 36 03 −41 02 35 135 37 −41 04 135 41 23 −39 00 02 127 19 52.2 +27 17 46.0 315.6817 −39.0309

If you take into account that I measure azimuth from the south but the other sources measure azimuth from the north (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 the other sources, and that the azimuth of the Sun from Pluto calculated by JPL HORIZONS differs by 180° from that of the other sources. I therefore suspect that Redshift 5 is in error for Jupiter. My azimuth for Pluto differs about 180° from that of JPL because I count azimuth from the south and JPL counts it from the north. The other sources are older than, or presumably haven't been updated for, the change in which of Pluto's poles is regarded as the north pole.

16. Derivation of \(Π\) and \(ε\)

Table 1 mentions \(Π\) en \(ε\). How do you calculate those?

The orientation of the orbit of a planet in space is expressed by three independent angles relative to a chosen fundamental plane and a chosen primary direction in that fundamental plane. These angles are part of the orbital elements. We use the inclination \(i\), the longitude \(Ω\) of the ascending node, and the argument \(ω\) of the perihelion as the three angles. The inclination is the angle between the orbit and the fundamental plane. The ascending node is one of the two intersections between the orbit and the fundamental plane. The longitude of the ascending node is the angle between the primary direction and the ascending node, as seen from the relevant focus of the orbit. Likewise, the argument of the perihelion is the angle between the ascending node and the perihelion of the orbit.

The calculation of the position of the Sun in the sky of a planet is easiest if the fundamental plane is chosen to be equal to the plane of the orbit of the planet, and the primary direction is chosen to be the intersection of the fundamental plane and the extension of the planet's equator into the sky. In the remainder of this text, we assume that that choice has been made. In the case of the Earth, the fundamental plane is the plane of the ecliptic, and the primary direction is the vernal equinox, and the coordinates are with respect to "the equator and equinox of the date".

Then the inclination of the orbit is zero, and the nodes aren't defined, and the ecliptic latitude of the Sun is zero. In that case the only angle needed to specify the orientation of the orbit is the angle between the primary direction and the perifocus of the orbit, measured in the plane of the orbit. That is the planet-based ecliptic longitude of the perifocus. We use symbol \(Π\) for that angle.

Once the planet-based ecliptic coordinates of the Sun are known, then we convert them into planet-based equatorial coordinates. For this we need to know the angle \(ε\) between the fundamental plane (chosen to be the plane of the orbit of the planet) and the equator of the planet. For the Earth, that angle is called the obliquity of the ecliptic.

So, to calculate the Sun's equatorial coordinates relative to the planet's equator, we need to know \(Π\) and \(ε\) of that planet.

How do we calculate \(ε\)? \(ε\) is the angle between the plane of the planet's orbit and the plane of the planet's equator. This is also the angle between the primary pole of the planet (perpendicular to the planet's equator) and the primary pole of the planet's orbit (perpendicular to the planet's orbit). The planet and its orbit each have two poles. One of those two has been chosen as the primary pole. Usually it is the north pole. See Chapter 21.

We calculate the direction of the planet's primary pole from its Earth-based equatorial coordinates, which are available on the internet. And we calculate the direction of the planet orbit's primary pole from its Earth-based orbital elements.

How do we calculate \(Π\)? \(Π\) is the angle between the planet's primary direction and the perifocus of the planet's orbit. We can calculate the direction of the perifocus of the planet in Earth-based coordinates from the widely available Earth-based orbital elements of the planet, so we need to calculate the Earth-based coordinates of the planet's primary direction. That primary direction is one of the intersections of the plane of the planet's orbit (the fundamental plane) and the plane of the extension of the planet's equator, so it is perpendicular to the primary pole of the planet and to the primary pole of the orbit of the planet.

Here are the mathematical details. We assume that the following quantities are given:

\({α_p}\) right ascension of north pole
\({δ_p}\) declination of north pole
\({Ω}\) longitude of ascending node
\({i}\) inclination of the orbit
\({ω}\) argument of perifocus
\({W_0}\) rotation angle at epoch
\({W_1}\) rotation speed at epoch
\({ε_0}\) obliquity of the ecliptic of Earth

All of these orbital elements except for \( W_0, W_1 \) are assumed specified with respect to the ecliptic (= orbit) or equator of the Earth, not those of the planet.

We use the following notation:

The following directions are of interest. They are defined in terms of lines and planes in space, which should be imagined extended out to infinity. Parallel lines in space then correspond to two opposite points on the celestial sphere, of which one is called the primary one. Parallel planes in space then correspond to a great circle on the celestial sphere. All lines perpendicular to those planes correspond to two opposite points on the celestial sphere, which are called the poles of the plane, of which one is called the primary pole of the planes. Intersecting planes correspond to two great circles on the celestial sphere, which intersect in two points, of which one is called the primary intersection,

The coordinate systems that are associated with the Earth have a lowercase letter for symbol, and the coordinate systems that are associated with the other planet have the corresponding uppercase letter for symbol. The coordinate systems are:

We derive formulas for converting between all of those coordinate systems and the Earth-based ecliptic coordinate system, which means we can convert between all of them. All of the described coordinate systems are right-handed.

The relevant quantities to begin with for all major planets and Pluto are listed in the following table. All angles are given in degrees. \( W_1 \) is given in degrees per Earth day (of 86400 seconds). \( α_p, δ_p, W_0, W_1 \) come from Report of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2009 and Erratum to: Reports of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2006 & 2009, found at //astrogeology.usgs.gov/groups/iau-wgccre. \( Ω, i, ω \) come from [Meeus] and are based on the VSOP87 model of Betagnon en Francou, except for Pluto. All values are given for the J2000.0 epoch (JD 2451545).

\({ α_p }\) \({ δ_p }\) \({ Ω }\) \({ i }\) \({ ω }\) \({ W_0 }\) \({ W_1 }\)
Mercury 281.0097 61.4143 48.330893 7.004986 29.125226 329.5469 6.1385025
Venus 272.76 67.16 76.679920 3.394662 54.883787 160.20 −1.4813688
Earth 0 90 174.873174 0 288.064174 190.147 360.9856235
Mars 317.68143 52.88650 49.558093 1.849726 286.502141 176.630 350.89198226
Jupiter 268.056595 64.495303 100.464441 1.303270 273.866868 284.95 870.5360000
Saturn 40.589 83.537 113.665524 2.488878 339.391263 38.90 810.7939024
Uranus 257.311 −15.175 74.005947 0.773196 98.999212 203.81 −501.1600928
Neptune 299.36 43.46 131.784057 1.769952 276.339634 253.18 536.3128492
Pluto 132.993 −6.163 110.307 17.140 113.768 302.695 56.3625225

For the rotation angle \( W \) of the planet at the time given by Julian Day Number \( J \) we have

\begin{equation} W = W_0 + W_1 × (J − J_{2000}) \end{equation}

As an example, we'll calculate \( Π, ε, θ_0 \) for Mars. I show all coordinates rounded to 7 digits after the decimal marker, and all angles to 4 digits after the decimal marker. For the obliquity of the ecliptic \( ε_0 \) I assume the value 23.4392911°.

  1. Fig. 1: Coordinate systems C, c, Q, q
    Fig. 1: Coordinate systems C, c, Q, q
    Figure 1 shows the Earth-based equatorial coordinate system \( q \) in fuchsia, the Earth-based ecliptic coordinate system \( c \) in red, the planet-based ecliptic coordinate system \( C \) in yellow, and the planet-based equatorial coordinate system \( Q \) in blue.

    The black circle represents the celestial sphere of infinite size, to which the Sun and stars and planets appear to be attached.

    All displayed ellipses are actually projections of circles on the celestial sphere, seen at an angle. Each circle represents the intersection of the celestial sphere and a flat plane of interest. The part of each circle on the near side of the celestial sphere is displayed as a solid curve, and the part on the far side of the celestial sphere is displayed as a dotted curve.

    All displayed straight lines from the center of the celestial sphere are projections of lines from the center of the sphere to a point on the surface of the sphere. Each of those lines represents the infinite extension of an axis or pole of interest.

    Actually, the circles and straight lines from the center represent entire families of parallel planes and lines in space. Because the celestial sphere is infinitely large, parallel lines in space seem to meet in a single point on the celestial sphere, just like parallel railroad tracks seem to meet at the horizon. Likewise, parallel planes in space seem to meet in a single great circle on the celestial sphere. In the current discussion we're not interested in differences in the origins (the zero points in space) of the coordinate systems, but only in their orientations in space. For the purposes of this discussion, we move the coordinate systems so that their origins coincide, without rotating them.

    The fuchsia curve represents the plane of the equator of the Earth. \( \vec{z}_{q} \) represents the primary pole of the Earth, which is the north pole. \( \vec{x}_q \) is the primary direction of the orbit of the Earth, one of the two intersections of the equator of the Earth and the orbit of the Earth. In the case of the Earth, it is called the vernal equinox. \( \vec{y}_q \) is perpendicular to both \( \vec{x}_q \) and \( \vec{z}_q \) and completes the Earth-based equatorial coordinate system.

    The red curve represents the plane of the orbit of the Earth, i.e., the ecliptic. \( \vec{z}_{c} \) represents the primary pole of the orbit of the Earth, i.e., the north pole of the ecliptic. \( \vec{y}_c \) completes the Earth-based ecliptic coordinate system. \( ε_0 \) is the angle over which to rotate the equatorial coordinate system \( q \) around its x axis to get the ecliptic coordinate system \( c \). In the case of the Earth, it is called the obliquity of the ecliptic.

    The yellow curve represents the plane of the orbit of the planet. \( \vec{z}_C \) represents the primary pole of the orbit of the planet. \( \vec{x}_C \) is the primary direction of the orbit of the planet, one of the two intersections of the equator of the planet and the orbit of the planet. \( \vec{y}_C \) completes the planet-based ecliptic coordinate system.

    The blue curve represents the plane of the equator of the planet. \( \vec{z}_Q \) represents the primary pole of the planet. \( \vec{x}_Q = \vec{x}_C \) is the primary direction, and \( \vec{y}_Q \) completes the planet-based equatorial coordinate system.

    \( ε_0 = ε_{cq} \) is the angle over which to rotate the Earth-based equatorial coordinate system \( q \) around its x axis to get the Earth-based ecliptic coordinate system \( c \). It is called the obliquity of the ecliptic.

    \( ε = ε_{CQ} \) is the angle over which to rotate the planet-based equatorial coordinate system \( Q \) around its x axis to get the planet-based ecliptic coordinate system \( C \). One might call it the planet-based obliquity of the ecliptic.

    Point \( Ω \) is the ascending node of the orbit of the planet (relative to the Earth's orbit). Angle \( Ω \) is the angle from the Earth-based primary direction (the vernal equinox) to the ascending node.

    \( i = ε_{Cc} \) is the inclination of the orbit of the planet (relative to the Earth's orbit). It is the angle over which to rotate the plane of the Earth's orbit around the ascending node to get the plane of the planet's orbit.

    Point \( P \) is the perifocus of the orbit of the planet: the point at which it is closest to the Sun.

    \( ω \) is the argument of the perifocus of the planet. It is the angle over which to rotate the ascending node around the z axis of the orbit of the planet to get \( P \).

    \( Π \) is the planet-based ecliptic longitude of the perifocus of the planet. It is the angle over which to rotate the planet-based primary direction around the z axis of the orbit of the planet to get \( P \).

  2. \( \vec{z}_{Qq} \) are the coordinates of the primary pole of the planet, expressed in Earth-based equatorial coordinates.

    \begin{equation} \vec{z}_{Qq} = \Matrix{ \cos α_p \cos δ_p \\ \sin α_p \cos δ_p \\ \sin δ_p } \end{equation}

    For Mars we find

    \[ \vec{z}_{Qq} = \Matrix{+0.4461587 \\ −0.4062376 \\ +0.7974418} \]

  3. Converting from Earth-based equatorial coordinates to Earth-based ecliptic coordinates goes via a rotation around the x axis by angle \( −ε_0 \):

    \begin{equation} C_{cq} = \Rotx(−ε_0) = \Matrix{ 1 \| 0 \| 0 \\ 0 \| \cos ε_0 \| \sin ε_0 \\ 0 \| −\sin ε_0 \| \cos ε_0 } \end{equation}

    and converting in the opposite direction, from Earth-based ecliptic coordinates to Earth-based equatorial coordinates goes via the opposite rotation:

    \begin{equation} C_{qc} = \Rotx(ε_0) = \Matrix{ 1 \| 0 \| 0 \\ 0 \| \cos ε_0 \| −\sin ε_0 \\ 0 \| \sin ε_0 \| \cos ε_0 } \end{equation}

    \( \vec{z}_{Qc} \) are the coordinates of the primary pole of the planet, expressed in Earth-based ecliptic coordinates.

    \begin{equation} \vec{z}_{Qc} = C_{cq} \vec{z}_{Qq} = \Matrix{ \cos α_p \cos δ_p \\ \sin ε_0 \sin δ_p + \cos ε_0 \sin α_p \cos δ_p \\ \cos ε_0 \sin δ_p − \sin ε_0 \sin α_p \cos δ_p } \end{equation}

    For Mars we find

    \[ \vec{z}_{Qc} = \Matrix{+0.4461587 \\ −0.0555116 \\ +0.8932306} \]

    which corresponds to Earth-based ecliptic longitude 352.9076° and latitude +63.2821°.

  4. \( \vec{z}_{Cc} \) are the coordinates of the primary pole of the planet orbit, expressed in Earth-based ecliptic coordinates.

    \begin{equation} \vec{z}_{Cc} = \Rotz(Ω) \Rotx(i) \Matrix{ 0 \\ 0 \\ 1 } = \Matrix{ \sin Ω \sin i \\ −\cos Ω \sin i \\ \cos i } \end{equation}

    where \( \Rotz \) represents rotation around the z axis, similar to \( \Rotx \) described before:

    \begin{equation} \Rotz(Ω) = \Matrix{ \cos Ω \| −\sin Ω \| 0 \\ \sin Ω \| \cos Ω \| 0 \\ 0 \| 0 \| 1 } \end{equation}

    For Mars we find

    \[ \vec{z}_{Cc} = \Matrix{+0.0245658 \\ −0.0209381 \\ +0.9994789} \]

    which corresponds to Earth-based ecliptic longitude 319.5581° and latitude +88.1503°.

  5. \( ε \) is the obliquity of the planet's orbit, which is the angle between the primary pole of the planet and the primary pole of the orbit of the planet. It has a value between 0 and 180°. Its cosine is equal to the inner product of the two poles:

    \begin{equation} \begin{split} \cos ε \| = \vec{z}_{Cc} ⋅ \vec{z}_{Qc} \\ \| = \sin i \left(−\cos Ω \left( \sin ε_0 \sin δ_p + \cos ε_0 \sin α_p \cos δ_p \right) + \sin Ω \cos α_p \cos δ_p\right) \\ \| + \cos i \left( \cos ε_0 \sin δ_p − \sin ε_0 \sin α_p \cos δ_p \right) \end{split} \end{equation}

    but going via the cosine leads to loss of precision for small values of \( ε \) (see the page about Distances in the Sky). For this, a better alternative is to go via

    \begin{equation} t^2 = \tan^2\left( \frac{1}{2}ε \right) = \frac{\left\lVert \vec{z}_{Cc} − \vec{z}_{Qc} \right\rVert^2}{\left\lVert \vec{z}_{Cc} + \vec{z}_{Qc} \right\rVert^2} \end{equation}

    If after that you are interested in \( \sin ε \) or \( \cos ε \), then they can be calculated from

    \begin{eqnarray} \sin ε \| = \frac{2t}{1 + t^2} \\ \cos ε \| = \frac{1 − t^2}{1 + t^2} \end{eqnarray}

    For Mars we find

    \begin{eqnarray*} \cos ε \| = \| \Matrix{+0.0245658 \\ −0.0209381 \\ +0.9994789} ⋅ \Matrix{+0.4461587 \\ −0.0555116 \\ +0.8932306} \\ \| = \| 0.0245658 × 0.4461587 + (−0.0209381) × (−0.0555116) + 0.9994789 × 0.8932306 = 0.9048877 \\ ε \| = \| \arccos(0.9048877) = 25.1918° \\ t^2 \| = \| \frac{\left\lVert \Matrix{+0.0245658 \\ −0.0209381 \\ +0.9994789} − \Matrix{+0.4461587 \\ −0.0555116 \\ +0.8932306} \right\rVert^2}{\left\lVert\Matrix{+0.0245658 \\ −0.0209381 \\ +0.9994789} + \Matrix{+0.4461587 \\ −0.0555116 \\ +0.8932306} \right\rVert^2} = \frac{\left\lVert \Matrix{−0.4215929 \\ +0.0345735 \\ +0.1062484} \right\rVert^2}{\left\lVert \Matrix{+0.4707245 \\ −0.0764497 \\ +1.8927095} \right\rVert^2} \\ \| = \| \frac{(−0.4215929)^2 + 0.0345735^2 + 0.1062484^2}{0.4707245^2 + (−0.0764497)^2 + 1.8927095^2} = \frac{0.1902247}{3.8097754} = 0.0499307 \\ ε \| = \| 2 \arctan(\sqrt{0.0499307}) = 25.1918° \\ \cos ε \| = \| \frac{1 − t^2}{1 + t^2} = \frac{1 − 0.0499307}{1 + 0.0499307} = 0.9048877 \\ \sin ε \| = \| \frac{2t}{1 + t^2} = \frac{2×\sqrt{0.0499307}}{1 + 0.0499307} = 0.4256504 \end{eqnarray*}

  6. \( \vec{x}_{Cc} \) are the coordinates of the primary direction of the planet, expressed in Earth-based ecliptic coordinates. The primary direction is perpendicular to the primary pole of the planet and also to the primary pole of the orbit of the planet. We calculate this using the vector cross product of the two poles. The cross product is perpendicular to both factors. The order of the two factors determines which of the two perpendicular directions is chosen.

    \begin{equation} \vec{x}_{Cc} = \frac{\vec{z}_{Qc} × \vec{z}_{Cc}}{\left\lVert \vec{z}_{Qc} × \vec{z}_{Cc} \right\rVert} = \frac{\vec{z}_{Qc} × \vec{z}_{Cc}}{|\sin ε|} \end{equation}

    The division by the length of the cross product or \( |\sin ε| \) is necessary because the length of this vector cross product isn't necessarily equal to 1.

    For Mars we find

    \[ \vec{x}_{Cc} = \frac{\Matrix{+0.4461587 \\ −0.0555116 \\ +0.8932306} × \Matrix{+0.0245658 \\ −0.0209381 \\ +0.9994789}}{0.4256504} = \frac{\Matrix{−0.0367801 \\ −0.4239833 \\ −0.0079780}}{0.4256504} = \Matrix{−0.0864092 \\ −0.9960834 \\ −0.0187432} \]

    which corresponds to Earth-based ecliptic longitude 265.0421° and latitude −1.0740°.

  7. Now that we have \( \vec{x}_{Cc} \) and \( \vec{z}_{Cc} \), we can complete the coordinate basis \( C_{Cc} \) for converting from Earth-based ecliptic coordinates to planet-based ecliptic coordinates. We combine the three vectors as rows of the matrix.

    \begin{eqnarray} \vec{y}_{Cc} \| = \vec{z}_{Cc} × \vec{x}_{Cc} \\ C_{Cc} \| = \Matrix{ \vec{x}'_{Cc} \\ \vec{y}'_{Cc} \\ \vec{z}'_{Cc} } \end{eqnarray}

    where the \( ' \) symbol represents the transposition operator which turns rows of a matrix into columns and columns into rows.

    The coordinate basis \( C_{cC} \) for converting from planet-based ecliptic coordinates to Earth-based ecliptic coordinates is the inverse of \( C_{Cc} \), which inverse is equal to

    \begin{equation} C_{cC} = C_{Cc}^\text{inv} = C_{Cc}' \end{equation}

    For Mars we find

    \[ \vec{y}_{Cc} = \Matrix{+0.0245658 \\ −0.0209381 \\ +0.9994789} × \Matrix{−0.0864092 \\ −0.9960834 \\ −0.0187432} = \Matrix{+0.9959568 \\ −0.0859037 \\ −0.0262788} \]

    which corresponds to Earth-based ecliptic longitude 355.0703° and latitude −1.5058°. Then

    \begin{eqnarray*} C_{Cc} \| = \| \Matrix{−0.0864092 \| −0.9960834 \| −0.0187432 \\ +0.9959568 \| −0.0859037 \| −0.0262788 \\ +0.0245658 \| −0.0209381 \| +0.9994789} \\ C_{cC} \| = \| \Matrix{−0.0864092 \| +0.9959568 \| +0.0245658 \\ −0.9960834 \| −0.0859037 \| −0.0209381 \\ −0.0187432 \| −0.0262788 \| +0.9994789} \end{eqnarray*}

  8. \( \vec{P}_c \) are the coordinates of the perifocus of the orbit of the planet, expressed in Earth-based ecliptic coordinates.

    \begin{equation} \vec{P}_c = \Rotz(Ω) \Rotx(i) \Rotz(ω) \Matrix{ 1 \\ 0 \\ 0 } = \Matrix{ \cos Ω \cos ω − \sin Ω \sin ω \cos i \\ \cos Ω \sin ω \cos i + \sin Ω \cos ω \\ \sin ω \sin i } \end{equation}

    For Mars we find

    \[ \vec{P}_c = \Matrix{+0.6486767 \| −0.7610641 \| 0 \\ +0.7610641 \| +0.6486767 \| 0 \\ 0 \| 0 \| +1} \Matrix{+1 \| 0 \| 0 \\ 0 \| +0.9994789 \| −0.0322782 \\ 0 \| +0.0322782 \| +0.9994789} \Matrix{+0.2840512 \| +0.9588091 \| 0 \\ −0.9588091 \| +0.2840512 \| 0 \\ 0 \| 0 \| +1} \Matrix{+1 \\ 0 \\ 0} = \Matrix{+0.9135923 \\ −0.4054519 \\ −0.0309486} \]

    which corresponds to Earth-based ecliptic longitude 336.0684° and latitude −1.7735°.

  9. \( \vec{P}_C \) are the coordinates of the perifocus of the orbit of the planet, expressed in planet-based ecliptic coordinates. We find it by converting \( \vec{P}_c \) from Earth-based ecliptic coordinates to planet-based ecliptic coordinates.

    \begin{equation} \vec{P}_C = C_{Cc} \vec{P}_c \end{equation}

    The perifocus is part of the orbit of the planet, so it lies in the planet-based plane of the ecliptic, so its planet-based ecliptic z coordinate and latitude are equal to 0.

    For Mars we find

    \[ \vec{P}_C = \Matrix{−0.0864092 \| −0.9960834 \| −0.0187432 \\ +0.9959568 \| −0.0859037 \| −0.0262788 \\ +0.0245658 \| −0.0209381 \| +0.9994789} \Matrix{+0.9135923 \\ −0.4054519 \\ −0.0309486} = \Matrix{+0.3255013 \\ +0.9455416 \\ 0} \]

    which corresponds to Earth-based ecliptic longitude 71.0041° and latitude 0°.

  10. \( Π \) is the planet-based ecliptic longitude of the perifocus.

    \begin{equation} Π = \arctan(P_{2C}, P_{1C}) \end{equation}

    For Mars we find

    \[ Π = \arctan(0.9455416, 0.3255013) = 71.0041° \]

  11. We can now finish the planet-based equatorial coordinate basis. We already had \( \vec{z}_{Qc} \), and the primary direction is shared between the planet-based equatorial coordinate system and the planet-based ecliptic coordinate system, so \( \vec{x}_{Qc} = \vec{x}_{Cc} \) and we already had \( \vec{x}_{Cc} \). We complete the coordinate basis \( C_{Qc} \).

    \begin{eqnarray} \vec{y}_{Qc} \| = \| \vec{z}_{Qc} × \vec{x}_{Qc} \\ C_{Qc} \| = \| \Matrix{ \vec{x}'_{Qc} \\ \vec{y}'_{Qc} \\ \vec{z}'_{Qc} } \\ C_{cQ} \| = \| C'_{Qc} \end{eqnarray}

    For Mars we find

    \begin{eqnarray*} \vec{y}_{Qc} \| = \| \Matrix{+0.4461587 \\ −0.0555116 \\ +0.8932306} × \Matrix{−0.0864092 \\ −0.9960834 \\ −0.0187432} = \Matrix{+0.8907726 \\ −0.0688209 \\ −0.4492080} \\ C_{Qc} \| = \| \Matrix{−0.0864092 \| −0.9960834 \| −0.01874317 \\ +0.8907726 \| −0.0688209 \| −0.4492080 \\ +0.4461587 \| −0.0555116 \| +0.8932306} \\ C_{cQ} \| = \| \Matrix{−0.0864092 \| +0.8907726 \| +0.4461587 \\ −0.9960834 \| −0.0688209 \| −0.0555116 \\ −0.0187432 \| −0.4492080 \| +0.8932306} \end{eqnarray*}

17. Derivation of \( θ_0 \)

Fig. 2: Coordinate systems C, E, Q, q, T
Fig. 2: Coordinate systems C, E, Q, q, T
Figure 2 shows the Earth-based equatorial coordinate system \( q \) in fuchsia, the planet-based equatorial coordinate system \( Q \) in blue, the planet-based ecliptic coordinate system \( C \) in yellow (all three as in Figure 1), the planet-based hybrid equatorial coordinate system \( E \) in olive, and the planet-based topographic coordinate system \( T \) in maroon.

Coordinate systems \( q \), \( Q \), and \( C \) were described in detail earlier.

\( \vec{x}_E \) represents the primary intersection of the equator of the Earth and the equator of the planet. This is the basis for the planet rotation angle \( W \) used by the IAU.

\( \vec{x}_T \) represents the primary location on the planet (with topographic latitude and longitude equal to zero). It rotates around the rotation axis of the planet (\( \vec{z}_Q = \vec{z}_E = \vec{z}_T \)) once per sidereal planet day. The planet's prime meridian is indicated by the dashed (front) and dotted (rear) ellipse passing through \( \vec{x}_T \) and \( \vec{z}_T \).

The sidereal time angle \( θ_0 \) says which planet-based right ascension passes through the zenith of an observer at the planet's prime location \( \vec{x}_T \) at the epoch (\( J_{2000} \)). This is equal to how far the planet's prime location \( \vec{x}_T \) has rotated beyond the planet-based primary direction \( \vec{x}_C = \vec{x}_Q \).

The angle \( υ = θ_0 − W_0 \) says how far the planet's primary intersection \( \vec{x}_E \) is beyond the planet's primary direction \( \vec{x}_C = \vec{x}_Q \).

How do you calculate \( θ_0 \)?

  1. \( \vec{x}_{Eq} \) are the coordinates of the primary intersection of the planet's equator with the Earth's equator, expressed in Earth-based equatorial coordinates, as defined by the IAU.

    \begin{equation} \vec{x}_{Eq} = \Matrix{ \cos(α_p + 90°) \cos 0 \\ \sin(α_p + 90°) \cos 0 \\ \sin 0 } = \Matrix{ −\sin α_p \\ \cos α_p \\ 0 } \end{equation}

    For Mars we find

    \[ \vec{x}_{Eq} = \Matrix{+0.6732522 \\ +0.7394129 \\ 0} \]

    which corresponds to Earth-based ecliptic longitude 47.6814° and latitude 0°.

  2. \( \vec{x}_{Ec} \) are the coordinates of the primary intersection expressed in Earth-based ecliptic coordinates.

    \begin{equation} \vec{x}_{Ec} = \Rotx(−ε_0) \vec{x}_{Eq} = \Matrix{ −\sin α_p \\ \cos ε_0 \cos α_p \\ −\sin ε_0 \cos α_p } \end{equation}

    For Mars we find

    \[ \vec{x}_{Ec} = \Matrix{+0.6732522 \\ +0.6783981 \\ −0.2941216} \]

    which corresponds to Earth-based ecliptic longitude 45.21813 and latitude −17.1049°.

  3. We now have the x axis \( \vec{x}_{Ec} \) of the planet-based hybrid equatorial coordinate system, in terms of Earth-based ecliptic coordinates. The z axis is the planet's primary pole

    \begin{equation} \vec{z}_{Ec} = \vec{z}_{Qc} \end{equation}

    Now that we have two of the three axes, we can complete the coordinate basis \( C_{Ec} \) for converting from Earth-based ecliptic coordinates to planet-based hybrid equatorial coordinates.

    \begin{eqnarray} \vec{y}_{Ec} \| = \vec{z}_{Ec} × \vec{x}_{Ec} \\ C_{Ec} \| = \Matrix{ \vec{x}'_{Ec} \\ \vec{y}'_{Ec} \\ \vec{z}'_{Ec} } \\ C_{cE} \| = C'_{Ec} \end{eqnarray}

    Normalization of the cross product isn't needed here because the two factors are perpendicular and of length 1.

    For Mars we find

    \begin{eqnarray*} \vec{y}_{Ec} \| = \| \Matrix{+0.4461587 \\ −0.0555116 \\ +0.8932306} × \Matrix{+0.6732522 \\ +0.6783981 \\ −0.2941216} = \Matrix{−0.5896388 \\ +0.7325944 \\ +0.3400465} \\ C_{Ec} \| = \| \Matrix{+0.6732522 \| +0.6783981 \| −0.2941216 \\ −0.5896388 \| +0.7325944 \| +0.3400465 \\ +0.4461587 \| −0.0555116 \| +0.8932306} \\ C_{cE} \| = \| \Matrix{+0.6732522 \| −0.5896388 \| +0.4461587 \\ +0.6783981 \| +0.7325944 \| −0.0555116 \\ −0.2941216 \| +0.3400465 \| +0.8932306} \end{eqnarray*}

  4. \( \vec{x}_{EQ} \) are the coordinates of the primary intersection of the equators of Earth and the planet, expressed in planet-based equatorial coordinates.

    \begin{equation} \vec{x}_{EQ} = C_{Qc} \vec{x}_{Ec} \end{equation}

    The z coordinate of that vector is equal to 0 by definition, because the primary intersection is on the celestial equator of the planet.

    For Mars we find

    \[ \vec{x}_{EQ} = \Matrix{−0.0864092 \| −0.9960834 \| −0.0187432 \\ +0.8907726 \| −0.0688209 \| −0.4492080 \\ +0.4461587 \| −0.0555116 \| +0.8932306} \Matrix{+0.6732522 \\ +0.6783981 \\ −0.2941216} = \Matrix{−0.7284035 \\ +0.6851484 \\ 0} \]

    which corresponds to Earth-based ecliptic longitude 136.7527° and latitude 0°.

  5. \( υ \) is the planet-based right ascension of the primary intersection \( \vec{x}_E \). This is the angle between the planet's primary direction and the primary intersection.

    \begin{equation} υ = \arctan(x_{2EQ}, x_{1EQ}) \end{equation}

    For Mars we find

    \[ υ = \arctan(0.6851484, −0.7284035) = 136.7527° \]

  6. Calculate \( θ_0 \) from \( W_0 \) and \( υ \):

    \begin{equation} θ_0 = W_0 + υ \end{equation}

    For Mars we find

    \[ θ_0 = 176.63° + 136.7527° = 313.3827° \]

18. Planet-based Topographic Coordinates

We derive the transformation of coordinates between planet-based topographic coordinate system \( T \) and Earth-based ecliptic coordinate system \( c \). Coordinate system \( T \) is the only one of the coordinate systems discussed on this page that rotates. The orientation of coordinate system \( T \) depends on the rotation angle \( W \).

\begin{equation} W = W_0 + W_1 × (J − J_{2000}) \end{equation}

The conversion from \( T \) to \( c \) goes as follows.

  1. \( \vec{z}_{Tc} \) are the coordinates of the primary pole of the planet, expressed in Earth-based ecliptic coordinates.

    \begin{equation} \vec{z}_{Tc} = \vec{z}_{Ec} \end{equation}

  2. \( \vec{x}_{Tc} \) are the coordinates of the primary location on the planet, expressed in Earth-based ecliptic coordinates. We find \( \vec{x}_{Tc} \) from \( \vec{x}_{Ec} \) by converting that to the planet-based equatorial coordinate system \( Q \), then rotating over angle \( W \) around the z axis of \( Q \) (the rotation axis of the planet), and then converting back to the Earth-based ecliptic coordinate system.

    \begin{equation} \vec{x}_{Tc} = C_{cQ} \Rotz(W) C_{Qc} \vec{x}_{Ec} = C_{cQ} \Matrix{ \cos W \| −\sin W \| 0 \\ \sin W \| \cos W \| 0 \\ 0 \| 0 \| 1 } C_{Qc} \vec{x}_{Ec} \end{equation}

  3. We find \( \vec{y}_{Tc} \) from \( \vec{x}_{Tc} \) and \( \vec{z}_{Tc} \), which are at right angles.

    \begin{equation} \vec{y}_{Tc} = \vec{z}_{Tc} × \vec{x}_{Tc} \end{equation}

  4. Now that we have all three coordinate axes we can construct the coordinate basis \( C_{Tc} \) for converting Earth-based ecliptic coordinates to planet-based topographic coordinates, and \( C_{cT} \) for the opposite conversion.

    \begin{eqnarray} C_{Tc} \| = \Matrix{ \vec{x}'_{Tc} \\ \vec{y}'_{Tc} \\ \vec{z}'_{Tc} } \\ C_{cT} \| = C'_{Tc} \end{eqnarray}

  5. \( C_{Tc} \) and \( C_{cT} \) vary with time, because they depend on \( W \) which varies with time.

For Mars for \( J = J_{2000} \) we find

\begin{eqnarray*} \vec{z}_{Tc} \| = \vec{z}_{Ec} = \vec{z}_{Qc} = \Matrix{ +0.4461587 \\ −0.0555116 \\ +0.8932306 } \\ \vec{x}_{Tc} \| = \Matrix{−0.0864092 \| +0.8907726 \| +0.4461587 \\ −0.9960834 \| −0.0688209 \| −0.0555116 \\ −0.0187432 \| −0.4492080 \| +0.8932306} \Matrix{−0.9982707 \| −0.0587837 \| 0 \\ +0.0587837 \| −0.9982707 \| 0 \\ 0 \| 0 \| +1} \Matrix{−0.0864092 \| −0.9960834 \| −0.01874317 \\ +0.8907726 \| −0.0688209 \| −0.4492080 \\ +0.4461587 \| −0.0555116 \| +0.8932306} \Matrix{+0.6732522 \\ +0.6783981 \\ −0.2941216} = \Matrix{−0.7067491 \\ −0.6341604 \\ +0.3136021} \\ \vec{y}_{Tc} \| = \Matrix{ +0.4461587 \\ −0.0555116 \\ +0.8932306 } × \Matrix{−0.7067491 \\ −0.6341604 \\ +0.3136021} = \Matrix{+0.5490429 \\ −0.7712063 \\ −0.3221690} \\ C_{Tc} \| = \Matrix{−0.7067491 \| −0.6341604 \| +0.3136021 \\ +0.5490429 \| −0.7712063 \| −0.3221690 \\ +0.4461587 \| −0.0555116 \| +0.8932306} \\ C_{cT} \| = \Matrix{−0.7067491 \| +0.5490429 \| +0.4461587 \\ −0.6341604 \| −0.7712063 \| −0.0555116 \\ +0.3136021 \| −0.3221690 \| +0.8932306} \end{eqnarray*}

19. Planet-Bound Horizontal Coordinates

From the planet-bound topographic coordinate system \( T \) we can deduce the planet-bound horizontal coordinate system \( H \) for location \( A \) at latitude \( Φ \) and longitude \( Λ \) on the planet, assuming the planet to be spherical. If the primary pole is the north pole, then \( Φ \) indicates north latitude and \( Λ \) indicates east longitude.

The z axis of \( H \) points to the zenith of an observer at that location. The y axis of \( H \) at that location points along the surface toward the primary pole. The x axis is perpendicular to the y and z axes and makes the coordinate system right-handed. If the north pole is the primary pole, then the x axis points to the east.

Coordinate system \( T \) is equal to coordinate system \( H \) of the primary location on the planet, except that the x axis of \( T \) is the z axis of \( H \), the z axis of \( T \) is the y axis of \( H \), and the y axis of \( T \) is the x axis of \( H \).

For location \( A \) on the planet, we find \( H \) from \( T \) as follows:

Combined this yields:

\begin{eqnarray} C_{HT} \| = \Matrix{ 0 \| 1 \| 0 \\ 0 \| 0 \| 1 \\ 1 \| 0 \| 0 } \Roty(Φ) \Rotz(−Λ) = \Matrix{ −\sin Λ \| \cos Λ \| 0 \\ −\cos Λ \sin Φ \| −\sin Λ \sin Φ \| \cos Φ \\ \cos Λ \cos Φ \| \sin Λ \cos Φ \| \sin Φ } \\ C_{TH} \| = C'_{HT} = \Matrix{ −\sin Λ \| −\cos Λ \sin Φ \| \cos Λ \cos Φ \\ \cos Λ \| −\sin Λ \sin Φ \| \sin Λ \cos Φ \\ 0 \| \cos Φ \| \sin Φ } \end{eqnarray}

Let's check. The direction from the center of the planet to location \( A \) at latitude \( Φ \) and longitude \( Λ \) is, expressed in topographic coordinates:

\begin{equation} \vec{r}_{AT} = \Matrix{ \cos Λ \cos Φ \\ \sin Λ \cos Φ \\ \sin Φ } \end{equation}

and that is indeed equal to the z axis (3rd column) of \( C_{TH} \), so equal to the direction of the zenith of point \( A \), expressed in topographic coordinates.

20. Coordinate Transformations Between All Systems

Fig. 3: Coordinate Transformations
Fig. 3: Coordinate Transformations
We now have a network of coordinate transformations with which we can convert coordinates from each of the coordinate systems \( q, c, Q, C, E, T \) to all other systems from that list. Figure 3 shows the coordinate systems. A line between two coordinate systems in that figure means we have described a direct transformation between the two systems, earlier on this page. By combining two direct transformations we can connect all other coordinate systems with each other. For example:

21. North Pole, South Pole, Primary Pole

The poles of a celestial object are the two intersections of the rotation axis of the object and the surface of the object. We call one of them the north pole and the other one the south pole, but which is which? Not everyone agrees about that. There are basically two ways to choose which pole is the north pole and which is the south pole:

  1. The north pole is the pole that points to the north side of a suitable fundamental plane.
  2. The north pole is the pole seen from above which the object appears to rotate counterclockwise around its axis.

The IAU uses the first of these definitions for the planets of the Solar System. The fundamental plane is then the Invariable Plane of the Solar System ― of which the poles are the average rotation axis for all orbital motion and rotation in the Solar System. The north side of that plane is the side to which the north pole of the Earth points. According to the IAU definition, the north poles of all planets of the Solar System point to the north side of the Invariable Plane of the Solar System.

For all other celestial objects, including dwarf planets and asteroids and comets in our Solar System, and for all celestial objects outside of the Solar System (including major planets), the IAU uses the second definition.

The great advantage of the second definition is that you need to look only at the celestial object itself, and do not need to search for some fundamental plane with which to compare the pole. The orientation of the fundamental plane might nog even be known very accurately.

For most planets in the Solar System both definitions yield the same result ― but not for Venus and Uranus. If you look down from above the IAU-defined north pole of Venus or Uranus, then those planets rotate clockwise around their axis.

Both definitions don't yield the same results for Pluto, either. Until 2006, Pluto was regarded as a small one of the major planets, to which the first definition applies, so the pole of Pluto that points to the north of the Invariable Plane of the Solar System was regarded as its north pole. Since 2006, Pluto is regarded as a dwarf planet, to which the second definition applies, so now the pole from above which Pluto appears to rotate in the counterclockwise direction is regarded as its north pole ― which is the opposite pole from what was regarded as its north pole until 2006.

Because it isn't always clear which pole is the north pole, I use the term "primary pole" instead in the above text, so that you cannot be led astray if your idea of which pole is the north pole were different from mine.



[AA]

languages: [en] [nl]

//aa.quae.nl/en/reken/zonpositie.html;
Last updated: 2021-07-19