Approximate astronomical positions

Sun rise and set, and twilight

[Root]

Contents

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

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

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

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;

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.

Spreadsheet

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

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;

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

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

'********************************************************************
'   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

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

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
keith@xylem.demon.co.uk