This page describes an algorithm that will enable you to calculate the time of Moon rise and set and Sun rise and set, as well as the twilight times, for any latitude on the Earth's surface. The algorithm presented here will warn you if the Moon or Sun remains above or below the horizon all day (a definite possibility above latitudes of 66.5 degrees North), or if the object does not rise or set within the day chosen (the Moon will exhibit this behaviour on one day in each lunar cycle). The algorithm is taken from Montenbruck and Pfleger's Astronomy on the Personal Computer second English edition, section 3.8.

If you just want the *results* for your observing location, then
I recommend the Web page
based programs provided by the US Naval Observatory. One of these
programs will give you a compact yearly calendar style printout of rise
and set times at a given location on a single sheet of paper. I use a
set of printouts in my own observation planning - its easier than
switching on the computer!

In the * implementation* of the algorithm on this page,
I have used approximate routines to calculate the coordinates of
the Sun and Moon. You can choose latitudes where the program will
predict a very short day when the Sun really remains below the
horizon all day - this inaccuracy affects latitudes close to
67.43 degrees North, and you can imagine a 'fuzzy zone' around
the Earth at this latitude (and a mirror image one in the South)
where the program may be in error. The situation for moonrise and
set is more complex as the 'lunar arctic circle' changes as the
Moon's orbit changes its orientation. I suspect that more
accurate position routines and better approximations for
topocentric parallax, refraction and the diameter of the Sun and
Moon would shrink this 'fuzzy zone', at the cost of increased
complexity. The program will give meaningful results *above* latitude
70 degrees where methods based on iteration cannot.

The algorithm is implemented in `QBASIC`

, a well
structured interpreted basic available free for Windows computers, and
easily adapted to Visual Basic syntax. The implementation is based very
closely on `PASCAL`

routines taken from Montenbruck and
Pfleger. I have also added a Javascript implementation.

The rising and setting of the Sun and Moon is defined here to be the
instant when the *upper limb* of the object is seen to cross a
mathematically flat horizon by an observer at sea level. Real observing
situations on land rarely approach this idealised model, and for my uses
calculating rise and set times to the nearest minute or so is fine.

The following factors need to be addressed

- The Sun and Moon vary in position against the stars over a period of one day - the Moon moves about 12 degrees per day on average.
- There are latitudes on Earth where the Sun and Moon can be above or below the horizon for days at a time. The algorithm has to detect and report these circumstances.
- The time of Moonset will be about 50 minutes later each night after New Moon. There will be at least one day in each lunation during which the Moon does not set, and one day where the Moon does not rise. The algorithm must report these days.
- You must calculate the altitude of the Sun and the Moon as seen by
the observer on the
*surface*of the Earth - allowing for parallax from the geocentric coordinates. - You must correct for the effects of atmospheric refraction (not relevant for the twilights).
- You must allow for the size of the Sun or Moon's disc (except for the twilights).

In Montenbruck and Pfleger's approximate implementation, the last three points above boil down to subtracting a single correction term from the geocentric altitude of the Moon or Sun when looking for rise or setting events.

- Sunrise occurs when the centre of the Sun's disc is at a geocentric altitude of -50 arc minutes.
- Moon rise occurs when the centre of the Moon's disc is at a geocentric altitude of +8 arc minutes.
- Nautical twilight is defined to be when the centre of the Sun is at a geocentric altitude of -12 degrees.

The large difference between the Moon and the Sun correction reflects the much larger parallax effect on Moon positions - roughly 1 degree for the Moon compared to 9 arc seconds for the Sun.

Possible algorithms include:-

- Use a single position and calculate as for a star
- Use a 'trial and improvement' or 'iterative' method
- Use a table based interpolation method

I used the 'single position and calculate as for star' method on a programmable calculator - the results for the Sun at moderate latitudes were 20 minutes or so wrong. Given the state of the local horizon, this was not much of a problem. This 'method' was very poor for the Moon, almost worse than guessing from the number of days since New Moon.

The 'trial and improvement' method with allowance for refraction, parallax and disc size is much better, and can give good results for the Sun for latitudes below 65N. The method will mostly converge after a few iterations - see the spreadsheet implementation on this Web site, or Paul Schlyter's page.

