A few days ago I sent a fourier approximation for EOT including an equation
for calculating the day number of the year. Unfortunately I had hastily cut
and pasted the day-number formula from a web page that I had found. On
later checking I found that it didn't work properly - Oops.
Below I have repeated my fourier approximation with a corrected day number
calculator. The formula are all in Excel Visual Basic so you can easily put
them in an Excel spreadsheet.
John later asks "what and when is the greatest 24 hour difference on the
EOT. Using this fourier approximation and simple differences [ EOT(n) -
EOT(n-1) ] I get 29.78 seconds on 23.5UTC of December. This is consistent
with the more precise values quoted by Gianni, David, and Fer. Note that
the approximation is for 2002. However, the method of simple differences
introduce a 12 hour phase error so we would be better off producing the
differential dEOT/dt. As the fourier approximation is linear this can be
done with high school calculus. I've included the differential function
below. Luckily this produces the same numerical result (to two dec. places)
as before except the date is now 23.0UTC Dec (as expected from the phase
argument).
Hope this is useful to the group
Many Regards
Hank
Const pi = 3.1415926535
REM Equation of time in minutes
Function EOT(ByVal j)
Dim t
t = 2 * pi * (j - 1) / 365.25
EOT = -0.000062368659 _
+ 7.3636768 * Cos(t + 1.5076785) _
+ 9.9351605 * Cos(2 * t + 1.9332853) _
+ 0.32924369 * Cos(3 * t + 1.8934472) _
+ 0.21208209 * Cos(4 * t + 2.3032843)
End Function
REM Day number of the Year
Function dn(ByVal y&, ByVal m&, ByVal d&)
dn = g(y&, m&, d&) - g(y&, 1, 1)
End Function
Function g(ByVal y&, ByVal m&, ByVal d&) As Long
If m& < 3 Then 'months are January or February
m& = m& + 12 'January=13, February=14
y& = y& - 1 'of the previous year
End If 'months are January or February
g = d& + (153 * m& - 457) \ 5 + (36525 * y&) \ 100 - (y& * 0.01) \ 100 +
(y& \ 400) + 1721118.5
End Function
REM differential of EOT in minutes per day
Function dEOT(ByVal j)
Dim t
t = 2 * pi * (j - 1) / 365.25
dEOT = -7.3636768 * Sin(t + 1.5076785) _
- 2 * 9.9351605 * Sin(2 * t + 1.9332853) _
- 3 * 0.32924369 * Sin(3 * t + 1.8934472) _
- 4 * 0.21208209 * Sin(4 * t + 2.3032843)
dEOT = dEOT * 2 * pi / 365.25
End Function
Hank de Wit
Regional Computer Manager
Bureau of Meteorology, SA
ph: 08 8366 2674
Email: [EMAIL PROTECTED]
-
