This page contains a recipe for finding the times of events such as
sunrise, sunset, and the times for various definitions of twilight.
The recipe is illustrated by direct calculation, and implemented as

a 'portable' spreadsheet using only 'standard' functions

as an Excel enhanced spreadsheet using VBA functions

as a QBASIC program

Sunrise and sunset are often defined as the instant when the upper
limb of the Sun's disk is just touching the observer's mathematical
horizon assuming a spherical earth, and allowing for the atmospheric
refraction. This corresponds to an altitude of -0.833 degrees for the
Sun. The various 'twilights' are usually defined in terms of the Sun's
altitude as follows;

The limit of astronomical twilight is defined as when the light from
the Sun scattered by the atmosphere is roughly the same as that the
combined light from the stars, the zodiacal light and the
gegenschein. In my inner city area, the sky brightness never
drops to such a low level, so I use the nautical twilight to indicate
observing times.

The method implemented on this page is taken from the
Explanatory Supplement to the Astronomical Almanac
section 9.33, and uses a simple iterative scheme which will converge
on the times of events. This method may not converge above latitudes
of 60 degrees North, or below latitudes of 60 degrees South. The times
produced are found to be within seconds of times from a planatarium
program, and Dr Ahmed Monsur's Mooncalc. A similar method is explained
by Paul
Schlyter on his excellent page.

The method described on this page falls into a number of steps.
Suppose we want to calculate the time of Sunrise one autumn morning.
The process goes as follows;

find the number of days from 0h UT on the day in question to
J2000.0, and divide this by 36525 to find the number of julian
centuries, t

Guess a figure for the UT at which the sunrise occurs

Calculate the Hour angle and declination at the Greenwich meridian
for the Sun at the time of the current guess,

Calculate a correction term and use the term to arrive at a new
guess for the time of the sunrise, ensuring that the times stay within
the range 0 to 24 hours,

Repeat steps 3 and 4 until there is no significant difference
between the succesive estimates of times for the sunrise

This time is the time of the sunrise.

You can calculate the rise and set times of the Sun for any
day in the past or future 500 years to reasonable accuracy using an
ordinary pocket calculator. However, there is a large calculational
effort at each repetition of steps 3 and 4, and most people will use a
programmable calculator, a spreadsheet or a program function to
calculate the times of events.

As a concrete example, I shall calculate the time of the Sunrise on
1998 October 25th at Birmingham UK, latitude 52.5N, longitude 1.9167W.

We start by finding the number of centuries since J2000.0, as all
the formulas for the position of the Sun depend on this. As outlined
above, we calculate the position of the Sun for 0h on the day in
question. We use this position to calculate the time of Sunrise (even
though the Sun will have moved a fraction of a degree in the sky), and
then recalculate the position of the Sun for the new time. A refined
time for the Sunrise can then be calculated.

In all the calculations below, we measure time in degrees (180 degs
= 12h), and we take a crude first guess at the time of sunrise as
12h UT on the day (180).

The 'low precision' formulas below give the Sun's position to an
accuracy of about 0.1 degree over a few centuries either side of
J2000.0. Taking the figure of t = -0.01186858316 for the number of
Julian centuries since J2000.0 as calculated above for our example
date;

Formula [1] gives us a way of calculating a better estimate of the
time of sunrise, given the current estimate (ut) and the hour angle of
the Sun for that time, and the correction term. The quantity GHA +
long gives us the Sun's hour angle on the local meridian. Our first
guess for the time of sunrise is ut = 180.

For setting events (sunset, end of twilight) we subtract the
correction term in formula [1];

(setting) UT' = UT - (GHA + long - correction) [1]

Formula [2] gives us the value of the correction term for each
successive estimate. As we will see in the spreadsheet example below,
the correction terms converge rapidly to a single value under most
circumstances. You have to calculate the correction term first, before
you can use the iteration formula [1] to find your next estimate. I've
never really understood why the textbooks list them in this order!

For our example date of 1998 October 25th, at Birmingham (phi = 52.5,
long = -1.9167), we have for the correction term [2];

As a 'column based' layout seems so natural for this kind of
calculation, I have devised a simple spreadsheet based on the formulas
above. I have tried to use 'standard' formulas and so if you follow
the instructions below, you should be able to build a working
spreadsheet using just about any spreadsheet program you may have.

