How to compute rise/set times and altitude above horizon

By Paul Schlyter, Stockholm, Sweden
email: pausch@stjarnhimlen.se or WWW: http://stjarnhimlen.se/

Break out of a frame

How to compute planetary positions
Tutorial with numerical test cases

1. Computing the Sun's altitude above the horizon

First we compute the Sun's RA and Decl for this moment, as outlined earlier. Now we need to know our Local Sidereal Time. We start by computing the sidereal time at Greenwich at 00:00 Universal Time, let's call this quantity GMST0:

    GMST0 = L + 180
L is the Sun's mean longitude, which we compute as: 
    L = M + w
Note that we express GMST in degrees here to simplify the computations. 360 degrees of course corresponds to 24 hours, i.e. each hour corresponds to 15 degrees.

Now we can compute our Local Sidereal Time (LST):

    LST = GMST0 + UT*15.0 + long
UT is the Universal Time, expressed in hours+decimals, the remaining quantities are expressed in degrees. To convert UT to degrees we must multiply it by 15 above. long is our local longitude in degrees, where east longitude counts as positive, and west longitude as negative. (this is according to the geographic standard, and the recent astronomical standard; if you prefer to use the older astronomical standard where west longitude counts as positive, then you must change the '+' in front of 'long' to a '-' above).

Next let's compute the Sun's Local Hour Angle (LHA), i.e. the angle the Earth has turned since the Sun last was in the south:

    LHA = LST - RA
A negative hour angle means the Sun hasn't been in the south yet, this day. The angle -10 degrees is of course the same as 350 degrees, i.e. adding or subtracting even multiples of 360 degrees does not change the angle.

We also need to know our latitude (lat), where north latitude counts as poisitive and south latitude as negative. Now we can compute the Sun's altitude above the horizon:

    sin(h) = sin(lat) * sin(Decl) + cos(lat) * cos(Decl) * cos(LHA)
We compute sin(h), and then take the arcsine of this to get h, the Sun's altitude above the horizon.

2. Computing the Sun's rise/set times.

This is really the inverse of the previous problem, where we computed the Sun's altitude at a specific moment. Now we want to know at which moment the Sun reaches a specific altitude.

First we must decide which altitude we're interested in:

h = 0 degrees: Center of Sun's disk touches a mathematical horizon
h = -0.25 degrees: Sun's upper limb touches a mathematical horizon
h = -0.583 degrees: Center of Sun's disk touches the horizon; atmospheric refraction accounted for
h = -0.833 degrees: Sun's upper limb touches the horizon; atmospheric refraction accounted for
h = -6 degrees: Civil twilight (one can no longer read outside without artificial illumination)
h = -12 degrees: Nautical twilight (navigation using a sea horizon no longer possible)
h = -15 degrees: Amateur astronomical twilight (the sky is dark enough for most astronomical observations)
h = -18 degrees: Astronomical twilight (the sky is completely dark)

As you can see, there are several altitides to choose among. In most countries an altitude of -0.833 degrees is used to compute sunrise/set times (Sun's upper limb touches the horizon; atmospheric refraction accounted for). One exception is the Swedish national alamancs, which use -0.583 degrees (Center of Sun's disk touches the horizon; atmospheric refraction accounted for) - however my own Swedish almanac Stjärnhimlen ("The Starry Sky") uses the international convention of -0.833 degrees.

When we've decided on some value for the altitude above the horizon, we start by computing the Sun's RA at noon local time. When the Local Sidereal Time equals the Sun's RA, then the Sun is in the south:

    LST  =  RA
which yields:

    GMST0 + UT*15.0 + long  = RA
Since we know GMST, long, and RA, it's now a simple matter to compute UT (GMST0 should also be computed at noon local time):

    UT_Sun_in_south = ( RA - GMST0 - long ) / 15.0
Now we're going to compute the Sun's Local Hour Angle (LHA) at rise/set (or at twilight, if we've decided to compute the time of twilight). This is the angle the Earth must turn from sunrise to noon, or from noon to sunset:

                sin(h) - sin(lat)*sin(Decl)
    cos(LHA) = -----------------------------
                   cos(lat) * cos(Decl)
If cos(LHA) is less than -1.0, then the Sun is always above our altitude limit. If we were computing rise/set times, the Sun is then aboute the horizon continuously; we have Midnight Sun. Or, if we computed a twilight, then the sky never gets dark (a good example is Stockholm, Sweden, at midsummer midnight: the Sun then only reaches about 7 degrees below the horizon: there will be civil twilight, but never nautical or astronomical twilight, at midsummer in Stockholm).

If cos(LHA) is greater than +1.0, then the Sun is always below our altitude limit. One example is midwinter in the arctics, when the Sun never gets above the horizon.

If cos(LHA) is between +1.0 and -1.0, then we take the arccos to find LHA. Convert from degrees to hours by dividing by 15.0

Now, if we add LHA to UT_Sun_in_south, we get the time of sunset. If we subtract LHA from UT_Sun_in_south, we get the time of sunrise.

Finally, we convert UT to our local time.

3. Higher accuracy: iteration

The method outlined above only gives an approximate value of the Sun's rise/set times. The error rarely exceeds one or two minutes, but at high latitudes, when the Midnight Sun soon will start or just has ended, the errors may be much larger. If you want higher accuracy, you must then iterate, and you must do a separate iteration for sunrise and sunset:

a) Compute sunrise or sunset as above, with one exception: to convert LHA from degrees to hours, divide by 15.04107 instead of 15.0 (this accounts for the difference between the solar day and the sidereal day. You should only use 15.04107 if you intend to iterate; if you don't want to iterate, use 15.0 as before since that will give an approximate correction for the Earth's orbital motion during the day).

