This page is based on page C24 of the 1996 Astronomical Almanac which provides a method for finding the position of the Sun in the sky to an accuracy of 0.01 degree between the years 1950 and 2050. Positions are given in the equatorial coordinate system used by astronomers, and in the horizon coordinate system for a location with given latitude and longitude. You can use these formulas to work out the bearing and elevation of the Sun in your sky at a given time of day - horizon coordinates are accurate to about 1.5 arcmin (or 1/20th of the Sun's diameter). To put this level of accuracy in perspective, the Sun moves through about a quarter of a degree in the sky every minute. If you are not familiar with astronomical coordinate systems and the words used to describe them, then I would strongly recommend a visit to Nick Strobel's astronomy without a telescope page.
The formulas are based on an elliptical orbit for the Earth, using mean orbital elements and a two term approximation for the 'equation of centre'. There is also an approximate allowance made for the change in obliquity of the ecliptic with time, needed when converting to right ascension and declination. The positions are thus apparent positions, they are referred to the mean ecliptic and equinox of date.
I compared the positions found using this low precision formula with values referred to the mean ecliptic and equinox of date from a more accurate program. The results (for the whole 1950 to 2050 range) are summarised below. I found the series to be accurate within 3 seconds of RA and 15 arc seconds in declination.
Note added March 2001: I have had a number of e-mails from people who want to calculate the the azimuth (3 figure bearing from North) and altitude (elevation) of the Sun at a specific time and place. This is easily done by finding the local sidereal time for the place, finding the hour angle of the Sun and then converting to horizon coordinates. For convenience, I have added these formulas and examples of their use to this page. Altitude seems good to about 1.5 minutes of arc, and Azimuth seems good to about half an arcmin, but I have not tested this as fully as the RA and DEC. At the risk of being obvious, the altitude angle here is referred to a 'mathematical' horizon on a spherical Earth. This is a good approximation to the observerd horizon at sea, but poor in most land situations where buildings or hills will raise or lower the true horizon. If your application is critical, you need to take time to understand the coordinate systems in use, and check on the effects of refraction near the horizon.
Below, I give the formulas from page C24 of the Astronomical Almanac, with modified notation. I have given the formulas together with numerical values for a specific day. The calculations were done on a normal scientific calculator with 8 figure accuracy, except for the sidereal time calculation where I used an HP48 with 12 figures of precision.
Position of the Sun at 11:00 UT on 1997 August 7th in
Birmingham UK (52.5 latitude, -1.91667 longitude. Longitude is
always taken as negative to the West)
1. Find the days before J2000.0 (d) from the table
or by using a formula based on year, month and day formulas
(see the QBASIC program for one simple example)
d = 11/24 + 212 + 7 - 1096.5 = -877.04167
2. Find the Mean Longitude (L) of the Sun on this day
L = 280.461 + 0.9856474 * d
= -583.99284 + 720
(add multiples of 360 to bring in range 0 to 360)
= 136.00716
3. Find the Mean anomaly (g) of the Sun for this day
g = 357.528 + 0.9856003 * d
= -506.88453 + 720
= 213.11547
4. Find the ecliptic longitude (lambda) of the sun
lambda = L + 1.915 * sin(g) + 0.020 * sin(2*g)
= 134.97925
(note that the sin(g) and sin(2*g) terms constitute an
approximation to the 'equation of centre' for the orbit
of the Sun)
beta = 0 (by definition as the Sun's orbit defines the
ecliptic plane. This results in a simplification
of the formulas below)
5. Find the obliquity of the ecliptic plane (epsilon)
epsilon = 23.439 - 0.0000004 * d
= 23.439351
6. Find the Right Ascension (alpha) and Declination (delta) of
the Sun at the day d is given by the formulas below
Y = cos(epsilon) * sin(lambda)
X = cos(lambda)
a = arctan(Y/X)
If X < 0 then alpha = a + 180
If Y < 0 and X > 0 then alpha = a + 360
else alpha = a
Y = 0.6489924
X = -0.7068507
a = -42.556485
alpha = -42.556485 + 180 = 137.44352 (degrees)
delta = arcsin(sin(epsilon)*sin(lambda))
= 16.342193 degrees
7. Find the Aziumuth (Z) and altitude (A) of the Sun, we first
find the local sidereal time, then the hour angle of the Sun.
Long is the longitude (west negative) of the place, and lat is
the latitude. LST and RA are taken in degrees. Note that the
calculation of LST needs about 12 places of precision.
LST = 280.46061837 + 360.98564736629 * d + long
bring into range 0 to 360
ha = LST - RA
bring into range 0 to 360
d = -877.04167
LST = -316320.911064 degs
= 119.088936 degs
ha = 119.088936 - 137.44352degs
= -18.354584 degs
8. Now we can use the formulas below to convert from ha and delta
to Z and A;
sin(A) = sin(delta) * sin(lat) + cos(delta) * cos(lat) * cos(ha)
from which A follows by taking the arcsin
Z is given from a tangent formula...
Y = -cos(delta) * cos(lat) * sin(ha)
X = sin(delta) - sin(lat) * sin(alt)
Z' = arctan(y/x)
If X < 0 then Z = Z' + 180
If Y < 0 and X > 0 then Z = Z' + 360
else Z = Z'
We have...
sin(delta) = 0.2813734
cos(delta) = 0.9595983
Sin(lat) = 0.7933533
cos(lat) = 0.6087614
sin(ha) = -0.3148968
cos(ha) = 0.9491259
And so...
sin(alt) = 0.2813734 * 0.7933533 + 0.9595983 * 0.6087614 * 0.9491259
= 0.7776760
A = 51.048280 = 51.05 degrees (SkyMap gives 51.061)
Y = -0.9595983 * 0.6087614 * -0.3148968 = 0.1839521
X = 0.2813734 - 0.7933533 * 0.7776760 = -0.3355984
Z' = atan(y/x) = -28.7285454 deg
As X < 0, we add 180 degrees to Z'
Z = 151.2714587 = 151.27 degrees (SkyMap gives 151.279)
Final result
Right ascension is usually given in hours of time, and both
figures need to be rounded to a sensible number of decimal
places.
alpha = 9.163 hrs or 9h 09m 46s
delta = +16.34 degrees or +16d 20' 32"
The Interactive Computer Ephemeris gives
alpha = 9h 09m 45.347s and
delta = +16d 20' 30.89"
'*********************************************************
' This program will calculate the position of the Sun
' using a low precision method found on page C24 of the
' 1996 Astronomical Almanac.
'
' The method is good to 0.01 degrees in the sky over the
' period 1950 to 2050.
'
' QBASIC program by Keith Burnett (kburnett@geocity.com)
'
'
' Work in double precision and define some constants
'
DEFDBL A-Z
pr1$ = "\ \#####.##"
pr2$ = "\ \#####.#####"
pr3$ = "\ \#####.###"
pi = 4 * ATN(1)
tpi = 2 * pi
twopi = tpi
degs = 180 / pi
rads = pi / 180
'
' Get the days to J2000
' h is UT in decimal hours
' 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 and a modified inverse
' tangent function
'
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 atn2 function below returns an angle in the range 0 to two pi
' depending on the signs of x and y.
'
DEF FNatn2 (y, x)
a = ATN(y / x)
IF x < 0 THEN a = a + pi
IF (y < 0) AND (x > 0) THEN a = a + tpi
FNatn2 = a
END DEF
'
' the function below returns the true integer part,
' even for negative numbers
'
DEF FNipart (x) = SGN(x) * INT(ABS(x))
'
' 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
'
' Find the ecliptic longitude of the Sun
'
DEF FNsun (d)
'
' mean longitude of the Sun
'
L = FNrange(280.461 * rads + .9856474# * rads * d)
'
' mean anomaly of the Sun
'
g = FNrange(357.528 * rads + .9856003# * rads * d)
'
' Ecliptic longitude of the Sun
'
FNsun = FNrange(L + 1.915 * rads * SIN(g) + .02 * rads * SIN(2 * g))
'
' Ecliptic latitude is assumed to be zero by definition
'
END DEF
'
'
'
CLS
'
' get the date and time from the user
'
INPUT " year : ", y
INPUT " month : ", m
INPUT " day : ", day
INPUT "hour UT : ", h
INPUT " minute : ", mins
INPUT " lat : ", glat
INPUT " long : ", glong
glat = glat * rads
glong = glong * rads
h = h + mins / 60
d = FNday(y, m, day, h)
'
' Use FNsun to find the ecliptic longitude of the
' Sun
'
lambda = FNsun(d)
'
' Obliquity of the ecliptic
'
obliq = 23.439 * rads - .0000004# * rads * d
'
' Find the RA and DEC of the Sun
'
alpha = FNatn2(COS(obliq) * SIN(lambda), COS(lambda))
delta = FNasin(SIN(obliq) * SIN(lambda))
'
' Find the Earth - Sun distance
'
r = 1.00014 - .01671 * COS(g) - .00014 * COS(2 * g)
'
' Find the Equation of Time
'
equation = (L - alpha) * degs * 4
'
' find the Alt and Az of the Sun for a given position
' on Earth
'
' hour angle of Sun
LMST = FNrange((280.46061837# + 360.98564736629# * d) * rads + glong)
hasun = FNrange(LMST - alpha)
'
' conversion from hour angle and dec to Alt Az
sinalt = SIN(delta) * SIN(glat) + COS(delta) * COS(glat) * COS(hasun)
altsun = FNasin(sinalt)
y = -COS(delta) * COS(glat) * SIN(hasun)
x = SIN(delta) - SIN(glat) * sinalt
azsun = FNatn2(y, x)
'
' print results in decimal form
'
PRINT
PRINT "Position of Sun"
PRINT "==============="
PRINT
PRINT USING pr2$; " days : "; d
PRINT USING pr1$; "longitude : "; lambda * degs
PRINT USING pr3$; " RA : "; alpha * degs / 15
PRINT USING pr1$; " DEC : "; delta * degs
PRINT USING pr2$; " distance : "; r
PRINT USING pr1$; " eq time : "; equation
PRINT USING pr1$; " LST : "; FNrange(LMST) * degs
PRINT USING pr1$; " azimuth : "; azsun * degs
PRINT USING pr1$; " altitude : "; altsun * degs
END
'*********************************************************
Below is the output from the old program when run for 11:00 UT on 1997
August 7. year : 1997
month : 8
day : 7
hour UT : 11
minute : 0
Position of Sun
===============
days : -877.04167
longitude : 134.98
RA : 9.163
DEC : 16.34
distance : 1.01408
eq time : -5.75
Below is the output from the program including Altitude and Azimuth in
Chicago, run for 15:00 UT on 2001 March 4th (09:00h Chicago time), compared with
the altitude and azimuth calculated from Chris Marriott's SkyMap Pro 6 demo. year : 2001
month : 3
day : 4
hour UT : 15
minute : 30
lat : 41.87
long : -87.64
Position of Sun
===============
days : 428.14583
longitude : 344.13
RA : 23.025
DEC : -6.24
distance : 0.99173
eq time : -11.68
azimuth : 134.56
altitude : 30.68
Skymap 6
--------
Right ascension: 23h 1m 29.93s = 23.025
Declination: -6° 14' 58.0" = -6.249
Altitude: 30° 42' 20" = 30.706
Azimuth: 134° 33' 41" = 134.561
The errors here are about 1.5 arcmin for altitude and much less for
Azimuth. Errors are higher than for the RA and Dec and this may be due to the
failure to allow for the TDT-UT time difference in the method given here. To put
this all in perspective, the Sun moves about 0.25 of a degree in the sky between
0900 and 0901 that morning!
You can implement the formulas above quite easily on a spreadsheet. You can download a spreadsheet in MS Excel 97 format as a small ZIP file. This spreadsheet can be loaded directly into Sun Star Office 5.2 onwards, as the screenshot shows. The main features of the spreadsheet implementation are
radians() function, degrees() function
and mod(x, 360) functions to convert angular arguments into
radian measure between 0 and 2pi If the year, month day, UT hour and minute are in cells B6 to
B10, then the formula below will return the days since J2000.0 for
any date since March 1900 until Dec 2099.
=367*B6-INT(7*(B6+INT((B7+9)/12))/4)+INT(275*B7/9)+B8+(B9+B10/60)/24-730531.5
If the days since J2000.0 are in cell B24, and the longitude in
degrees (West negative) is in cell B11, then the mean local
sidereal time is given in degrees as
= MOD(280.46061837 + 360.98564736629 * B24 + B11, 360)
And you can convert that to radians as
=RADIANS(B25)You could concatenate the formulas so that cells always gave degrees in the results, like (say) this
=SIN(RADIANS(MOD(280.46061837 + 360.98564736629 * B24 + B11, 360)))
to give the sine of the local sidereal time, but I find it easier to keep column B for degrees, Column C for radians, D for sines, E for cosines and so on. Each row is a new step in the calculation.
Having said that, the Sun has ecliptic longitude of zero to an accuracy of arcseconds (the ecliptic plane being defined by the plane of the Sun-Earth orbit) so a one cell formula for the ecliptic longitude of the Sun in degrees might be useful.
=MOD(280.461 + 0.9856474 * B24, 360) +1.915 * SIN(RADIANS(MOD(357.528 + 0.9856003 * B24, 360))) +0.02 * SIN(2 * RADIANS(MOD(357.528 + 0.9856003 * B24, 360)))Where cell
B24 must contain the days since J2000.0. This sort
of thing can get inefficient, notice that the mean anomaly is being calculated
twice for the same instant in the one cell formula above. With modern
spreadsheets and computers, this minor duplication might matter less than the
presentation advantage gained by using a more compact layout.
I modified the QBASIC program above to produce a file of
positions for days from -20,000 to +20,000 - a 106 year period centred on
J2000.0. The RA and DEC figures were rounded to 4 places of decimals in this
file. I used Planeph to generate a similar file of positions for the Sun,
referred to the mean ecliptic and equinox of date. I then loaded both files into
a spreadsheet, and found the errors in seconds of time (RA) and arcseconds
(DEC). The maximum and minimum errors are shown in the table below for various
ranges of time about J2000.0
Sun error
RA sec DEC arcsec
Max within 3 year 0.6 8.9
Min within 3 year -2.1 -8.2
Max within 10 year 0.6 10.9
Min within 10 year -2.6 -12.5
Max within 50 year 1.0 16.8
Min within 50 year -2.9 -16.1
Error = C24 low precision method - Planeph
Note: Planeph was set to give output referred to mean
ecliptic and equinox of date.
[ Root ]
Last Modified 6th May 2001
Keith Burnett