The sunrise.zip file contains
copies of this spreadsheet in .WKS and .SLK formats,
as well as the full MS Excel version. One of these should load into
most spreadsheet programs.

The main limitations of working in a 'portable' spreadsheet format
are;

all the trig functions work in radians!
I have converted the coefficients in the formulas to
give angles in radian measure

I can't make assumptions about being able to 'name' certain cells,
so I have to use 'absolute cell referencing'

the =mod(x,y) function changes implementation from
program to program, but this does not change the final answers.

To build the spreadsheet, I have an 'input area' with labels in A5
to A12 and corresponding values in B5 to B12, and an 'output area'
with labels in D5:D8 and the results of the calculations in E5:E8. The
main calculation consists of four successive approximations to the UT
of the event, and I put these calculations in the cells B17 to E28,
with labels for each term in A17 to A28. It might help if you look at
a screenshot of the basic spreadsheet
layout.

These formulas will give the number of days from J2000.0 in cell
E7and the number of Julian centuries from J2000.0 in
E8. I get -433.5 and -0.0118685832 in cells E7 and E8 using
the 1998 October 25th test date.

Notice the use of 'absolute cell referencing' for $E$7, $E$8, and a
few other cells. Note how I have kept a large number of decimal places
in the coeficients for the mean longitude (L) calculation. Any
rounding here causes trouble. Using our example data for Sunrise on
1998 October 25th, at Birmingham, 52.5 N, 1.9167 W, I get the
following
values in cells B17 to B28 (output formatted to 6 decimal places);

The next set of formulas in cells C17 to C28 calculates the next
approximation
to the time of Sunrise, these formulas are slightly different to the
first column
in that they pick up the new centuries figure from B28 and the new
guess for
UT from cell B27;

I put the time of Sunrise in cell E5 with the following formula to
convert the time
from radians to hours;

=E27*57.29577951/15
=E5 + B10

Cell E8 now holds the time of the event in decimal hours in the time
zone entered in cell B10 originally. For Birmingham on 1998 October
25th, I get 6.843 UT as the decimal hours of Sunrise, which
corresponds to 0651 hrs.

The 'portable' spreadsheet above will work on almost any program,
but it is inconvenient if you want to make a table of sunset or
sunrise times. Most modern business spreadsheets have a macro language
built into them, and Microsoft Excel comes with Visual Basic for
Applications. The VBA code below defines a range of 'user functions'
including;

degree based trig functions

days2000(year, month, day, hour, minute, second, optional
greg) which returns the number of days since J2000.0 for a given
instant. This function is valid for negative years and far into the
future, unlike the simple one line function used in the 'portable'
spreadsheet above. The optional argument 'greg' uses the Julian
calendar if set to 0, and the gregorian if set to 1 or not used. This
allows for countries like England and Sweeden, who did not adopt the
Gregorian calendar in 1582.

sunrise(days, glat, glong, index, optional altitude)
which returns the time of a rising or setting event. Index = +1 gives
a rising event, index = 2 gives a setting event. The optional argument
altitude allows you to specify the altitude you want the Sun to have,
i.e. a value of -12 for altitude and +1 for index will give you the
time when the night brightens to astronomical twilight.

An Excel 95 spreadsheet which contains a simple example of the use of
Sunrise() is available for
download. There are instructions about how to add the VBA code
below to an Excel spreadsheet after the code listing.

In the browser window, highlight all the text between the rows of
* below

Select Edit | Copy from the browser Edit menu

open Excel 95 with a blank workbook

click on the Insert menu and select Insert | Macro |
Module.
A new module window should appear

click in the new VBA module window, and select Edit | Paste
from the
Excel edit menu

Save the new workbook

Switch to a blank worksheet, and type the formula
=day2000(1998,10,25,0,0,0)/36525
into a cell, then press enter.
The cell value should be -0.01186858316, i.e. the number
of centuries since J2000.0 for 0h 1998 October 25th

Definitive reference on all aspects of the ephemeris and
associated calculations. No hostages taken, no example calculations
and modern vectoral notation used throughout. Relevant sections are
9.311 (p484) and 9.33 (p486).