User:Airman72/sandbox

Calculating the equation of time

As noted in the introduction to this article, for many purposes the equation of time is obtained by looking it up in a published table of values or on a graph. For dates in the past, such tables are produced either from measurements done at the time, or by calculations. Of course for dates in the future, tables can only be prepared from calculations. Also, in devices such as computer-controlled heliostats, the computer is often programmed to calculate the equation of time whenever it is needed, instead of looking it up. Two types of calculation are currently in use, numerical and analytical. The former are based on numerical integration of the differential equations of motion that include all significant gravitational interactions and relativistic effects. The results are accurate to better than 1 second of time and form the basis for Almanac data. The later are based on the known solution of the equations of motion that include only the gravitational interaction between the Sun and Earth. These produce results that are not as accurate as the former, but they can be understood more easily and the equations illustrate the functional dependence on the governing parameters. Their accuracy can be improved by including small corrections, but in the process concurrent increases in complexity are incurred.

The following discussion describes a simple yet reasonably accurate (agreeing with Almanac data to within 3 seconds of time over a wide range of years) algorithm for calculating the equation of time that is well known to astronomers. It also shows how to obtain a simple approximate formula (accurate to within 1 minute of time over a large time interval), that can be easily evaluated with a calculator and provides the simple explanation of the phenomonon that was utilized previously in this article.

Mathematical description

The precise definition of the Equation of Time is

${\mbox{EOT = GHA - GMHA}}$

The quantities occurring in this equation are

• GHA is the Greenwich hour angle of the apparent (actual) Sun,
• GMHA = UT - Offset, is the Greenwich hour angle of the mean (fictitious) Sun,
• EOT is the time difference between apparent solar time and mean solar time.

Here time and angle are quantities that are related by conversion factors such as; 2π radian = 360° = 1 day = 24 hour. The difference, EOT, is measurable since GHA is an angle that can be measured and Universal Time, UT, is a scale for the measurement of time. The Offset = π = 180° = 12 hour is needed because UT is zero at mean midnight while GMHA = 0 at mean noon. Universal Time varies from 0 to 24 hours and is discontinuous at mean midnight, another quantity day number N, an integer, is required in order to form the continuous quantity time, t = (N + UT/24 hr) day. Both GHA and GMHA, like all physical angles, have a mathematical, but not a physical discontinuity at their respective (apparent and mean) noon. Despite the mathematical discontinuities of its components, EOT is defined as a continuous function by adding (or subtracting) 24 hours in the small time interval between the discontinuities in GHA and GMHA.

According to the definitions of the angles on the celestial sphere GHA = GAST - α (see hour angle) where

• GAST is the Greenwich apparent sidereal time (the angle between the apparent vernal equinox and the meridian in the plane of the equator); it is a known function of UT,
• α is the right ascension of the apparent Sun (the angle between the apparent vernal equinox and the actual Sun in the plane of the equator).

On substituting into the equation of time, it is

${\mbox{EOT = GAST}}-\alpha -{\mbox{UT}}+{\mbox{Offset}}$

Like the formula for GHA above, one can write GMHA = GAST - α_M, where the last term is the right ascension of the mean Sun. The equation is often written in these terms as

${\mbox{EOT}}=\alpha _{M}-\alpha$

where α_M = GAST - UT + Offset. In this formulation a measurement or calculation of EOT at a certain value of time depends on a measurement or calculation of α at that time. Both α and α_M vary from 0 to 24 hours during the course of a year. The former has a discontinuity at a time that depends on the value of UT, while the later has its at a slightly later time. As a consequence, when calculated this way EOT has two, artificial, discontinuities. They can both be removed by subtracting 24 hours from the value of EOT in the small time interval after the discontinuity in α and before the one in α_M. The resulting EOT is a continuous function of time.

Another definition, denoted E to distinguish it from EOT, is

${\mbox{E = GMST}}-\alpha -{\mbox{UT}}+{\mbox{Offset}}$

Here GMST = GAST - eqeq, is the Greenwich mean sidereal time (the angle between the mean vernal equinox and the mean Sun in the plane of the equator). Therefore GMST is an approximation to GAST (and E is an approximation to EOT); eqeq is called the equation of the equinoxes and is due to the wobbling, or nutation of the Earth's axis of rotation about its precessional motion. Since the amplitude of the nutational motion is only about 1.2 sec of time (18 arcsec of longitude) the difference between EOT and E can be ignored unless one is interested in subsecond accuracy.