Tral and improvement methods get fiddly in situations where the object remains above or below the horizon for a number of days, or where the object does not set one day or the next.

A table based method could be very simple - just calculate the altitude of the Moon for each minute of the day and find when the altitude changes sign. There is your rising or setting event.

This simple table based method could detect the 'always above the horizon' or 'always below the horizon' situations as the (corrected) altitude would never change sign over the day. The method could also cope with the Moon missing a rise or set each lunar cycle, but 'grazing' situations, where the Sun never quite dips below the horizon, would be more difficult as there would be no sign change in the altitudes. A real program based on this system would be really slow, requiring at worst case 1440 calculations of altitude for each object for each day.

Interpolation based methods use a shorter table of altitudes - in my implementation we use a table with a maximum of 25 rows - and fit a quadratic function to each set of three altitudes in turn. We then use some GCSE algebra to find the times when the quadratic crosses the X axis - has zero altitude - and these times are candidates for a rising or setting event. This system can also deal with the situation when the rising and setting occur within the hour - the quadratic then has two zeros. Finally, if the two zeros are the same, then there must be a 'grazing' event in that interval, and the algorithm can detect the event.

The steps in the procedure are

- Calculate the sine of the altitude (corrected for refraction, parallax and the limb radius) for the first three hours in the day. Using the sine of the altitude saves some calculation and means that the values will always be within the range +1 to -1.
- Save sign of the altitude at 0h
- Take the values of the sine of the altitude for each set of three consecutive hours, correct for the refraction, parallax and fit a quadratic curve through these three points. See Montenbruck and Pfleger for algebra.
- Calculate the 'discriminant' for the quadratic, and classify as
- No roots within the 2 hour interval - move onto the next 2 hour range and start from step 1
- Two roots, one root within the interval - classify as rising or
setting event and appropriate flag
- If two flags set, i.e. rising and setting events both found, go to step 5
- If one flag set, then return to step 1 for next two hour interval, to search for remaining event.

- Two roots, both within the interval - set the rising and setting flags and go to step 5.
- One root (i.e. both roots same) within interval - set both rising and setting flags and go to step 5.

- If neither rise nor set flag is set, then warn user that object is 'always above' or 'always below' horizon depending on the sign of the altitude in step 2, and move onto next object or quit.
- Print rise and set times, or message if one event missing in current day.
- Move onto next object or quit.

You can download the tested `QBASIC`

code below. The
`QBASIC`

editor insists on adding functions and subroutines
to the end of the listing in alphabetical order - which does not make
reading the code any easier. I have added comments in the 'main' program
that should help you link the code to the description of the algorithm
above.

- Link to QBASIC code (12 K) (some long lines)
- Code in ZIP file format (4 K)

Below are some samples of program input and output, together with accurate results from Chris Marriott's SkyMap Pro 6.0 and with a year 2000 calendar for the Sun, Moon and Nautical Twilight generated from the USNO web based program.

Rise and set times for today in Birmingham UK - simple case with no
ambiguities. Below is the dialogue with the QBASIC program. Times of
events are given in form `hhmm`

, and are UT times.

Rise and set for Sun and Moon ============================= Year (yyyy) - - - - - - - - :2000 Month (mm) - - - - - - - - :1 Day (dd) - - - - - - - - :3 Time zone (East +) - - - - :0 Longitude (w neg, decimals) :-1.91667 Latitude (n pos, decimals) :52.5 Moon 501 1409 Sun 818 1606 Nautical twilight 653 1731

Skymap Pro gives the following times for the Sun and Moon on 2000 Jan 3 - agreement is fine to nearest minute, except for Moonrise that rounds the wrong way to give a result one minute different from the USNO tables.

Sunrise and twilight Date Rise Set Naut Twi E Naut Twi B 03 Jan 2000 08:18:11 16:05:50 17:30:56 06:53:06 Moonrise Rise: 5h 0m 23s Set: 14h 9m 18s

`QBASIC`