b) Re-do the computation but compute the Sun's RA and Decl, and also GMST0, for the moment of sunrise or sunset last computed.

c) Iterate b) until the computed sunrise or sunset no longer changes significantly. Usually 2 iterations are enough, in rare cases 3 or 4 iterations may be needed.

d) Make sure you iterate towards the sunrise or sunset you want to compute, and not a sunrise or sunset one day earlier or later. If the computed rise or set time is, say, -0.5 hours local time, this means that it really happens at 23:30 the day before. If you get a value exceeding 24 hours local time, it means it happens the day after. If this is what you want, fine, otherwise add or subtract 24 hours. This only becomes a problem when there soon will be, or just has been, Midnight Sun.

e) Above the arctic circle, there are occasionally two sunrises, or two sunsets, during the same calendar day. Also there are days when the Sun only sets, or only rises, or neither rises nor sets. Pay attention to this if you don't want to miss any sunrise or sunset in your computations.

f) If you compute twilight instead of rise/set times, e) applies to the "twilight arctic circle". The normal arctic circle is situated at 66.7 deg latitude N and S (65.9 deg if one accounts for atmospheric refraction and the size of the solar disk). The "twilight arctic circle" is situated 6, 12, 15 or 18 deg closer to the equator, i.e. at latitude 60.7, 54.7, 51.7 or 48.7 degrees, depending on which twilight you're computing.

4. Computing the Moon's rise/set times.

This is really done the same way as the Sun's rise/set times, the only difference being that you should compute the RA and Decl for the Moon and not for the Sun. However, the Moon moves quickly and its rise/set times may change one or even two hours from one day to the next. If you don't iterate the Moon's rise/set times, you may get results which are in error by up to an hour, or more.

Another thing to consider is the lunar parallax, which affects the Moon's rise/set time by several minutes or more. One way to deal with the lunar parallax is to always use the Moon's topocentric RA and Decl. Another, simpler, way is to use the Moon's geocentric RA and Decl and instead adjust h, the rise/set altitude, by decreasing it by m_par, the lunar parallax. If you want to compute rise/set times for the Moon's upper limb rather than the center of the Moon's disk, you also need to compute m_sd, the semi-diameter or apparent radius of the Moon's disk in the sky. Note that the Moon's upper limb may for some lunar phases and circumstances be on the dark part of the Moon's disk

Thus you choose your h for Moon rise/set computation like this:

h = -m_par: Center of Moon's disk touches a mathematical horizon
h = -(m_par+m_sd): Moon's upper limb touches a mathematical horizon
h = -0.583 degrees - m_par: Center of Moon's disk touches the horizon; atmospheric refraction accounted for
h = -0.583 degrees - (m_par+m_sd): Moon's upper limb touches the horizon; atmospheric refraction accounted for

Yet another thing to consider: the Sun is always in the south near 12:00 local time, but the Moon may be in the south at any time of the day (or night). This means you must pay more attention that you're really iterating towards the rise or set time you want, and not some rise/set time one day earlier, or later.

Since the Moon rises and sets on the average 50 minutes later each day, there usually will be one day each month when the Moon never rises, and another day when it never sets. You must have this in mind when iterating your rise/set times, otherwise your program may easily get caught into an infinite loop when it tries to force e.g. a rise time between 00:00 and 24:00 local time on a day when the Moon never rises.

At high latitudes the Moon occasionally rises, or sets, twice on a single calendar day. This may happen above the "lunar arctic circle", which moves between 61.5 and 71.9 deg latitude during the 18-year period of the motion of the lunar nodes. You may want to pay attention to this.

Yes, computing the Moon's rise/set times is unfortunately messy, much due to its quick orbital motion around the Earth.

5. Computing rise/set times for other celestial bodies.

This is done the same way as for the Sun, with some differences:

a) Compute the RA and Decl for that body instead of for the Sun. If the body is a star, get its RA and Decl from a suitable star catalog.

b) GMST0 is still needed, so you should compute the Sun's mean longitude.

c) Always use 15.04107 instead of 15.0 when converting LHA from degrees to hours.

Since the planets move much slower than the Moon, and the stars hardly move at all, one need not iterate. If one wants high accuracy, one may find it worthwhile to iterate the rise/set times for Mercury, Venus and Mars (these are the planets that move most quickly).