A third definition, denoted Δt to distinguish it from EOT and E, and now called the Equation of Ephemeris Time (prior to the distinction that is now made between EOT, E, and Δt the latter was known as the Equation of Time) is

$\Delta t=\Lambda -\alpha$

here Λ is the ecliptic longitude of the mean Sun (the angle from the mean vernal equinox to the mean Sun in the plane of the ecliptic).

The difference Λ - [GMST - UT + Offset], is 1.3 seconds of time from 1960 to 2040. Therefore over this restricted range of years Δt is an approximation to EOT whose error is in the range 0.1 to 2.5 sec depending on the longitude correction in the equation of the equinoxes; for many purposes, for example correcting a sundial, this accuracy is more than good enough.

Right ascension calculation

The right ascension, hence the equation of time, can be calculated based on Newton's 2 body theory of celestial motion, in which the bodies (earth and sun) describe elliptical orbits about their common mass center. Using this theory the equation of time becomes

$\Delta t=M+\lambda _{p}-\alpha$

where the new angles that appear are

• M = 2π(t - t_p) /t_Y, is the mean anomaly, the angle from the periapsis of the elliptical orbit to the mean Sun; its range is from 0 to 2π as t increases from t_p to t_p + t_Y,
• t_Y = 365.2596358 day is the length of time in an anomalistic year; the time interval between two successive passages of the periapsis,
• λ_p = Λ-M, is the ecliptic longitude of the periapsis,
• t is dynamical time, the independent variable in the theory. Here it is taken to be identical with the continuous time based on UT (see above), but in more precise calculations (of E or EOT) the small difference between them must be accounted for as well as the distinction between UT1 and UTC.

To complete the calculation three additional angles are required; they are

• E the Sun's eccentric anomaly (note that this is different than E),
• ν the Sun's true anomaly,
• λ = ν + λ_p the Sun's true longitude on the ecliptic.

The celestial sphere and the Sun's elliptical orbit as seen by a geocentric observer looking normal to the ecliptic showing the 6 angles (M, λ_p, α, ν, λ, E) needed for the calculation of the equation of time. For the sake of clarity the drawings are not to scale.

All these angles are shown in the figure on the right, which shows the celestial sphere and the Sun's elliptical orbit seen from the Earth (the same as the Earth's orbit seen from the Sun). In this figure ε is the obliquity, while e = [1 − (b/a)²]^1/2 is the eccentricity of the ellipse.

Now given a value of 0≤M≤2π, one can calculate α(M) by means of the following, well known, procedure:

First, given M, calculate E from Kepler's equation

$M=E-e\sin E$

Although this equation cannot be solved exactly in closed form, values of E(M) can be obtained from infinite (power or trigonometric) series, graphical, or numerical methods. Alternatively, note that for e = 0, E = M, and by iteration, E ~ M + e sin M. This approximation can be improved, for small e, by iterating again, E ~ M + e sin M + (1/2) e² sin 2M, and continued iteration produces successively higher order terms of the power series expansion in e. For small values of e (much less than 1) two or three terms of the series give a good approximation for E; the smaller e, the better the approximation.

Next, knowing E, calculate the true anomaly ν from an elliptical orbit relation

$\nu =2\tan ^{-1}\left[{\sqrt {\frac {1+e}{1-e}}}\tan {\frac {E}{2}}\right]$

The correct branch of the multiple valued function tan⁻¹x to use is the one that makes ν a continuous function of E(M) starting from ν(E=0) = 0. Thus for 0≤ E < π use tan⁻¹x = Tan⁻¹x, and for π < E ≤ 2π use tan⁻¹x = Tan⁻¹x + π. At the specific value E = π for which the argument of tan is infinite, use ν = E. Here Tan⁻¹x is the principal branch, |Tan⁻¹x| < π/2; the function that is returned by calculators and computer applications. Alternatively, this function can be expressed in terms of its Taylor series in e, the first three terms of which are, ν ~ E + e sin E + (1/4) e² sin 2E. For small e this approximation (or even just the first two terms) is a good one. Combining the approximation for E(M) with this one for ν(E) produces

$\nu \sim M+2e\sin M+(5/4)e^{2}\sin 2M$

The relation ν(M) is called the Equation of the center; the expression written here is a second order approximation in e. For the small value of e that characterises the Earth's orbit this gives a very good approximation for ν(M).

Next knowing ν calculate λ from its definition above

$\lambda =\nu +\lambda _{p}$

