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
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 degreesThe 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)  sin(h) - sin(phi) * sin(delta) cosc = --------------------------------  cos(phi) * cos(delta) If cosc > 1 then correction = 0 Else If cosc < -1 then correction = 180 Else Correction = acos(cosc)Formula  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 ;
(setting) UT' = UT - (GHA + long - correction) Formula  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  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 ;
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.436916and this leads to the refined estimate for the time of Sunrise;
UT' = UT - (GHA + long + correction)  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.
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.
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
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.
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)/36525into a cell, then press enter.
-0.01186858316, i.e. the number of centuries since J2000.0 for 0h 1998 October 25th
QBASIClisting 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
QBASICor 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
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).
Paul Schlyter's page contains a restatement of a method similar to
the one used here, but applicable to his planetary and lunar positions
Last Modified 1999 May 30th
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.