Depending on one's need for accuracy, there are any number of ways to calculate the Equation of Time. Wikipedia's EoT page and its Analemma page give a number of these. The simplest rely on the fact that both the Eccentricity and Obliquity effects are almost sine curves. Others are Fourier-based solutions. Finally there are those that use sound astronomical methods. Here you will find an example each method.
WHY BOTHER TO CALCULATE AT ALL...
On September one, trust the Sun.
Come Halloween, subtract sixteen.
On Christmas Day, you're OK.
For your Valentine true, add a dozen and two.
The mid of month four, add no more
At the mid of May, take four away.
On June fourteen, don't add a bean.
When August begins, add seven little mins.
The rest is easy: for any date, All you do is interpolate.
a poem by Tad Dunne
VERY SIMPLE METHOD
A 1970 letter from my Uncle John Wigham Richardson explaining how to calculate the Equation of Time
Better METHODS
Two routines are provided below, both require the calculation of the Number of Days since Noon UTC 1st January 2000. In both, the date/time are in local Standard Time : the Zone parameter conveys to UTC.
In Blue the functions below are in Visual Basic, which can be pasted directly into a Module in Excel's Visual Basic.
In Red, the functions are in Python
Function DY2k(Year, Month, Day, Hour, Zone)
If Month <= 2 Then
Year = Year - 1
Month = Month + 12
End If
AAAA = Int(Year / 100)
BBBB = 2 - AAAA + Int(AAAA / 4)
CCCC = Int(365.25 * Year)
DDDD = Int(30.6001 * (Month + 1))
DY2k = BBBB + CCCC + DDDD + Day + (Hour - Zone) / 24# - 730550.5
End Function
If Month <= 2 Then
Year = Year - 1
Month = Month + 12
End If
AAAA = Int(Year / 100)
BBBB = 2 - AAAA + Int(AAAA / 4)
CCCC = Int(365.25 * Year)
DDDD = Int(30.6001 * (Month + 1))
DY2k = BBBB + CCCC + DDDD + Day + (Hour - Zone) / 24# - 730550.5
End Function
def DY2k(Year,Month,Day,Hour,Zone):
if Month <= 2 :
Year = Year - 1
Month = Month + 12
AAAA = int(Year / 100)
BBBB = 2 - AAAA + int(AAAA / 4)
CCCC = int(365.25 * Year)
DDDD = int(30.6001 * (Month + 1))
JD = BBBB + CCCC + DDDD + Day + (Hour-Zone)/24.- 730550.5
return JD
if Month <= 2 :
Year = Year - 1
Month = Month + 12
AAAA = int(Year / 100)
BBBB = 2 - AAAA + int(AAAA / 4)
CCCC = int(365.25 * Year)
DDDD = int(30.6001 * (Month + 1))
JD = BBBB + CCCC + DDDD + Day + (Hour-Zone)/24.- 730550.5
return JD
ASTRONOMICAL ALMANAC METHOD
This method is adapted from Section C5 of the Astronomical Almanac and EoT is accurate to +/- 3.5 seconds during this Century.
Function EoT(Year, Month, Day, Hour, Zone)
Days_Since_Epoch = DY2k(Year, Month, Day, Hour, Zone)
Mean_Long_Sun_d = 280.46 + 0.9856474 * Days_Since_Epoch
Mean_Long_Sun_d = Mean_Long_Sun_d - Fix(Mean_Long_Sun_d / 360#) * 360#
Mean_Anomaly_d = 357.528 + 0.9856003 * Days_Since_Epoch
Mean_Anomaly_d = Mean_Anomaly_d - Fix(Mean_Anomaly_d / 360) * 360#
Mean_Anomaly_r = Mean_Anomaly_d * WorksheetFunction.Pi / 180#
Ecliptic_Long_d = Mean_Long_Sun_d + 1.915 * Sin(Mean_Anomaly_r) + 0.02 * Sin(2 * Mean_Anomaly_r)
Ecliptic_Long_r = Ecliptic_Long_d * WorksheetFunction.Pi / 180#
Obliquity_d = 23.439 - 4E-07 * Days_Since_Epoch
Obliquity_r = Obliquity_d * WorksheetFunction.Pi / 180#
Right_Ascension_r = WorksheetFunction.Atan2(Cos(Ecliptic_Long_r), Cos(Obliquity_r) * Sin(Ecliptic_Long_r))
Right_Ascension_d = Right_Ascension_r * 180# / WorksheetFunction.Pi
If Right_Ascension_d < 0 Then Right_Ascension_d = Right_Ascension_d + 360#
Right_Ascension_h = Right_Ascension_d / 15#
Declination_r = WorksheetFunction.Asin(Sin(Obliquity_r) * Sin(Ecliptic_Long_r))
Declination_d = Declination_r * 180# / WorksheetFunction.Pi
EoT_d = Mean_Long_Sun_d - Right_Ascension_d
EoT_m = 4 * EoT_d
If EoT_m > 50 Then EoT_m = EoT_m - 1440 'Astronomical EoT
EoT = -EoT_m 'Gnomoc EoT
End Function
Days_Since_Epoch = DY2k(Year, Month, Day, Hour, Zone)
Mean_Long_Sun_d = 280.46 + 0.9856474 * Days_Since_Epoch
Mean_Long_Sun_d = Mean_Long_Sun_d - Fix(Mean_Long_Sun_d / 360#) * 360#
Mean_Anomaly_d = 357.528 + 0.9856003 * Days_Since_Epoch
Mean_Anomaly_d = Mean_Anomaly_d - Fix(Mean_Anomaly_d / 360) * 360#
Mean_Anomaly_r = Mean_Anomaly_d * WorksheetFunction.Pi / 180#
Ecliptic_Long_d = Mean_Long_Sun_d + 1.915 * Sin(Mean_Anomaly_r) + 0.02 * Sin(2 * Mean_Anomaly_r)
Ecliptic_Long_r = Ecliptic_Long_d * WorksheetFunction.Pi / 180#
Obliquity_d = 23.439 - 4E-07 * Days_Since_Epoch
Obliquity_r = Obliquity_d * WorksheetFunction.Pi / 180#
Right_Ascension_r = WorksheetFunction.Atan2(Cos(Ecliptic_Long_r), Cos(Obliquity_r) * Sin(Ecliptic_Long_r))
Right_Ascension_d = Right_Ascension_r * 180# / WorksheetFunction.Pi
If Right_Ascension_d < 0 Then Right_Ascension_d = Right_Ascension_d + 360#
Right_Ascension_h = Right_Ascension_d / 15#
Declination_r = WorksheetFunction.Asin(Sin(Obliquity_r) * Sin(Ecliptic_Long_r))
Declination_d = Declination_r * 180# / WorksheetFunction.Pi
EoT_d = Mean_Long_Sun_d - Right_Ascension_d
EoT_m = 4 * EoT_d
If EoT_m > 50 Then EoT_m = EoT_m - 1440 'Astronomical EoT
EoT = -EoT_m 'Gnomoc EoT
End Function
def EoT(Year,Month,Day,Hour,Minute,Second,Zone,DST):
Days_since_Epoch = DY2k(Year,Month,Day,Zone)
Mean_Long_Sun_d = (280.46 + 0.9856474 * Days_since_Epoch) % 360.
Mean_Anomaly_d = (357.528 + 0.9856003 * Days_since_Epoch) % 360.
Mean_Anomaly_r = radians(Mean_Anomaly_d)
Ecliptic_Long_d = Mean_Long_Sun_d + 1.915 *sin(Mean_Anomaly_r) + 0.020 * sin(2*Mean_Anomaly_r)
Ecliptic_Long_r = radians(Ecliptic_Long_d)
Obliquity_d = 23.439 - 0.0000004 * Days_since_Epoch
Obliquity_r = radians(Obliquity_d)
Days_since_Epoch = DY2k(Year,Month,Day,Zone)
Mean_Long_Sun_d = (280.46 + 0.9856474 * Days_since_Epoch) % 360.
Mean_Anomaly_d = (357.528 + 0.9856003 * Days_since_Epoch) % 360.
Mean_Anomaly_r = radians(Mean_Anomaly_d)
Ecliptic_Long_d = Mean_Long_Sun_d + 1.915 *sin(Mean_Anomaly_r) + 0.020 * sin(2*Mean_Anomaly_r)
Ecliptic_Long_r = radians(Ecliptic_Long_d)
Obliquity_d = 23.439 - 0.0000004 * Days_since_Epoch
Obliquity_r = radians(Obliquity_d)
Right_Ascension_r = atan2(cos(Obliquity_r) * sin(Ecliptic_Long_r),cos(Ecliptic_Long_r))
Right_Ascension_d = (degrees(Right_Ascension_r)) % 360.
Right_Ascension_h = Right_Ascension_d/15.
Declination_r = asin(sin(Obliquity_r) * sin(Ecliptic_Long_r))
Declination_d = degrees(Declination_r)
EoT_d = Mean_Long_Sun_d - Right_Ascension_d
if EoT_d > 50 : EoT_d -= 360
EoT_m = -EoT_d * 4 # Gnomonical EoT
Right_Ascension_d = (degrees(Right_Ascension_r)) % 360.
Right_Ascension_h = Right_Ascension_d/15.
Declination_r = asin(sin(Obliquity_r) * sin(Ecliptic_Long_r))
Declination_d = degrees(Declination_r)
EoT_d = Mean_Long_Sun_d - Right_Ascension_d
if EoT_d > 50 : EoT_d -= 360
EoT_m = -EoT_d * 4 # Gnomonical EoT
return EoT_m
FOURIER METHOD
This method provides the quickest way to find the Equation of Time. This was developed by many EoT calculations (every 6 hours every day over this Century) from NASA/JPL Horizons application (see below for details of this application). Amongst the many variables that Horizons can calculate is the Hour Angle of the Sun, from which it is easy to deduce the EoT. The EoT values thus provided were subjected to rigorous Fourier analysis. This gives results that are within +/- 13 seconds of the value predicted by Horizons over this Century.
Function EoT_F(Year, Month, Day, Hour, Zone)
Days_Since_Epoch = D2000(Year, Month, Day, Hour, Zone)
Days_Since_Epoch = D2000(Year, Month, Day, Hour, Zone)
Cycle = 4 * Days_Since_Epoch
Cycle = Cycle - Fix(Cycle / 1461#) * 1461#
Cycle = Cycle - Fix(Cycle / 1461#) * 1461#
Theta = Cycle * 0.004301
EoT1 = 7.353 * Sin(1 * Theta + 6.209)
EoT2 = 9.927 * Sin(2 * Theta + 0.37)
EoT3 = 0.337 * Sin(3 * Theta + 0.304)
EoT4 = 0.232 * Sin(4 * Theta + 0.715)
EoT_F = 0.019 + EoT1 + EoT2 + EoT3 + EoT4
End Function
EoT1 = 7.353 * Sin(1 * Theta + 6.209)
EoT2 = 9.927 * Sin(2 * Theta + 0.37)
EoT3 = 0.337 * Sin(3 * Theta + 0.304)
EoT4 = 0.232 * Sin(4 * Theta + 0.715)
EoT_F = 0.019 + EoT1 + EoT2 + EoT3 + EoT4
End Function
def EoT_Fourier(Year,Month,Day,Hour,Minute,Second,Zone,DST):
UTC_hrs = Hour + Minute/60. + Second/3600. - Zone - DST
UTC_hrs = Hour + Minute/60. + Second/3600. - Zone - DST
JD_Date = Julian(Year,Month,Day,0,0,0)
JD_2000 = Julian(2000,1,1,12,0,0)
JD_days = JD_Date - JD_2000 + UTC_hrs / 24.
JD_2000 = Julian(2000,1,1,12,0,0)
JD_days = JD_Date - JD_2000 + UTC_hrs / 24.
Cycle = (4 * JD_days) % 1461.
Theta = Cycle * 0.004301
EoT1 = 7.353 * sin(1 * Theta + 6.209)
EoT2 = 9.927 * sin(2 * Theta + 0.37)
EoT3 = 0.337 * sin(3 * Theta + 0.304)
EoT4 = 0.232 * sin(4 * Theta + 0.715)
EoT = 0.019 + EoT1 + EoT2 + EoT3 + EoT4
return EoT
Theta = Cycle * 0.004301
EoT1 = 7.353 * sin(1 * Theta + 6.209)
EoT2 = 9.927 * sin(2 * Theta + 0.37)
EoT3 = 0.337 * sin(3 * Theta + 0.304)
EoT4 = 0.232 * sin(4 * Theta + 0.715)
EoT = 0.019 + EoT1 + EoT2 + EoT3 + EoT4
return EoT
FURTHER DETAILS ON CALCULATING THE ELLIPTICITY EFFECT
MEEUS' ALGORITHMS
The most complete 'non-professional' algorithms for calculation of the Equation of Time are provided in by Jan Meeus - Astronomical Algorithms (1998), 2nd ed - ISBN 0-943396-61-1.
You can view the author's Javascript implementation of Meeus' algorithms for the Sun here. If you would like to use these algorithms, contact the author, who would be happy to send you a text file.
OTHER SOURCES
For more information of the astronomical background of the Equation of Time, see the article below which was published in NASS Compendium : Vol 25 Nos 3 & 4, Sept & Dec 2018. (Including some corrections for the published text)
MICA
From the US Naval Observatory, this is a relatively cheap program that provides everything that an serious amateur astronomer, gnomonist or navigator might need. N,B. The Mac version of MICA does NOT work under recent Apple Mac OSX.
HORIZONS
This is the program, in whose user guide, it states that the user should consult the web-master if using the program for manned planetary landings. It is east to use and quite straight-forward. It is lightening fast, free and astronomically the best that there is.
The above images gives an indication of how input to Horizons is set up. Note the time must be UTC. The quantities, under Table Settings, will output Right Ascension & Declination, Azimuth & Altitude, and Local Solar Time (subtract the UTC value to get the gnomonical convention topographical Equation of Time plus the longitude correction)