The value of λ varies non-linearly with M because the orbit is elliptical and not circular. From the approximation for ν, λ ~ M + λ_p + 2e sin M +(5/4)e² sin 2M.

Finally, knowing λ calculate α from a relation for the right triangle on the celestial sphere shown above

$\alpha =\tan ^{-1}[\cos \varepsilon \,\tan \lambda ]$

Note that the quadrant of α is the same as that of λ, therefore reduce λ to the range 0 to 2π and write α = Tan⁻¹[cos ε tan λ] + kπ, where k is 0 if λ is in quadrant 1, it is 1 if λ is in quadrants 2 or 3 and it is 2 if λ is in quadrant 4. For the values at which tan is infinite, α = λ.

Although approximate values for α can be obtained from truncated Taylor series like those for ν, it is more efficatious to use the equation

$\alpha =\lambda -\sin ^{-1}[y\sin(\alpha +\lambda )]$

where y = tan²(ε/2). Note that for ε = y = 0, α = λ and iterating twice, α ~ λ - y sin 2λ + (1/2)y² sin 4λ.

Equation of time

The Equation of Time is obtained by substituting the result of the right ascension calculation into an equation of time formula. Here Δt(M) = M + λ_p - α[λ(M)] is used; in part because small corrections (of the order of a second of time), that would justify using E, are not included, and in part because the goal is to obtain a simple analytical expression. Using two term approximations for λ(M) and α(λ), allows Δt to be written as an explicit expression of two terms, which is designated Δt_ey because it is a first order approximation in e and in y.

$\Delta t_{ey}=-2e\sin M+y\sin(2M+2\lambda _{p})=[-7.659\sin M+9.863\sin(2M+3.5932)]{\mbox{min}}$

This equation was first derived by Milne, who wrote it in terms of Λ = M + λ_p. The numerical values written here result from using the orbital parameter values, e = 0.016709, ε = 23.4393° = 0.409093, and λ_p = 282.9381° = 4.938201 that correspond to the epoch 1 January 2000 at 12 noon. When evaluating the numerical expression for Δt_ey as given above, a calculator must be in radian mode to obtain correct values because the value of 2λ_p - 2π in the argument of the second term is written there in radians. Higher order approximations can also be written, but they necessarily have more terms. For example, the second order approximation in both e and y consists of five terms

$\Delta t_{e^{2}y^{2}}=\Delta t_{ey}-(5/4)e^{2}\sin 2M+ey\sin M\cos(2M+2\lambda _{p})-(1/2)y^{2}\sin(4M+4\lambda _{p})$

This approximation has the potential for high accuracy, however in order to achieve it over a wide range of years, the parameters e, ε, and λ_p must be allowed to vary with time. This creates additional calculational complications. Other approximations have been proposed, for example, Δt_e which uses the first order equation of the center but no other approximation to determine α, and Δt_e² which uses the second order equation of the center.

The time variable, M, can be written either in terms of, n, the number of days past perihelion, or, D, the number of days past a specific date and time (epoch)

$M={\frac {2\pi }{t_{Y}}}n{\mbox{ day}}=M_{D}+{\frac {2\pi }{t_{Y}}}D{\mbox{ day}}=6.24004077+0.01720197D$

Here M_D is the value of M at the chosen date and time. For the values given here, in radians, M_D is that measured for the actual Sun at the epoch, 1 January 2000 at 12:00 noon, and D is the number of days past that epoch. At periapsis M = 2π, so solving gives D = D_p = 2.508109. This puts the periapsis on 4 Jan 2000 at 11 min and 41 seconds past midnight while the actual periapsis is, according to results from the Multiyear Interactive Computer Almanac (abbreviated as MICA), on 3 Jan 2000 at 5 hr, 17 min and 30 seconds past midnight. This large discrepancy happens because the difference between the orbital radius at the two locations is only 1 part in a million; in other words, radius is a very weak function of time near periapsis. As a practical matter this means that one cannot get a highly accurate result for the equation of time by using n and adding the actual periapsis date for a given year. However, high accuracy can be achieved by using the formulation in terms of D.

Curves of Δt and Δt_ey along with symbols locating the daily values at noon (at 10 day intervals) obtained from the Multiyear Interactive Computer Almanac vs d for the year 2000.

