# Sun rise and set, and twilight

[ Root ]
```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.
```

## Overview and method

[Top]

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;

```       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;

1. 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
2. Guess a figure for the UT at which the sunrise occurs
3. Calculate the Hour angle and declination at the Greenwich meridian for the Sun at the time of the current guess,
4. 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,
5. Repeat steps 3 and 4 until there is no significant difference between the succesive estimates of times for the sunrise
6. 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.

## Formulas and example calculation

[Top]

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.

### Centuries since J2000.0

[Top]

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

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
```

### Sun's position

[Top]

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
```

### Correction term and new estimate of time

[Top]

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.436916
```
and 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.
```

### Recalculation for second approximation

[Top]

The above is not a bad approximation to the true time that day, but we can calculate the next approximation by;

• re-calculating the number of centuries since J2000.0 using the new estimate of UT
• recalculating the position of the Sun for the new t value,
• calculating a new correction term 
• using the iteration formula  to calculate the new time
Below is a summary of this second iteration, the formulas for the Sun's position are identical to those above.
```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.

[Top]

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.

### Input area

[Top]

Put the following labels in A5 to A12;

```Year
Month
Day
Latitude
Longitude
time zone
Required alt
Rise or set
```
and 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.

### Output area part one

[Top]

Put the following labels into cells D5 to D8;

```Event UT
Event zone time
Days since J2000
Centuries
```
and 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/36525
```
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.

### Calculation area, first approximation

[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))/36525
```
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);
```3.720725
5.066172
-0.031580
3.689146
0.069088
0.409096
-0.208593
0.069088
0.251446
1.316622
1.789335
-0.011861
```

### Calculation area, further approximations

[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)/36525
```
And 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.011861
```
You 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.011861
```
Four 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.

### Output area part 2

[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 + 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.

## VBA function sunrise()

[Top]

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.
```'   define some numerical constants - these are not

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
'   from the Excel order, so the order matches most
'   textbooks

Function DegSin(x As Double) As Double
End Function

Function DegCos(x As Double) As Double
End Function

Function DegTan(x As Double) As Double
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
```

• 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

## QBASIC program

[Top] For those who prefer a program to a spreadsheet, below is the `QBASIC` listing for a simple program which will prompt the user for;
• latitude (decimal degrees, North positive) of observer
• longitude (decimal degrees, West negative)
• time zone offset from Greenwich, West negative
• rising or setting events (+1 for rise, -1 for set)
• desired altitude for Sun (i.e. -12 is nautical twilight)
• Year
• Month
• Day
and print the time of the selected event. The program prints the time in hhmm 24 hour format using the function FNdegmin\$, re-cycled from my Burnham precession program.

The easiest way to transfer the listing into QBASIC would be to;

• highlight (select) all the code betwen the lines of `*******`
• copy this text to the clipboard (cntrol key and C)
• paste the code into Notepad
• save, setting the file type to `All files (*.*)`, and using the file extension `.BAS`
• Check the file for broken lines, which occur if the browser window is set less wide than the width of some of the wider lines of BASIC
• Load resulting file into `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
'
'   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
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
```

## References

[Top]

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).

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 ]