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.
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.
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.
Position of the Sun at 11:00 UT on 1997 August 7th 1. Find the days before J2000.0 (d) from the table d = 11/24 + 212 + 7 - 1096.5 = -877.04167 2. Find the Mean Longitude (L) of the Sun 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 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 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 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 (email@example.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 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 ' ' 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 END '*********************************************************Below is the output from the 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
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.
Last Modified 7th August 1997