Computer code to calculate Equation of Time and Solar Declination
The Equation of Time is the difference between the time that is shown by a sundial and the time that is shown by a clock that is set to local mean time. Essentially, the Equation of Time shows the longitude of the sun in the sky, relative to its mean position. The Solar Declination is the sun's latitude. Both of these quantities must be known for many purposes related to solar energy.
The following subroutine, here written in QBasic but easily translatable to other languages, calculates the Equation of Time and the Solar Declination on any day of the year. It is quite accurate. Its Root-Mean-Square error for the Equation of Time is only 3.7 seconds. The Declination is calculated with errors that are always small compared with the angular radius of the sun as seen from the earth (about 0.25 degrees).
Anyone who wants further information should follow the following link: [Link] The program that includes the routine is ETIMSDEC. There are instructions how to run it. The article titled "The Latitude and Longitude of the Sun", which I wrote several years ago, describes the astronomical logic behind the routine.
DOwenWilliams 10:19, 18 May 2011 (PDT) David Williams
FUNCTION ET.Dec (D, F%) STATIC
' Calculates equation of time, in minutes, or solar declination,
' in degrees, on day number D of year. (D = 0 on January 1.)
' F% selects function: True (non-zero) for Equation of Time,
' False (zero) for Declination.
' STATIC means variables are preserved between calls of function
IF PI = 0 THEN ' first call, initialize constants
PI = 4 * ATN(1)
W = 2 * PI / 365 ' earth's mean orbital angular speed in radians/day
DR = 180 / PI ' degree/radian factor
C = -23.45 / DR ' reverse angle of earth's axial tilt in radians
ST = SIN(C) ' sine of reverse tilt
CT = COS(C) ' cosine of reverse tilt
E2 = 2 * .0167 ' twice earth's orbital eccentricity
SP = 12 * W ' 12 days from December solstice to perihelion
D1 = -1 ' holds last D. Saves time if D repeated for both functions
END IF
IF D <> D1 THEN ' new value of D
A = W * (D + 10) ' Solstice 10 days before Jan 1
B = A + E2 * SIN(A - SP)
D1 = D
END IF
IF F% THEN ' equation of time calculation
C = (A - ATN(TAN(B) / CT)) / PI
ET.Dec = 720 * (C - INT(C + .5)) ' this is value of equation of time
' in 720 minutes, earth rotates PI radians relative to sun
ELSE ' declination calculation
C = ST * COS(B)
ET.Dec = ATN(C / SQR(1 - C * C)) * DR ' this is value of declination
' arcsine of C in degrees. ASN not directly available in QBasic
END IF
END FUNCTION