When D > D_p, M is greater than 2π and one must subtract a multiple of 2π (that depends on the year) from it to bring it into the range 0 to 2π. Likewise for years prior to 2000 one must add multiples of 2π. For example, for the year 2010, D varies from 3653 on 1 January at noon to 4017 on 31 December at noon, the corresponding M values are 69.0789468 and 75.3404748 and are reduced to the range 0 to 2π by subtracting 10 and 11 times 2π respectively. One can always write D = n_Y + d, where n_Y is the number of days from the epoch to noon on 1 January of the desired year, and 0≤d≤ 364 (365 if the calculation is for a leap year).

The result of the computations is usually given as either a set of tabular values, or a graph of the equation of time as a function of d. A comparison of plots of Δt, Δt_ey, and results from MICA all for the year 2000 is shown in the figure on the right. The plot of Δt_ey is seen to be close to the results produced by MICA, the absolute error, Err = |Δt_ey − MICA2000|, is less than 1 minute of time throughout the year; its largest value is 43.2 sec and occurs on day 276 (October 3). The plot of Δt is indistinguishable from the results of MICA, the largest absolute error between the two is 2.46 sec on day 324 (November 20).

Secular effects

The difference between the MICA and Δt results was checked every 5 years over the range from 1960 to 2040. In every instance the maximum absolute error was less than 3 seconds of time, the largest difference, 2.91 seconds occurred on 22 May 1965 (day 141). However, in order to achieve this level of accuracy over this range of years it is necessary to account for the secular change in the orbital parameters with time. The equations that describe this variation are

$e=1.6709\cdot 10^{-2}-4.193\cdot 10^{-5}\left({\frac {D}{36525}}\right)-1.26\cdot 10^{-7}\left({\frac {D}{36525}}\right)^{2}$

$\varepsilon =\left[23.4393-0.013\left({\frac {D}{36525}}\right)-2\cdot 10^{-7}\left({\frac {D}{36525}}\right)^{2}+5\cdot 10^{-7}\left({\frac {D}{36525}}\right)^{3}\right]{\mbox{ deg}}$

$\lambda _{p}=\left[282.93807+1.7195\left({\frac {D}{36525}}\right)+3.025\cdot 10^{-4}\left({\frac {D}{36525}}\right)^{2}\right]{\mbox{ deg}}$

According to these relations, in 100 years (D = 36525), λ_p increases by about 1/2 percent (1.7 degrees), e decreases by about 1/4 percent, and ε decreases by about 1/20 percent.

As a result the number of calculations required for any of the higher order approximations of the equation of time requires a computer to complete them, if one wants to achieve their inherent accuracy over a wide range of time. In this event it is no more difficult to evaluate Δt using a computer than any of its approximations.

In all this note that Δt_ey as written above is easy to evaluate, even with a calculator, is accurate enough (better than 1 minute of time over the 80 year range) for correcting sundials, and has the nice physical explanation as the sum of two terms, one due to obliquity and the other to eccentricity that was used previously in the article. This is not true either for Δt considered as a function of M or for any of its higher order approximations.

References

• "Approximate Solar Coordinates", "Naval Oceanography Portal".
• Burington R S 1949 Handbook of Mathematical Tables and Formulas (Sandusky, Ohio: Handbook Publishers)
• "Computing Greenwich Sidereal Time", "Naval Oceanography Portal".
• Duffett-Smith P 1988 Practical Astronomy with your Calculator Third Edition (Cambridge: Cambridge University Press)
• Heilbron J L 1999 The Sun in the Church, (Cambridge Mass: Harvard University Press|isbn=0-674-85433-0)
• Hinch E J 1991 Perturbation Methods, (Cambridge: Cambridge University Press)
• Hughes D W, et al. 1989, The Equation of Time, Monthly Notices of the Royal Astronomical Society 238 pp 1529-1535
• Meeus, J 1997 Mathematical Astronomy Morsels, (Richmond, Virginia: Willman-Bell)
• Milne R M 1921, Note on the Equation of Time, The Mathematical Gazette 10 (The Mathematical Association) pp 372–375
• Moulton F R 1970 An Introduction to Celestial Mechanics, Second Revised Edition, (New York: Dover)
• Muller M 1995, "Equation of Time - Problem in Astronomy", Acta Phys Pol A 88 Supplement, S-49.
• Roy A E 1978 Orbital Motion, (Adam Hilger|ISBN=0-85274-228-2)
• U.S.Naval Observatory, 2010 Multiyear Interactive Computer Amanac 1800 - 2050, (Richmond, Virginia: Willman-Bell).
• Whitman A M 2007, "A Simple Expression for the Equation of Time", Journal Of the North American Sundial Society 14 pp 29–33