program gives a 'day' lasting 8
minutes:
Rise and set for Sun and Moon ============================= Year (yyyy) - - - - - - - - :1999 Month (mm) - - - - - - - - :12 Day (dd) - - - - - - - - :25 Time zone (East +) - - - - :0 Longitude (w neg, decimals) :0 Latitude (n pos, decimals) :67.43 Moon 1747 1206 Sun 1156 1204 Nautical twilight 747 1613

Skymap Pro shows that the Sun remains below the horizon for the whole
of that day at the latitude given. SkyMap Pro gives the altitude of the
(centre of) Sun at transit for this day as roughly 53 minutes below the
horizon. This implies an error of roughly 30 minutes in Sun position at
that time for the `QBASIC`

program. This seems rather large
and I will be puggling about with this a bit more.

The Moon times and times for twilight are within a minute of those
produced by SkyMap Pro. You might treat output from the
`QBASIC`

program presented here with caution for locations
with latitudes between 67 and 68 North.

For Narvik in Norway (further North), the program gives correct results give or take a minute:-

Rise and set for Sun and Moon ============================= Year (yyyy) - - - - - - - - :2000 Month (mm) - - - - - - - - :1 Day (dd) - - - - - - - - :3 Time zone (East +) - - - - :1 Longitude (w neg, decimals) :17.42 Latitude (n pos, decimals) :68.43 Moon 629 1157 Sun always below horizon Nautical twilight 743 1607SkyMap Pro gives

Moon Rise: 6h 27m 59s Set: 11h 57m 47s Date Rise Set Naut Twi E Naut Twi B 03 Jan 2000 --:--:-- --:--:-- 16:07:06 07:42:21

**Dec 2002 : **I have added a Javascript page with an
implementation of the algorithm. I have 'wrapped up' the main
calculation loop in a function and this makes the program easier
to handle. I have separated the moonrise and the sun related
events (twilight and sunrise) in different functions so as to
allow different tabulation intervals.

Sunrise times change relatively slowly, and a weekly listing is fine for a year. Moonrise times change by 50 minutes a day on average and so a daily listing is needed.

The Javascript page works fine on Netscape 4.7x, MS Internet Explorer 6 and Opera 7.

**Jan 2003 : **I have added a VBA module with an implementation of the
algorithm for sunrise / set and the various twilights. The module
makes various 'user defined functions' available in an Excel
spreadsheet (tested on Excel 97). See the VBA comments for
details.

This VBA application is still a 'work in progress' so let me know if you find any logic errors.

Oliver Montenbruck and Thomas Pfleger,
*Astronomy on the Personal Computer*

Springer

1994, 3rd edition

ISBN 3-540-63521-1

Section 3.8 of this book is the basis of the `QBASIC`

listing here. I have translated the `PASCAL`

routines into
`QBASIC`

one by one. This book does not contain an exhaustive
explanation of the quadratic fit method, and the test data the authors
provide will not test all branches of the code. This book comes with a
software disc with all the `PASCAL`

routines - but does not
cover all the ground in, say, Meeus' Astronomical
Algorithms.

Sky and Telescope's software pages contain listings in 'spaghetti
basic' that implement a similar method for the Sun and the Moon - there
is a clever fix used to avoid re-calculating positions for the Sun and
Moon up to 12 times each day at the cost of some increase in the 'fuzzy
zone' when the programs may predict a Sunrise and set when there is
none. The links below are to the `BASICA`

listings.

http://www.skypub.com/resources/software/basic/programs/sunup.bas

http://www.skypub.com/resources/software/basic/programs/moonup.bas

Seidelmann, P. Kenneth (ed)
*Explanatory Supplement to the Astronomical Almanac*

University Science Books

1992, completely revised

ISBN 0-935702-68-7

Definitive reference on all aspects of the ephemeris and associated calculations. No hostages taken, no example calculations and modern vectoral notation used throughout. Sections 9.311 (p484) and 9.33 (p486) present the iterative algorithm for Sunrise for moderate latitudes.

Duffett-Smith, Peter

Practical Astronomy with your calculator

Cambridge University Press

3rd edition 1988

ISBN 0-521-35699-7

Good treatment of the basics of coordinate systems. Presents the formulas for rise and set time of a star, and presents an interesting linear interpolation scheme for Moonrise.

[ Root ]

Last Modified 2002 Jan 2nd

Keith Burnett