1st Sept 01 People have reported errors when using the spreadsheet, so treat this implementation with caution. There is a method that works for both Sun and Moon at any latitude that you might find more reliable.
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
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;
Astronomical twilight altitude -18 degrees
Nautical twilight altitude -12 degrees
Civil twilight altitude -6 degrees
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;
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.
We find the number of days since J2000.0 (in this case a negative number) and then divide by 36525 to find the number of (julian) centuries. The table below can be used to find the days since J2000.0;
Calculating the centuries from J2000
------------------------------------
The tables below can be used to calculate the number of days and
the fraction of a day since the epoch J2000. If you need the
number of Julian centuaries, then just divide the 'day number'
by 36525.
Table A | Table B
Days to beginning of | Days since J2000 to
month | beginning of each year
|
Month Normal Leap | Year Days | Year Days
year year | |
| |
Jan 0 0 | 1998 -731.5 | 2010 3651.5
Feb 31 31 | 1999 -366.5 | 2011 4016.5
Mar 59 60 | 2000 -1.5 | 2012 4381.5
Apr 90 91 | 2001 364.5 | 2013 4747.5
May 120 121 | 2002 729.5 | 2014 5112.5
Jun 151 152 | 2003 1094.5 | 2015 5477.5
Jul 181 182 | 2004 1459.5 | 2016 5842.5
Aug 212 213 | 2005 1825.5 | 2017 6208.5
Sep 243 244 | 2006 2190.5 | 2018 6573.5
Oct 273 274 | 2007 2555.5 | 2019 6938.5
Nov 304 305 | 2008 2920.5 | 2020 7303.5
Dec 334 335 | 2009 3286.5 | 2021 7669.5
Worked Example
To find the number of days from J2000.0 for 1998 October 25th,
1. find from table A the number of days to the beginning of
October from the start of the year, here 273 days
4. write down the day number within the month, here 25 above
5. find from table B the days since J2000.0 to the beginning of
the year, here -731.5
6. add these three numbers.
For the date above;
273 + 25 + -731.5 = -433.5 days from J2000.0
7. to find the number of centuries t, divide by 36525,
t = -433.5 / 36525 = -0.01186858316
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;
L = 280.460 + 36000.770 * t (Mean longitude including aberration)
= -146.81813257 + 360
= 213.1818674
G = 357.528 + 35999.050 * t (Mean anomaly)
= -69.7297186 + 360
= 290.2702814
ec = 1.915 * sin(G) + 0.020 * sin(2*G) (eq centre correction)
= -1.7964017 - 0.0129997
= -1.8094014
lambda = L + ec (ecliptic longitude of Sun)
= 213.1818674 - 1.8094014
= 211.3724660
E = -ec + 2.466 * sin (2 * lambda) - 0.053 * sin(4 * lambda)
= 1.8094014 + 2.1922164 - 0.0431536
= 3.9584642
GHA = UTo - 180 + E,
where UTo is the current estimate of the time of Sunrise expressed in
degrees. For this first iteration, this gives
GHA = 180 - 180 + E (Hour angle of Sun on Greenwich meridian)
= 3.9584642
We also need the Sun's declination (delta), for which we need the
obliquity of the ecliptic (tilt of the Earth's axis);
Obl = 23.4393 - 0.0130 * t
= 23.43945429
delta = asin ( sin (obl) * sin(lambda)) (Sun's declination)
= -11.9515161
The guts of this iterative algorithm are the two formulas below;
UT' = UT - (GHA + long + correction) [1]
sin(h) - sin(phi) * sin(delta)
cosc = -------------------------------- [2]
cos(phi) * cos(delta)
If cosc > 1 then correction = 0
Else
If cosc < -1 then correction = 180
Else
Correction = acos(cosc)
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];
sin(h) - sin(phi) * sin(delta)
cosc = --------------------------------
cos(phi) * cos(delta)
= -0.1453808 - 0.79335234 * -0.20708391
---------------------------------------
0.608761429 * 0.978323186
= 0.14975242
------------ = 0.2514458
0.59556542
acos(cosc) = correction = 75.436916
and this leads to the refined estimate for the time of Sunrise; UT' = UT - (GHA + long + correction) [1]
UT' = 180 - (3.9584642 + -1.9167 + 75.436916)
= 102.52132 degrees
= 6.83475 hrs UT.
The above is not a bad approximation to the true time that day, but we can calculate the next approximation by;
New figure for centuries;
t = (-433.5 + 102.52132/360) /36525
= -0.0118607863
Sun's modified position
L = 213.4625604
G = 290.5509609
ec = - 1.7931300 - 0.0131480 = -1.806278
lambda = 211.6562823
E = 1.806278 + 2.2032968 - 0.0425355
= 3.9670393
GHA = ut' - 180 + E
= 102.52132 - 180 + 3.9670393
= - 73.5116407
New obliquity and declination
Obl = 23.43945419 (not very significant change!)
delta = -12.04991163
Correction term
cosc = -0.1453808 - 0.79335234 * -0.208763698
---------------------------------------
0.608761429 * 0.977966113
(notice that only the values of cos and sin of delta
actually change, the other numbers stay the same)
correction = acos(cosc) = 75.298919
New estimate of UT
ut'' = ut' - (GHA + long + correction)
ut'' = 102.52132 - (-73.5116407 - 1.9167 + 75.298919)
= 102.650741
= 6.84338 hrs UT
= 6 h 50 m 36 sec
When I did this calculation myself (using a basic scientific calculator and rounding answers in a convenient if unsystematic way) I used a column layout, with a new column for the second iteration.
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;
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.
Put the following labels in A5 to A12;
Year Month Day Latitude Longitude time zone Required alt Rise or setand put the 'test values' in B5 to B12;
1998 10 25 52.5 -1.9167 0 -0.833 (-6, -12 or -18 for twilights) 1 (-1 for setting events instead of rising events)These values correspond to 1998 October 25th, Birmingham, no zone offset, finding Sunrise (upper limb in contact with mathematical horizon), and a rise event.
[Top]
Put the following labels into cells D5 to D8;
Event UT Event zone time Days since J2000 Centuriesand put the following formulas into E7 and E8,
=367*B5-INT(7*(B5+INT((B6+9)/12))/4)+INT(275*B6/9)+B7-730531.5 =E7/36525These 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.
[Top]
Now put the following labels into cells A17 to A28, these names should tie in with the example calculations above;
L G ec lambda E obl delta GHA cosc correction utnew new centuries
Now we can put the first column of formulas into cells B17 to B28. You should be able to copy and paste the formulas into a spreadsheet as one block;
=MOD(4.8949504201433+628.331969753199*$E$8,6.28318530718) =MOD(6.2400408+628.3019501*$E$8,6.28318530718) =0.033423*SIN(B18)+0.00034907*SIN(2*B18) =B17+B19 =0.0430398*SIN(2*B20) - 0.00092502*SIN(4*B20) - B19 =0.409093-0.0002269*$E$8 =ASIN(SIN(B22)*SIN(B20)) =3.14159265358979 - 3.14159265358979 + B21 =(SIN(0.017453293*$B$11) - SIN(0.017453293*$B$8)*SIN(B23))/(COS(0.017453293*$B$8)*COS(B23)) =ACOS(B25) =3.14159265358979 - (B24+0.017453293*$B$9 + $B$12*B26) =($E$7 + B27/(6.28318530718))/36525Notice 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);
3.720725 5.066172 -0.031580 3.689146 0.069088 0.409096 -0.208593 0.069088 0.251446 1.316622 1.789335 -0.011861
[Top]
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;
=MOD(4.8949504201433+628.331969753199*B28,6.28318530718) =MOD(6.2400408+628.3019501*B28,6.28318530718) =0.033423*SIN(C18)+0.00034907*SIN(2*C18) =C17+C19 =0.0430398*SIN(2*C20) - 0.00092502*SIN(4*C20) - C19 =0.409093-0.0002269*B28 =ASIN(SIN(C22)*SIN(C20)) =B27 - 3.14159265358979 + C21 =(SIN(0.017453293*$B$11) - SIN(0.017453293*$B$8)*SIN(C23))/(COS(0.017453293*$B$8)*COS(C23)) =ACOS(C25) =B27 - (C24+0.017453293*$B$9 + $B$12*C26) =($E$7 + C27/6.28318530718)/36525And using the test values, I got the following values in cells C17 to C28;
3.725625 5.071071 -0.031525 3.694099 0.069238 0.409096 -0.210311 -1.283020 0.253776 1.314214 1.791594 -0.011861You can now copy the formulas from C17:C28 into D17:D28 and into E17:E28 to provide the next iterations in the calculation. I get the following values in D17:D28, and E17:E28 when doing this;
3.725631 3.725631 5.071077 5.071077 -0.031525 -0.031525 3.694105 3.694105 0.069238 0.069238 0.409096 0.409096 -0.210313 -0.210313 -1.280761 -1.280758 0.253779 0.253779 1.314211 1.314211 1.791597 1.791597 -0.011861 -0.011861Four iterations are usually enough to get a reliable result, and as you can see the results in C27 and D27 (1.791597) are identical to 6 dp in this case.
[Top]
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 + B10Cell 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;
' define some numerical constants - these are not
' accessible in the spreadsheet.
Public Const pi As Double = 3.14159265358979
Public Const tpi As Double = 6.28318530717958
Public Const degs As Double = 57.2957795130823
Public Const rads As Double = 1.74532925199433E-02
' The trig formulas working in degrees. This just
' makes the spreadsheet formulas a bit easier to
' read. DegAtan2() has had the arguments swapped
' from the Excel order, so the order matches most
' textbooks
Function DegSin(x As Double) As Double
DegSin = Sin(rads * x)
End Function
Function DegCos(x As Double) As Double
DegCos = Cos(rads * x)
End Function
Function DegTan(x As Double) As Double
DegTan = Tan(rads * x)
End Function
Function DegArcsin(x As Double) As Double
DegArcsin = degs * Application.Asin(x)
End Function
Function DegArccos(x As Double) As Double
DegArccos = degs * Application.Acos(x)
End Function
Function DegArctan(x As Double) As Double
DegArctan = degs * Atn(x)
End Function
Function DegAtan2(y As Double, x As Double) As Double
' this returns the angle in the range 0 to 360
' instead of -180 to 180 - and swaps the arguments
' This format matches Meeus and Duffett-Smith
DegAtan2 = degs * Application.Atan2(x, y)
If DegAtan2 < 0 Then DegAtan2 = DegAtan2 + 360
End Function
Private Function range2pi(x)
'
' returns an angle x in the range 0 to two pi rads
' This function is not available in the spreadsheet
'
range2pi = x - tpi * Int(x / tpi)
End Function
Private Function range360(x)
'
' returns an angle x in the range 0 to 360
' used to prevent the huge values of degrees
' that you get from mean longitude formulas
'
' this function is private to this module,
' you won't find it in the Function Wizard,
' and you can't use it on a spreadsheet.
' If you want it on the spreadsheet, just remove
' the 'private' keyword above.
'
range360 = x - 360 * Int(x / 360)
End Function
Function degdecimal(d, m, s)
' converts from dms format to ddd format degdecimal = d + m / 60 + s / 3600
End Function
Function day2000(year As Integer, month As Integer, day As Integer,
hour As Integer, _
min As Integer, sec As Double, Optional greg) As Double
' returns days before J2000.0 given date in gregorian calender (greg=1)
' or julian calendar (greg = 0). If you don't provide a value for greg,
' then assumed Gregorian calender
Dim a As Double
Dim b As Integer
If (IsMissing(greg)) Then greg = 1
a = 10000# * year + 100# * month + day
If (month <= 2) Then
month = month + 12
year = year - 1
End If
If (greg = 0) Then
b = -2 + Fix((year + 4716) / 4) - 1179
Else
b = Fix(year / 400) - Fix(year / 100) + Fix(year / 4)
End If
a = 365# * year - 730548.5
day2000 = a + b + Fix(30.6001 * (month + 1)) + day + (hour + min /60 + sec / 3600) / 24
End Function
Function Sunrise(ByVal day As Double, _
glat As Double, glong As Double, index As Integer, _
Optional altitude) As Double
' Takes the day number of 0h UT on a given day, and returns the
' time of sunrise or set according to index 1 or 2. The optional
' argument can be used to give an arbitrary value for the altitude
' of the Sun. The default value corresponds to the top limb of the
' Sun touching the mathematical horizon, allowing for refraction
' at sea level. Method from 9.311 and 9.33 of the Explanatory
' Suppliment to the Astronomical Almanac (Seidelmann, 1992).
Dim sinalt As Double, gha As Double, lambda As Double, delta As Double, _
t As Double, c As Double, days As Double, utold As Double, utnew As Double, _
sinphi As Double, cosphi As Double, _
L As Double, G As Double, E As Double, obl As Double, signt As Double, _
act As Double
If (IsMissing(altitude)) Then altitude = -0.833
utold = 180 'initial value finds first position at 12h UT on day
utnew = 0
sinalt = Sin(altitude * rads) 'go for the sunrise/sunset altitude first
sinphi = DegSin(glat)
cosphi = DegCos(glat)
If index = 1 Then signt = 1 Else signt = -1
Do While Abs(utold - utnew) > 0.1
utold = utnew 'carry forward the iteration
days = day + utold / 360 'update the arguments for sun's position
t = days / 36525
' find sun's coordinates
L = range360(280.46 + 36000.77 * t)
G = 357.528 + 35999.05 * t
lambda = L + 1.915 * DegSin(G) + 0.02 * DegSin(2 * G)
E = -1.915 * DegSin(G) - 0.02 * DegSin(2 * G) + 2.466 * DegSin(2 * lambda) - 0.053 * DegSin(4 * lambda)
obl = 23.4393 - 0.13 * t
' update the hour angle and declination of the Sun
gha = utold - 180 + E
delta = DegArcsin(DegSin(obl) * DegSin(lambda))
' calculate the correction term for this loop
c = (sinalt - sinphi * DegSin(delta)) / (cosphi * DegCos(delta))
act = DegArccos(c)
If c > 1 Then act = 0
If c < -1 Then act = 180
' find the newly corrected ut for sunrise/set
utnew = range360(utold - (gha + glong + signt * act))
Loop
Sunrise = utnew
End Function
To add these user defined functions to your spreadsheet, you should;
=day2000(1998,10,25,0,0,0)/36525 into a cell, then press
enter.-0.01186858316, i.e. the
number of centuries since J2000.0 for 0h 1998 October 25th QBASIC listing for a simple program which will prompt the user for;
The easiest way to transfer the listing into QBASIC would be to;
*******
All files (*.*), and using the
file extension .BAS
QBASIC or the FirstBas basic
compiler. '********************************************************************
' Sun rise/set /twilight
' From Explanatory Supplement
' definitions and functions
DEFDBL A-Z
pi = 4 * ATN(1)
tpi = 2 * pi
twopi = tpi
degs = 180 / pi
rads = pi / 180
'
' the function below returns the true integer part,
' even for negative numbers
'
DEF FNipart (x) = SGN(x) * INT(ABS(x))
'
' FNday only works between 1901 to 2099 - see Meeus chapter 7
'
DEF FNday (y, m, d, h) = 367 * y - 7 * (y + (m + 9) \ 12) \ 4 + 275 * m \ 9 + d - 730531.5 + h / 24
'
' define some arc cos and arc sin functions
'
DEF FNacos (x)
s = SQR(1 - x * x)
FNacos = ATN(s / x)
END DEF
DEF FNasin (x)
c = SQR(1 - x * x)
FNasin = ATN(x / c)
END DEF
'
' the function below returns an angle in the range
' 0 to two pi
'
DEF FNrange (x)
b = x / tpi
a = tpi * (b - FNipart(b))
IF a < 0 THEN a = tpi + a
FNrange = a
END DEF
'
' the function below returns the time in 24 hour
' notation - hhmm - given a decimal hours figure.
DEF FNdegmin$ (d)
c = ABS(d)
a = INT(c)
b = 60 * (c - a)
d$ = STR$(INT(a * 100 + b + .5))
d$ = RIGHT$(d$, LEN(d$) - 1)
WHILE LEN(d$) < 4
d$ = "0" + d$
WEND
FNdegmin$ = d$
END DEF
CLS
'
' get the date and geographical position from user
'
INPUT " lat : ", glat
INPUT " long : ", glong
INPUT " zone : ", zone
INPUT "rise/set : ", riset
INPUT " Sun alt : ", altitude
INPUT " year : ", y
INPUT " month : ", m
INPUT " day : ", d
day = FNday(y, m, d, 0)
PRINT " day no : "; day
utold = pi
utnew = 0
sinalt = SIN(altitude * rads) 'go for the sunrise/sunset altitude first
sinphi = SIN(glat * rads)
cosphi = COS(glat * rads)
glong = glong * rads
DO WHILE ABS(utold - utnew) > .001
utold = utnew
days = day + utold / tpi
t = days / 36525
' get arguments of Sun's orbit
L = FNrange(4.8949504201433# + 628.331969753199# * t)
G = FNrange(6.2400408# + 628.3019501# * t)
ec = .033423 * SIN(G) + .00034907# * SIN(2 * G)
lambda = L + ec
E = -1 * ec + .0430398# * SIN(2 * lambda) - .00092502# * SIN(4 * lambda)
obl = .409093 - .0002269# * t
delta = FNasin(SIN(obl) * SIN(lambda))
GHA = utold - pi + E
cosc = (sinalt - sinphi * SIN(delta)) / (cosphi * COS(delta))
SELECT CASE cosc
CASE cosc > 1
correction = 0
CASE cosc < -1
correction = pi
CASE ELSE
correction = FNacos(cosc)
END SELECT
utnew = FNrange(utold - (GHA + glong + riset * correction))
LOOP
PRINT " UT : "; FNdegmin$(utnew * degs / 15)
PRINT " zone : "; FNdegmin$(utnew * degs / 15 + zone)
END
'******************************************************************
Below is the input and output of the program given the example case used throughout this page, sunrise on 1998 October 25th at Birmingham UK.
lat : 52.5 long : -1.9167 zone : 0 ise/set : 1 Sun alt : -0.833 year : 1998 month : 10 day : 25 day no : -433.5 UT : 0651 zone : 0651
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. Relevant sections are 9.311 (p484) and 9.33
(p486).
http://193.12.249.96/pausch/comp/ppcomp.html
Paul Schlyter's page contains a restatement of a method similar to the one
used here, but applicable to his planetary and lunar positions method.
[Root
]
Last Modified 1999 May 30th
Keith Burnett
kburnett@btinternet.com