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 (firstname.lastname@example.org) ' ' ' 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.75Below 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.561The 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
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
B10, then the formula below will
return the days since J2000.0 for any date since March 1900 until
If the days since J2000.0 are in cell
the longitude in degrees (West negative) is in cell
B11, then the mean local sidereal time is given in
= 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
B24must 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