Astronomy Answers: Phenomena in the Horizontal Coordinate System

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

\( \def\floorratio#1#2{\left\lfloor \frac{#1}{#2} \right\rfloor} \def\ceilratio#1#2{\left\lceil \frac{#1}{#2} \right\rceil} \)

As an example, we'll make our calculations for the Moon, for a location with north latitude \( φ \) = 52° and west longitude \( l_{\text{w}} \) = −5° (i.e., 5 degrees east longitude), based on the following positions (for 00:00 Central European Time at the beginning of the mentioned dates = 23:00 UTC of the preceding date):

	`α`	`δ`
2007-01-08	10h43m27s	+8°33′44″
2007-01-09	11h26m31s	+2°55′33″
2007-01-10	12h08m29s	−2°44′44″
2007-01-11	12h50m28s	−8°17′43″
2007-01-12	13h33m33s	−13°34′25″

To start with, we'll translate the right ascension \( α \) and declination \( δ \) to degrees (with 1 hour = 15°). For example, 10h43m27s = 10 hours + 43 minutes + 27 seconds = 10 + 43/60 + 27/3600 hours = 10.72417 hours = 10.72417 × 15° = 160.8625°, and +8°33′44″ = 8 degrees + 33 arc minutes + 44 arc seconds = 8 + 33/60 + 44/3600 degrees = 8.56222°. We also calculate the corresponding sidereal time \( θ \) (for example as explained on the Calculation Page about Sidereal Time), also measured in degrees.

	`α`	`δ`	`θ`
2007-01-08	160.8625	+8.5622	97.1266
2007-01-09	171.6292	+2.9258	98.1122
2007-01-10	182.1208	−2.7456	99.0979
2007-01-11	192.6167	−8.2953	100.0835
2007-01-12	203.3875	−13.5736	101.0692

1. Horizontal Coordinates

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 \cos δ \sin φ - \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. We measure it here in degrees, just like the right ascension \( α \) and the sidereal time \( θ \). Sometimes people measure them in hours instead of degrees, with 15° = 1 hour. If \( H \) is between 0° and 180°, then \( A \) must be between 0° and 180°, too.

To calculate equatorial coordinates from horizontal coordinates, you can use the following formulas:

\begin{align} \sin H \cos δ \| = \sin A \cos h \\ \cos H \cos δ \| = \cos A \cos h \sin φ + \sin h \cos φ \\ \sin δ \| = \sin h \sin φ - \cos h \cos φ \cos A \label{eq:δ} \\ α \| = θ - H \label{eq:αH} \\ H \| = \arctan(\sin A, \cos A \sin φ + \tan h \cos φ) \label{eq:H} \end{align}

The Calculation Page about Sidereal Time explains how you can calculate the sidereal time \( θ \) for a given instant, and the other way around.

For example, where is the Moon in the sky at 00:00 hours Central European Time at the beginning of 2007 January 9 (from \( φ = 52° \) and \( l_{\text{w}} = −5° \))? We found earlier, for that moment and location, \( θ = 98.1122° \). From table 1 we learn that \( α = 171.6292° \) and \( δ = +2.9258° \) at the desired instant. Formula \ref{eq:Hα} then yields

\[ H = 98.1122° - 171.6292° = −73.5169° \]

Formula \ref{eq:h} yields

\begin{align*} \sin h \| = \sin(2.9258°) \sin(52°) + \cos(2.9258°) \cos(52°) \cos(−73.5169°) \\ \| = 0.05104 × 0.78801 + 0.99870 × 0.61566 × 0.28373 = 0.21468 \end{align*}

\[ h = \arcsin(0.21468) = 12.397° \]

Formula \ref{eq:A} yields

\begin{align*} A \| = \arctan(\sin(−73.5169°), \cos(−73.5169°) \sin(52°) - \tan(2.9258°) \cos(52°)) \\ \| = \arctan(−0.95890, 0.19212) = −78.671° \end{align*}

So, we find that at 00:00 CET at the beginning of 9 January 2007 as seen from 52° north latitude and 5° east longitude, the Moon is 12.4° above the horizon, at 78.7° east of south.

Now we transform back from horizontal coordinates \( h \) en \( A \) to equatorial coordinates \( α \) and \( δ \). Formula \ref{eq:δ} yields

\begin{align*} \sin δ \| = \sin(12.397°) \sin(52°) - \cos(12.397°) \cos(52°) \cos(−78.671°) = 0.05104 \\ δ \| = \arcsin(0.05104) = 2.9258° \end{align*}

Equation \ref{eq:H} yields

\begin{align*} H \| = \arctan(\sin(−78.671°), \cos(−78.671°) \sin(52°) + \tan(12.397°) \cos(52°)) \\ \| = \arctan(−0.98051. 0.29013) = −73.5169° \end{align*}

and then equation \ref{eq:αH} yields

\[ α = 98.1122° - (−73.5169°) = 171.6292° \]

These values are equal to the values that we started with, within the precision of the calculations.

2. Transit

The transit is the moment when the object passes through the celestial meridian, which runs from the north point through the zenith to the south point. The object reaches its greatest altitude above the horizon then. The hour angle \( H \) is then equal to zero. If \( α \) is the object's right ascension, then we have at the moment \( t_\text{transit} \) when the object reaches its transit:

\begin{align} t \| = \frac{θ - θ_p - L_0 - L_1 ∆J_d}{θ_1} = \frac{θ}{θ_1} + t_d \pmod{t_\text{s}} \label{eq:lokaal} \\ t_d \| = -\frac{θ_p + L_0 + L_1 ∆J_d}{θ_1} \pmod{t_\text{s}} \\ θ_p \| ≈ θ_1 t_z - l_\text{w} \pmod{t_\text{s}} \\ L_0 \| = 99.967794687° \\ L_1 \| = 0.98564736628603° \\ θ_1 \| ≈ 15.04106864026192° \\ t_\text{s} \| = \frac{360°}{θ_1} ≈ 23.93446959 \end{align}

if we ignore that \( θ_p \), \( θ_1 \), hence also \( t_\text{s} \) depend slightly on the time. With equation \ref{eq:door} this leads to

\begin{equation} t_\text{transit} = \frac{α}{θ_1} + t_d \pmod{t_\text{s}} \label{eq:doortijd} \end{equation}

If \( α \) depends on the time (i.e, is not constant), then \( α \) in equation \ref{eq:doortijd} must be calculated for the instant defined by \( ∆J_d \) and \( t_\text{transit} \), so then we need to know \( t_\text{transit} \) before we can calculate \( t_\text{transit} \). This is certainly the case for the Sun, Moon, and planets. (For very accurate calculations, \( θ_p \), \( θ_1 \) and \( t_\text{s} \) also depend on the time, so then you must know the time before you can calculate the time, even if \( α \) is constant.)

If \( α \) changes only slowly with time, then you can pretend that it is constant. You'll get a wrong answer then, but the mistake will be small.

At what time \( t_\text{transit} \) (Central European Time, so \( t_\text{z} = −1 \)) does the Moon transit on 2007 January 9, seen from 52° north latitude and 5° east longitude? We fill in the values we already know into equation \ref{eq:doortijd} and find

\begin{equation} t_\text{transit} = \frac{α}{15.04106864} + 17.4115 \pmod{23.93446959} \label{eq:9jan}\end{equation}

Because \( α \) changes with time and we don't know the time \( t_\text{transit} \) yet, we now have the problem that was already described above. We can find estimates by assuming that \( α \) is constant. To also estimate how big of an error we then make, we'll make two such estimates, with the \( α \) from the beginning of January 9 and from the beginning of January 10.

With the \( α \) of January 9 (171.6292°) we find from equation \ref{eq:9jan}

\[ t_\text{transit} = 28.8222 \bmod 23.9345 = 4.8877 = \text{04:53} \]

and with the \( α \) of January 10 (182.1208°) we find from equation \ref{eq:9jan}

\[ t_\text{transit} = 29.5197 \bmod 23.9345 = 5.5852 = \text{05:35} \]

so the transit happens somewhere between 04:53 and 05:35 CET. By assuming \( α \) to be constant we can evidently make a mistake of about three quarters of an hour. This mistake is relatively large for the Moon because the Moon moves relatively fast between the stars in the sky. The mistake would be a lot smaller for the Sun or the planets.

Chapter 3 describes two ways in which you can handle cases where the sought quantity (here time) cannot be reduced to just one instance in the equation.

3. Methods for Searching

for a certain function \( f \). How can you find that \( x \) ― if it even exists (which is not the case for every function \( f \))?

As an example, we'll look for solutions to equation \ref{eq:xfx} for function \( f_1(x) ≡ \frac{x^2}{100} - 1 \).

3.1. The Repetition Method

If function \( f(x) \) near \( x = x_0 \) rises or falls more slowly than function \( x \) does (in other words: if \( |f'(x)| < 1 \)), then you can get ever closer to \( x_0 \) by starting near it and then keep putting the result back into function \( f \). So: \( x_2 = f(x_1) \), and then \( x_3 = f(x_2) \), and then \( x_4 = f(x_3) \), and so on until the \( x \) that you put in and the \( x \) that comes out are not noticeably different any more: then you have found a good estimate for \( x_0 \).

For function \( f_1 \) we choose 0 as the starting value for \( x \). We find (to four decimal places): \( f_1(0.0000) = −1.0000 \), \( f_1(−1.0000) = −0.9900 \), \( f_1(−0.9900) = −0.9902 \) and then \( f_1(−0.9902) = −0.9902 \), so we've found a solution already (to four decimal places).

If you make a graph of function \( f_1 \) then you find that there is a second solution, near \( x = 101 \), but the function turns out to rise too quickly there for the repetition method. If we start with that value, then we find: 101.01, 101.03, 101.07, 101.15, 101.32, 101.66, 102.34, 103.74, 106.62, 112.69, 125.99, 157.73, and ever faster to infinity, though we started out with much smaller steps than for the other solution.

If function \( f(x) \) rises or falls near \( x = x_0 \) more quickly than function \( x \), then you cannot apply the repetition method to function \( f \), but you can then apply that method to its inverse function \( f^\text{inv} \), if you can find it (which is not always possible). For a function \( f \) and its inverse \( f^\text{inv} \) we have

if you already have a reasonable estimate \( x \) for \( f \), then that same estimate is usually a good starting value for \( f^\text{inv} \).

If the function sometimes rises or falls fast and sometimes slowly, then it is important to choose a good starting value \( x \) that is sufficiently close to the desired value \( x_0 \).

There may be more than one value of \( x \) for which equation \ref{eq:xfx} holds. The solution that you find using the repetition method then depends on the starting value that you chose.

We would like to find the second solution for \( f_1 \), near \( x = 101 \). For \( f_1 \) for \( x > −1 \) we can find the inverse function: \( f_1^\text{inv}(x) = \sqrt{100 x + 100} \). For example, \( f_1(2) = −0.96 \) and \( f_1^\text{inv}(−0.96) = 2 \), and likewise for other values greater than −1. For \( f_1^\text{inv} \) with starting value 101, and using the repetition method, we then find, to four decimal places: 100.9950, 100.9926, 100.9914, 100.9908, and so on until 100.9902.

We have now found both solutions of \( f_1(x) ≡ \frac{x^2}{100} - 1 = x \), to four decimal places, namely \( x = −0.9902 \) en \( x = 100.9902 \).

3.2. The Restriction Method

then the solutions of equation \ref{eq:g} (the roots of function \( g \)) are the same as those of equation \ref{eq:xfx}.

We call \( g(x) \) the "deviation". We seek an \( x \) for which the deviation is equal to 0. The method works as follows:

Also for this method it is important to find good starting values, but it does not matter how fast the function rises or falls, as long as there are no singularities in the chosen interval.

If \( g(x) \) is a linear function (so that \( g \) rises or falls at a constant rate when \( x \) increases), then the first \( x_3 \) (from equation \ref{eq:schat}) is already the desired solution. This was not the case for the repetition method.

You can use formula \ref{eq:schat} even if \( x_1 \) and \( x_2 \) are on the same side of the desired \( x \), but then there is no guarantee that you will find a solution. If \( g_1 g_2 > 0 \), then replace with \( x_3 \) the \( x \) that has the greatest associated \( |g| \).

To our old function \( f_1(x) = \frac{x^2}{100} - 1 \) corresponds \( g(x) = \frac{x^2}{100} - x - 1 \). We choose as our starting values \( x_1 = −1 \) and \( x_2 = 0 \), and then find (to four decimal places) \( g_1 = 1.0400 \), \( g_2 = −1.0000 \). We use equation \ref{eq:schat} to find the middle value, and find \( x_3 = \frac{−1 × −2 - 1.0400 × 0}{−1 - 1.0400} = −0.9804 \) en \( g_3 = −0.0100 \).

We see that \( g_3 \) has the same sign as \( g_2 \), so we replace \( x_2 \) with \( x_3 \) and \( g_2 \) with \( g_3 \). We then have \( x_1 = −2 \) and \( x_2 = −0.9804 \).

We now repeat the procedure and find \( x_3 = −0.9901 \) and \( g_3 = −0.0001 \). After one more repetition we find \( x_3 = −0.9902 \) and \( g_3 = 0.0000 \), so we've found a solution, which turns out to be identical to the one found earlier with the repetition method.

We now try it with the starting values \( x_1 = 100 \) and \( x_2 = 102 \) and find

\({x_1}\)	\({g_1}\)	\({x_2}\)	\({g_2}\)	\({x_3}\)	\({g_3}\)
100	−1.0000	102	1.0400	100.9804	−0.0100
100.9804	−0.0100	102	1.0400	100.9901	−0.0001
100.9901	−0.0001	102	1.0400	100.9902	0.0000

so we've quickly found the second solution, at \( x = 100.9902 \), the same as we found with the repetition method.

3.3. Application to the Transit

Now we'll use these methods to calculate the transit of the Moon for which we found rather inaccurate estimates above.

First we try the repetition method. We had already calculated above that the transit time for the \( α \) of January 9 was 04:54. This is our first estimate.

Now we must calculate \( α \) for that instant. We assume that \( α \) changes linearly with time. The increase in \( α \) is 182.1208° − 171.6292° = 10.4917° in 24 hours, so the increase in 4:53 = 4.8877 hours is equal to 10.4917° × 4.8877 / 24 = 2.1367°, so at 04:53 we find \( α = 171.6292° + 2.1367° = 173.7658° \). With this \( α \) we find \( t_\text{transit} = \left( \frac{173.7658}{15.04106864} + 17.4115 \right) \bmod 23.9345 = 5.0298 = \text{05:02} \). If we repeat this once more then we find \( t_\text{transit} = 5.0339 \) en then \( t_\text{transit} = 5.0340 \), and after that it doesn't change anymore (to four decimal places), so we have found the answer: the Moon transits at 5.0340 uur = 05:02. Apparently the variation of \( α \) with time is slow enough that the repetition method leads to a solution. We assumed that \( α \) varies linearly with time, which is not entirely correct, so we still make a mistake, but that mistake is now a lot smaller than before.

Now we do it again, but using the restriction method. We measure the deviation as the difference between \( t_\text{transit} \) that we calculate from equation \ref{eq:doortijd} and the instant for which we calculated the \( α \), because when that difference is equal to 0 then we have found the right instant. For 00:00 hours CET on January 9 (i.e., \( x_1 = 0 \), measured in hours since that instant) we already found \( t_\text{transit} = 4.8877 \), so \( g_1 = 4.8877 - 0 = 4.8877 \). For 00:00 hours CET on January 10 (\( x_2 = 24 \) in the same time scale as \( x_1 \)) we already found \( t_\text{transit} = 5.5852 \), so \( g_2 = 5.5852 - 24 = −18.4148 \).

Equation \ref{eq:schat} then yields

\[ x_3 = \frac{−18.4148 × 0 - 4.8877 × 24}{−18.4148 - 4.8877} = 5.0340 \]

the same (within the calculation precision) as we found above using the repetition method, but here we find the right answer without any repetition (because we assumed ― just like for the repetition method ― that \( α \) varies linearly with time).

4. A Particular Altitude

When does the object reach a particular altitude \( h_0 \)? Suppose that we know the coordinates \( α \) and \( δ \) for the moment defined by \( ∆J_d \) and \( t \). Formula \ref{eq:h} then yields the hour angle \( H \):

\begin{align} q \| = \frac{\sin h_0 - \sin φ \sin δ}{\cos φ \cos δ} = \frac{\sin h_0}{\cos φ \cos δ} - \tan φ \tan δ \label{eq:q} \\ H \| = \arccos(q) \label{eq:H_h} \end{align}

This \( H \) can only be calculated if \( −1 ≤ q ≤ +1 \), i.e., if \( |q| ≤ 1 \). If \( |q| \) is greater than 1, then an object with that declination \( δ \) never reaches altitude \( h_0 \). If \( q > 1 \), then the object remains at greater altitudes than \( h_0 \), and if \( q < −1 \), then the object remains at smaller altitudes than \( h_0 \).

If we have an \( H \), then an estimate for the time \( t_\text{rise} \) at which the object passes upwards through altitude \( h_0 \) is

\begin{align} H_\text{rise} \| = -H = -\arccos(q) \\ θ_\text{rise} \| = α + H_\text{rise} = α - \arccos(q) \pmod{360°} \\ t_\text{rise} \| = t_d + \frac{α - \arccos(q)}{θ_1} = t_\text{transit} - \frac{\arccos(q)}{θ_1} \pmod{24} \end{align}

and for the time \( t_\text{set} \) at which the object passes downwards through altitude \( h_0 \)

\begin{align} H_\text{set} \| = +H = +\arccos(q) \\ θ_\text{set} \| = α + H_\text{set} = α + \arccos(q) \pmod{360°} \\ t_\text{set} \| = t_d + \frac{α + \arccos(q)}{θ_1} = t_\text{transit} + \frac{\arccos(q)}{θ_1} \pmod{24} \end{align}

These times are only correct if the equatorial coordinates \( α \) and \( δ \) do not vary with time. If they do (as they do for the Sun, Moon, and planets), then we can find the correct times using one of the search methods explained in chapter 3. You can then use as the deviation the difference between the calculate altitude (from equation \ref{eq:h}) and the desired altitude \( h_0 \), or else the difference between the times calculated using the above formulas and the time for which the equatorial coordinates \( α, δ \) were calculated.

When near 9 January 2007 does the Moon reach an altitude of 30°? Then \( h_0 = 30° \). With the coordinates \( α, δ \) for January 9 we find from equation \ref{eq:q}

\[ q = \frac{\sin(30°) - \sin(52°) × \sin(+2.9258°)}{\cos(52°) × \cos(+2.9258°)} = \frac{0.5 - 0.7880 × 0.0510}{0.6157 × 0.9987} = 0.7478 \]

This value is between −1 and +1 so we can really calculate a time.

Equation \ref{eq:H_h} yields \( H = \arccos(0.7478) = 41.6018° \). We remember from above that \( t_d = 17.4115 \) for this date and location, so we find

\[ t_\text{rise} = \left( 17.4115 + \frac{171.6292 - 41.6018}{15.04106864} \right) \bmod 24 = 26.0563 \bmod 24 = 2.1218 = \text{02:07} \]

and

\[ t_\text{set} = \left( 17.4115 + \frac{171.6292 + 41.6018}{15.04106864} \right) \bmod 24 = 31.5881 \bmod 24 = 7.6536 = \text{07:39} \]

so our first estimate is that the Moon passes upwards through 30° altitude around 02:07 and passes downwards through 30° again around 07:39.

With the coordinates of January 10 we find in a similar fashion

\begin{align*} q \| = 0.8744 \\ H \| = 29.0202° \\ t_\text{rise} \| = \left( 17.4115 + \frac{182.1208 - 29.0202}{15.04106864} \right) \bmod 24 = 3.6558 = \text{03:39} \\ t_\text{set} \| = \left( 17.4115 + \frac{182.1208 + 29.0202}{15.04106864} \right) \bmod 24 = 7.5146 = \text{07:31} \end{align*}

We now use the restriction method to find more accurate estimates. We begin with \( t_\text{rise} \). We have \( x_1 = 0 \) and \( x_2 = 24 \) hours, and \( g_1 = 2.1218 \) and \( g_2 = 3.6558 - 24 = −20.3442 \), and then find

\[ t_3 = \frac{−20.3442 × 0 - 2.1218 × 24}{−20.3442 - 2.1218} = 2.2667 = \text{02:16} \]

If we assume that the coordinate vary linearly with time, then we have now found the instant at which the Moon rises through the 30° altitude level. Similarly, we find for the time when the Moon sets through that level 7.6095 = 07:37.

5. A Particular Direction

When is the object in the direction defined by azimuth \( A_0 \)? We can calculate that as follows (according to Equations \ref{eq:D1}ff in Section 7):

\begin{align} s \| = \cos A_0 \cos δ \\ c \| = -\sin A_0 \cos δ \sin φ \\ a \| = -\sin A_0 \sin δ \cos φ \\ D \| = s^2 + c^2 - a^2 = \cos^2(δ) - \sin^2(A_0) \cos^2(φ) \\ p \| = a + c \\ q_1 \| = s + \sqrt{D} \\ q_2 \| = s - \sqrt{D} \\ f_1 \| = 2pq_1 \\ g_1 \| = p^2 - q_1^2 \\ f_2 \| = 2pq_2 \\ g_2 \| = p^2 - q_2^2 \\ H_1 \| = \arctan(f_1, g_1) \\ H_2 \| = \arctan(f_2, g_2) \end{align}

If \( D \lt 0 \) or if \( f = g = 0 \) then there are no solutions. If \( D = 0 \) and not \( f = g = 0 \) then there is a single solution. If \( D \gt 0 \) then there are two solutions, of which perhaps one is bad, if for that one \( f = g = 0 \) or \( f\sin(A_0) \lt 0 \).

If \( |δ| < |φ| \), then the object can stand in any direction, and stands in a given direction about once a day. If \( |δ| ≥ |φ| \), then the object cannot stand in all directions. The azimuth of the object must then be sufficiently close to 0° or 180°.

Figure 1 shows that in such a case any azimuth that can occur does occur about twice every day, once at greater altitude (in the upper half of the path) and once at lesser altitude (in the lower half of the path).

When on 9 January 2007 is the Moon exactly due east, seen from 52° north latitude and 5° east longitude? Then we must have \( A_0 = 270° \). The declination of the Moon is always less than \( φ = 52° \), so the Moon can stand at any azimuth as seen from that location.

On 9 January at 00:00 hours, the Moon was at \( δ = 2.9258° \). Then

\begin{align*} s \| = \cos(A_0)\cos(δ) = 0×0.9986965 = 0 \\ c \| = -\sin(A_0)\cos(δ)\sin(φ) = -(−1)×0.9986965×0.7880108 = 0.786984 \\ a \| = -\sin(A_0)\sin(δ)\cos(φ) = -(−1)×0.05104265×0.6156615 = 0.031425 \\ D \| = s^2 + c^2 - a^2 = 0.618356 \\ p \| = a + c = 0.818409 \\ q_1 \| = s + \sqrt{D} = 0 + \sqrt{0.818409} = 0.786356 \\ q_2 \| = s - \sqrt{D} = 0 - \sqrt{0.818409} = −0.786356 \\ f_1 \| = 2pq_1 = 2×0.818409×0.786356 = 1.28712 \\ g_1 \| = p^2 - q_1^2 = 0.818409^2 - 0.786356^2 = 0.051437 \\ f_2 \| = 2pq_2 = 2×0.818409×−0.786356 = −1.28712 \\ g_2 \| = p^2 - q_2^2 = 0.818409^2 - (−0.786356)^2 = 0.051437 \\ H_1 \| = \arctan(f_1, g_1) = 87.7115° \\ H_2 \| = \arctan(f_2, g_2) = −87.7115° \end{align*}

\( f_2\sin(A_0) \gt 0 \) so \( H = H_2 = −87.7115° \) is the correct solution. The right ascension of the Moon on 9 January at 00:00 hours was \( α = 171.6292° \), so the sidereal time \( θ \) at which the Moon (with those equatorial coordinates) would stand exactly in the east is

\[ θ = H + α = −87.7115° + 171.6292° = 83.9177° \]

and the corresponding clock time is

\[ t = \frac{θ}{θ_1} + t_d = \frac{83.9177}{15.04106864} + 17.4115 = 22.9907 = \text{22:59} \pmod{t_\text{s} = 23.93446959} \]

With the \( α \) and \( δ \) of the Moon for 10 January 00:00 hours we find

\( H = −92.1472° \), \( θ = 89.9736° \), \( t = 23.3934 = \text{23:23} \pmod{t_\text{s} = 23.93446959} \).

Now we use the restriction method. We take for \( x \) the time in hours since 9 January 00:00, so then \( x_1 = 0 \) and \( x_2 = 24 \). We take for \( g \) the difference between the estimated east time and the clock time on 9 January, in hours, so \( g_1 = 22.9907 - 0 = 22.9907 \) and \( g_2 = 23.3934 - 24 = −0.6066 \). Using equation \ref{eq:schat} we then find

\[ x_3 = \frac{−0.6066 × 0 - 22.9907 × 24}{−0.6066 - 22.9907} = 23.3830 = \text{23:22} \]

Our estimate is therefore that on 9 January 2007 as seen from 52° north latitude and 5° east longitude, the Moon stands due east at 23:22 CET.

6. All Moments

We saw earlier how to calculate an instant at which the object reaches a certain altitude or direction, given one or two initial estimates. If you want to know all instants at which the object reaches that altitude or direction within a specified period of time, then you cannot just take initial estimates with some randomly chosen fixed interval, because the objects move relative to the stars so the number of instants at which they show the desired behavior within one such interval won't always be the same. For example, the Moon and planets have one transit in most calendar days, but can sometimes have two or zero in one calendar day.

The solution to this problem is to tailor the initial estimates to the motion of the object, so that the collection of initial estimates on average stays in step with the motion of the object.

We saw before that to calculate the time of a horizontal phenomenon one always begins by calculating the time of the nearest transit. The phenomenon of interest cannot happen more than about 12 hours before or after the corresponding transit, because beyond that are the domains of the preceding or following transit. So, our first task is to find all transits during the time period.

We want to calculate \( H \) in such a way that its value is exactly 360° greater for each next transit, because then you can figure out how many transits there were between two instants from the \( H \) of those instants. How can we do that?

Suppose that you have a series of measurements \( x_i \) and a series \( m_i \) and you want to derive a series \( y_i \) such that

\begin{equation} m_i - \frac{1}{2}P ≤ y_{i+1} - y_i < m_i + \frac{1}{2}P \end{equation}

In other words: the expected increase of the \( y_i \) lies between \( m_i - \frac{1}{2}P \) and \( m_i + \frac{1}{2}P \). In practice, often all \( m_i \) are the same.

If you seek the smallest possible differences \( y_{i+1} - y_i \), then set \( m_i = 0 \): then the differences are between \( -\frac{1}{2}P \) and \( +\frac{1}{2}P \).

If you know that the \( y_i \) can only increase in value, then you get the smallest possible differences if \( m_i = \frac{1}{2}P \).

For example, assume that some interesting event happened at the following numbers of hours since the beginning of Monday: \( z_i \) = 21, 47, 73, 95, 117. Unfortunately, our assistant has written down the clock time but not the date, so the measurements that we got from him were: \( x_i \) = 21, 23, 1, 23, 21, which deviate by multiples of \( P = 24 \) (a whole number of days) from the \( z_i \), so we know the numbers \( \mod 24 \). We assume that our assistant wrote down the numbers in the order in which the events occurred. We can then calculate \( d_1 = x_2 - x_1 = 23 - 21 = 2 \) and likewise for the other values. The results are shown in Table 2.

Because these numbers are clock times, we know that the "real" difference between the instants was always greater than 0 (because clocks only run forwards). We try \( m_i = \frac{1}{2}P = 12 \), which means that we assume that the difference between successive "real" values is always between 0 and 24. We then find that \( e_1 = \left[ \frac{d_1 - m_1}{P} \right] = \left[ \frac{2 - 12}{24} \right] = 0 \) and so on. See the \( e_i \) column in the \( m_i = 12 \) part of table 2. Then we find \( s_1 = 0 \), \( s_2 = s_1 + e_1 = 0 + 0 = 0 \) and so on, and then \( y_1 = x_1 - P s_1 = 21 - 24×0 = 21 \) and 23, 25, 47, 69. The values that we find for \( y_i \) are indeed equal to the corresponding \( x_i \) except for multiples of 24, and indeed never decrease, but yet they are not all equal to the \( z_i \), because our assumption was incorrect that the difference between the successive real values was not more than 24 hours.

Our assistant tells us that there was always "about a day" between the events, so let's try again but now with \( m_i = 24 \), which means that we assume that the difference is between 12 and 36 hours (with 24 hours in the middle). The results are shown in table 2. Now we do find the correct values. Of course, in real life we won't know beforehand what the correct values are, otherwise we would not need to do any of these calculations in the first place, so it is important to find reliable information about which differences are possible.

The table also shows results for \( m_i = 0 \). These yield the smallest differences (positive or negative) that are possible \( \bmod 24 \).

Table 2: Examples of the Reconstruction of Mod-Values

		\({m_i = 0}\)			\({m_i = 12}\)			\({m_i = 24}\)
\({x_i}\)	\({d_i}\)	\({e_i}\)	\({s_i}\)	\({y_i}\)	\({e_i}\)	\({s_i}\)	\({y_i}\)	\({e_i}\)	\({s_i}\)	\({y_i}\)
21	2	0	0	21	0	0	21	−1	0	21
23	−22	−1	0	23	−1	0	23	−2	−1	47
1	22	1	−1	25	0	−1	25	0	−3	73
23	−2	0	0	23	−1	−1	47	−1	−3	95
21			0	21		−2	69		−4	117

6.1. Till How Far Does This Work?

If you have a good \( m_i \) for your \( x_i \) then you can find the right \( y_i \) even if there are many periods \( P \) between successive \( y_i \). However, this cannot be extended indefinitely.

Suppose that we know of \( y \) that it increases with time at a rate that is between \( μ - δ \) and \( μ + δ \) per unit time, and suppose that we know \( x = y \mod P \) for a number of moments spaced at a time interval of \( ∆t \). Then we know that the corresponding successive \( y \) grow at a rate of between \( (μ - δ)∆t \) and \( (μ + δ)∆t \) per measurement. We can then use the reconstruction method with \( m_i = μ∆t \), which means that we assume that the increase of \( y \) lies between \( μ∆t - \frac{1}{2}P \) and \( μ∆t + \frac{1}{2}P \). This assumption fails if \( δ∆t ≥ \frac{1}{2}P \), so if we set \( m_i = μ∆t \) then we must make sure that

This estimate is based on the assumption that the deviations with respect to \( μ \) are all in the same direction. If the deviations are periodic or random, then they'll partly cancel out, so then it will take longer until the total deviation grows too large. The above estimate is a lower limit: if you keep to that estimate then you'll not go wrong, but greater \( ∆t \) may be OK in practice.

For example, the right ascension of the Moon (relative to the equinox of the date) between the years 2000 and 4737 grows at a rate of between 10.37° and 17.39° per day, so \( μ = \frac{10.37 + 17.39}{2} = 13.88 \) per day and \( δ = \frac{17.39 - 10.37}{2} = 3.51 \) per day. If we want to apply the reconstruction method to those right ascensions, with \( P = 360 \), and if we set \( m_i = μ∆t \), then we must have \( ∆t < \frac{360}{2×3.51} = 51.3 \) days. If we have right ascensions for moments that are more than 51 days apart, then we should not expect that the reconstruction method will yield the correct reconstruction.

If we want keep things simple and set \( m_i = 0 \), then we must have \( ∆t < \frac{360}{17.39} = 20.7 \). In other words, if we set \( m_i = 0 \), then the time interval between the measurement or calculation of successive right ascensions may be up to 20 days.

6.2. Application to Transits

For the Sun, Moon, and planets the hour angle \( H \) always increases with time. The average increase of the hour angle per day is close to 360°. If you have one measurement or calculation of \( H \) after each interval of \( ∆t \) days, then the expected increase of \( H \) per measurement is \( 360° × ∆t \), so then you can use the reconstruction method with \( m_i = 360° × ∆t \).

Based on calculations for the years 2000 − 4737, if you set \( m_i = 360° × ∆t \), then \( ∆t \) must not be greater than is shown in the following table. The table also shows the increase \( ∆α \) of the right ascension in degrees per day, and the increase \( ∆H \) of \( H \) per day.

	\({∆α}\)	\({∆H}\)	\({∆t}\)
	min	max	min	max
Sun	0.88	1.11	359.89	360.11	2800
Mercury	−1.50	2.42	358.56	362.48	140
Venus	−0.70	1.38	359.60	361.69	200
Moon	10.37	17.39	343.58	350.60	21
Mars	−0.44	0.84	360.14	361.44	240
Jupiter	−0.15	0.26	360.72	361.15	300
Saturn	−0.09	0.14	360.83	361.08	320
Uranus	−0.05	0.07	360.94	361.04	340
Neptune	−0.03	0.04	360.97	361.04	350

So the problem is largest for the Moon, which was to be expected because the Moon moves fastest with respect to the stars. If we have at least one \( H \) per 21 days then we can use the reconstruction method (with \( m_i = 360° × ∆t \)) to figure out how many transits there were.

\begin{align} n \| = n_2 - n_1 + 1 \\ n_1 \| = ⌈H_1⌉ \\ n_2 \| = ⌊H_2⌋ \end{align}

where \( ⌊ \, ⌋ \) indicates rounding down to the nearest whole number, and \( ⌈ \, ⌉ \) indicates rounding up to the nearest whole number.

From the \( α \) en \( θ \) from table 1 we calculate for 00:00 hours CET \( H \mod 360° \) and then using the reconstruction method (\( P = 360° \), \( m_i = 360 \)) the "real" \( H \):

	\({α}\)	\({θ}\)	\({H \bmod 360}\)	\({H}\)
2007-01-08	160.8625	97.1266	296.2641	296.2641
2007-01-09	171.6292	98.1122	286.4831	646.4831
2007-01-10	182.1208	99.0979	276.9770	996.9770
2007-01-11	192.6167	100.0835	267.4668	1347.4668
2007-01-12	203.3875	101.0692	257.6817	1697.6817

We find \( n_1 = \ceilratio{H(t_1)}{360°} = 1 \) and \( n_2 = \floorratio{H(t_2)}{360°} = 4 \), and hence for the number of transits between \( t_1 \) = 8 January 00:00 CET and \( t_2 \) = 12 January 00:00 CET: \( n = 4 - 1 + 1 = 4 \),

You can now find first estimates for the times of transit in different ways. One way is to approximate the increase in \( H \) as being linear. If we use \( t_1 \), \( t_2 \), \( H_1 = H(t_1) \), and \( H_2 = H(t_2) \) for that, then we find

\begin{align} H(t) \| = a t + b \\ a \| = \frac{H_2 - H_1}{t_2 - t_1} \\ b \| = \frac{H_1 t_2 - H_2 t_1}{t_2 - t_1} \end{align}

The transits happen when \( H = 360° k \) with \( k \) a whole number. If we set \( H(t) = 360° k \) then we find

\begin{align} t \| = c k + d \\ c \| = \frac{360°}{a} = 360° \frac{t_2 - t_1}{H_2 - H_1} \\ d \| = -\frac{b}{a} = -(t_2 - t_1)\frac{H_1 t_2 - H_2 t_1}{H_2 - H_1} \end{align}

For our example we measure \( t \) in days since 0 January, 00:00 hours CET (so \( t_1 = 8 \) and \( t_2 = 12 \)) and we find

\begin{align*} a \| = \frac{1697.6817 - 296.2641}{12 - 8} = 350.3544 \\ b \| = \frac{296.2641 × 12 - 1697.6817 × 8}{12 - 8} = −2506.5711 \\ c \| = \frac{360}{350.3544} = 1.0275 \\ d \| = \frac{2506.5711}{350.3544} = 7.1544 \end{align*}

For our first estimate for the moments of transit, we've now found \( t = 1.0275 k + 7.1544 \). For \( k \) from \( n_1 = 1 \) through \( n_2 = 4 \) we then find \( t \) equal to 8.1819, 9.2094, 10.2370, 11.2645. That 8.1819 means 0.1819 days after the beginning of January 8, which is 24×0.1819 = 4.3660 hours after the beginning of January 8, i.e., 04:22 hours CET on January 8. Likewise, the other estimates are 05:02 CET on January 9, 05:41 CET on January 10, and 06:21 CET on January 11. We can improve these estimates using the methods described before.

If the time difference \( t_2 - t_1 \) is so large that you are afraid that a linear approximation would deviate too much from the real hour angle (if the time difference is at least as great as \( ∆t \) from table 3), then you can use the linear approximation to approximate the first transit (for \( k = 1 \), so \( t = c + d = \frac{360° - b}{a} \)), then use that approximation to find the true time of the transit, and then add \( c \) to that true time to find the first approximation for the next transit, and so on.

7. Derivation: From Azimuth to Hour Angle

7.1. Introduction

Given the equatorial coordinates of a celestial object and the geographical coordinates of the observer, what is the hour angle \( H \) at which that object reaches azimuth \( A \)? We cannot use equation \ref{eq:H} for that, because it requires that we already know the height \( h \) of the object above the horizon. Equation \ref{eq:A} contains only values that we know already, so that provides a better basis. That formula was

\begin{equation} A = \arctan(\sin H, \cos H \sin φ - \tan δ \cos φ) \label{eq:A1} \end{equation}

which means that \( \sin A \) and \( \sin H \) have the same sign, so \( \sin(A)\sin(H) ≥ 0 \). Equation \ref{eq:A1} leads to

\begin{equation} \tan A = \frac{\sin H}{\cos H \sin φ - \tan δ \cos φ} \label{eq:A2} \end{equation}

\begin{equation} \sin H - \tan A \sin φ \cos H = - \tan A \tan δ \cos φ \label{eq:AtoHprobleem}\end{equation}

7.2. General Linear Combination of Sine and Cosine

of which we want to know the solutions for \( x \), given \( s \), \( c \), and \( a \). We substitute

\begin{align} t \| = \tan\left( \frac{1}{2}x \right) \\ \cos x \| = \frac{1 - t^2}{1 + t^2} \\ \sin x \| = \frac{2t}{1 + t^2} \end{align}

\begin{equation} \frac{2st}{1 + t^2} + \frac{c(1 - t^2)}{1 + t^2} = a \end{equation}

However,\( t = \arctan\left( \frac{1}{2}x \right) \) can be infinitely large, which can be awkward. \( \sin(x) \) and \( \cos(x) \) never exceed 1, and from them we can find \( x \) using the two-argument arctan function.

\begin{align} x \| = \arctan(\sin x, \cos x) \\ \| = \arctan\left( \frac{2t}{1 + t^2}, \frac{1 - t^2}{1 + t^2} \right) \end{align}

Because \( 1 + t^2 \) is always positive, we can multiply the two arguments of the arctan function with it without changing its outcome.

\begin{align} x \| = \arctan(2t, 1 - t^2) \\ \| = \arctan\left( 2\frac{s ± \sqrt{D}}{a + c}, 1 - \left( \frac{s ± \sqrt{D}}{a + c} \right)^2 \right) \\ \| = \arctan\left( 2\frac{s ± \sqrt{D}}{a + c}, \frac{(a + c)^2 - \left( s ± \sqrt{D} \right)^2}{(a + c)^2} \right) \end{align}

Now the division by \( a + c \) can still be a problem, because it might be equal to 0 and then we get infinite results, or even 0/0 which has no solution. \( (a + c)^2 \) is never negative, so we can multiply the two arguments of the arctan function by it without changing its outcome.

\begin{equation} x = \arctan\left( 2(a + c)\left( s ± \sqrt{D} \right), (a + c)^2 - \left( s ± \sqrt{D} \right)^2 \right) \end{equation}

Now we're not troubled by infinities anymore, except if \( s \), \( c \), or \( a \) themselves can be infinitely large.

For example, suppose that we seek solutions \( x \) of

\[ \sin x + 3 \cos x = 2 \]

Then \( s = 1 \), \( c = 3 \), and \( a = 2 \), and then

\begin{align*} D \| = 1^2 + 3^2 - 2^2 = 6 \\ t \| = \frac{1 ± \sqrt{6}}{5} \\ x \| = 2 \arctan t \end{align*}

So \( t ≈ −0.289898 \) or \( t ≈ 0.689898 \), from which follow \( x ≈ 1.207828 \text{ rad} = 69.20343° \) or \( x ≈ 5.718859 \text{ rad} = 327.6665° \).

These values of \( x \) satisfy the given equation.

With the infinity avoiding version we get

\begin{align*} x \| = \arctan\left( 2×(2 + 3)×(1 ± \sqrt{6}), (2 + 3)^2 - (1 ± \sqrt{6})^2 \right) \\ \| = \arctan( 10 ± 10\sqrt{6}, 18 ∓ 2\sqrt{6} ) \end{align*}

\[ x = \arctan( 10 + 10\sqrt{6}, 18 - 2\sqrt{6} ) ≈ 1.207828 \text{ rad} = 69.20343° \]

\[ x = \arctan( 10 - 10\sqrt{6}, 18 + 2\sqrt{6} ) ≈ 5.718859 \text{ rad} = 327.6665° \]

just like before.

7.3. Application to the Problem

\begin{align} x \| = H \\ s \| = 1 \\ c \| = -\tan A \sin φ = -\frac{\sin A \sin φ}{\cos A} \\ a \| = -\tan A \cos φ \tan δ = -\frac{\sin A \cos φ \sin δ}{\cos A \cos δ} \end{align}

The divisions by \( \cos A \) and \( \cos δ \) are undesirable, because those can be equal to 0 and that gives problems. If we multiply \( s \), \( c \), and \( a \) all by the same (non-zero) factor, then the solutions \( x \) remain the same. If we multiply all of them by \( \cos A \cos δ \) then we're rid of the divisions. Then

\begin{align} x \| = H \\ s \| = \cos A \cos δ \\ c \| = -\sin A \cos δ \sin φ \\ a \| = -\sin A \sin δ \cos φ \end{align}

\begin{align} D \| = s^2 + c^2 - a^2 = \cos^2(δ) - \sin^2(A)\cos^2(φ) \label{eq:D1} \\ p \| = a + c \\ q_1 \| = s + \sqrt{D} \\ q_2 \| = s - \sqrt{D} \\ f_1 \| = 2pq_1 \\ g_1 \| = p^2 - q_1^2 \\ f_2 \| = 2pq_2 \\ g_2 \| = p^2 - q_2^2 \\ H_1 \| = \arctan( f_1, g_1 ) \\ H_2 \| = \arctan( f_2, g_2 ) \end{align}

If \( D \gt 0 \) then we find two solutions for \( H \). Do both of them satisfy equation \ref{eq:A1}? One way to find out is to enter each of them into \ref{eq:A1} and see if that leads to a true statement. But we can reject potential solutions earlier. The sign of \( \sin H \) must be equal to the sign of \( \sin A \), so if \( f \sin A \lt 0 \) then that \( H \) is wrong.

If \( |δ| = 90° \), and \( |φ| = 90° \) or \( \sin A = 0 \), then \( s = c = a = 0 \) so \( D = f = g = 0 \) and then \( H \) cannot be calculated ― but that is appropriate because then the object is always at the same location in the sky (in the zenith or the nadir) and not just for a particular value of \( H \).

If \( |φ| \lt |δ| = 90° \) and \( \sin A ≠ 0 \) then \( s = c = 0 \) and \( a ≠ 0 \) so \( D \lt 0 \) so then there is no solution. That is appropriate, because if \( |δ| = 90° \) then the object stands in a celestial pole, which always has \( \sin A = 0 \).

If \( |φ| ≤ |δ| \lt 90° \), then the object tracks circles around the nearest celestial pole, and then it cannot stand at every \( A \), but never gets further from that celestial pole in azimuth then a certain distance \( A_x ≤ 90° \)

\begin{equation} A_x = \arcsin\left( \frac{\cos^2 δ}{\cos^2 φ} \right) \end{equation}

For \( A \) at less than \( A_x \) from the azimuth of the nearest celestial pole, we have \( D \gt 0 \) and then there are two \( H \) that are both appropriate. For \( A \) at exactly \( A_x \) from the azimuth of the nearest celestial pole we have \( D = 0 \) and then there is one \( H \), which is then appropriate. For \( A \) at more than \( A_x \) from the azimuth of the nearest celestial pole we have \( D \lt 0 \) and then there are no solutions.

If \( |δ| \lt |φ| \lt 90° \), then the object orbits obliquely (and possibly very widely) around the zenith and the nadir and then \( A \) can reach any value. In that case always \( D \gt 0 \) so there are two solutions, and also \( f_1f_2 \lt 0 \) so only one of those two solutions has the same sign as \( \sin(A) \) and is therefore useful, and that is the solution for which \( f\sin(A) \gt 0 \), if \( \sin(A) ≠ 0 \).

\( \sin(A) = 0 \) is a special case: If \( \sin A = 0 \) (when the object is due north or due south) and \( |δ| \lt 90° \), then \( c = a = 0 \), and \( s = \cos(δ) \) (if \( A = 0° \)) or \( s = -\cos(δ) \) (if \( A = 180° \)), so \( D \gt 0 \) and \( f = 0 \), and \( g = 0 \) or \( g = −2\cos^2(δ) \lt 0 \). \( f = g = 0 \) yields no solution, so then only \( g \lt f = 0 \) remains, so \( H = \arctan(0,−1) = 180° \).

In summary, if \( D \lt 0 \) or \( f = g = 0 \), then there is no solution. If \( D = 0 \) and not \( f = g = 0 \) then there is one solution and that solution is appropriate. If \( D \gt 0 \) then there are two solutions of which only the one for which \( f\sin(A) \gt 0 \) or not \( f = g = 0 \) is appropriate.

Let's check. Let \( A = 133° \), \( φ = −64° \) en \( δ = 17° \). Then, according to equations \ref{eq:D1}ff:

\begin{align*} s \| = \cos(133°)\cos(17°) \\ \| = −0.6819984×0.9563048 = −0.6521983 \\ c \| = -\sin(133°)\cos(17°)\sin(−64°) \\ \| = −0.7313537×0.9563048×−0.898794 = 0.6286139 \\ a \| = -\sin(133°)\sin(17°)\cos(−64°) \\ \| = −0.7313537×0.2923717×0.4383711 = −0.09373564 \\ D \| = (−0.6521983)^2 + 0.6286139^2 - (−0.09373564)^2 = 0.8117316 \\ p \| = −0.09373564 + 0.6286139 = 0.5348782 \\ q_1 \| = −0.6521983 + \sqrt{0.8117316} = 0.2487632 \\ q_2 \| = −0.6521983 - \sqrt{0.8117316} = −1.5531598 \\ f_1 \| = 2×0.5348782×0.2487632 = 0.2661161 \\ g_1 \| = (−0.09373564)^2 - 0.2487632^2 = 0.2242116 \\ f_2 \| = 2×0.5348782×−1.5531598 = −1.661503 \\ g_2 \| = (−0.09373564)^2 - (−1.5531598)^2 = −2.126211 \\ H_1 \| = \arctan(0.2661161, 0.2242116) = 49.88475° \\ H_2 \| = \arctan(−1.661503, −2.126211) = −141.9946° \end{align*}

Are both solutions correct? \( \sin(A) \gt 0 \) so only solutions with \( f \gt 0 \) can be correct. That means that the solution \( H = −141.9946° \) is wrong, and only \( H = 49.88475° \) can be correct. Let's check by entering each one into Equation \ref{eq:A1}:

\begin{align*} A_1 \| = \arctan(\sin(−141.9946°), \cos(−141.9946°)\sin(−64°) - \tan(17°)\cos(−64°)) \\ \| = \arctan(−0.6157363, −0.7879523×−0.898794 - 0.3057307×0.4383711) \\ \| = \arctan(−0.6157363, 0.5741833) = −47° \\ A_2 \| = \arctan(\sin(49.88475°), \cos(49.88475°)\sin(−64°) - \tan(17°)\cos(−64°)) \\ \| = \arctan(0.7647500, 0.6443271×−0.898794 - 0.3057307×0.4383711) \\ \| = \arctan(0.7647500, −0.7131409) = 133° \end{align*}

As expected, \( H = −141.9946° \) yields the wrong \( A \), and \( H = 49.88475° \) yields the right \( A \).

Astronomy AnswersPhenomena in the Horizontal Coordinate System