Skip to content

The Equation of Time

Binder

The "Equation of Time" refers to the observed difference between a given time, say noon UTC, and the time at which the sun is directly overhead. This phenomenon is due to two effects; the tilt of the Earth's axis \(\alpha\), and the elliptical nature of the Earth's orbit around the sun owing to a non-zero eccentricity \(e\).

import numpy as np
import matplotlib.pyplot as plt
from analemma import orbit, geometry

earth = orbit.PlanetParameters.earth()

t = earth.daily_noons()
psi = earth.rotation_angle(t)

pi = np.pi

# hour angle with circular orbit (zero eccentricity, linear in t)
sigma_alpha_only = np.mod(pi + earth.rho + earth.om_y*t, 2*pi)
mu_alpha_only = np.mod(geometry.hour_angle(earth.alpha, sigma_alpha_only, psi), 2*pi)

# hour angle with zero planetary rotation axis tilt
sigma_eot = np.mod(pi + earth.rho + orbit.orbital_angle(orbit.spinor_time(t)), 2*pi)
mu_e_only = np.mod(geometry.hour_angle(0.0, sigma_eot, psi), 2*pi)

# hour angle with non-zero eccentricity and axis tilt
mu_eot = np.mod(geometry.hour_angle(earth.alpha, sigma_eot, psi), 2*pi)

def s2d(seconds):
    "Convert seconds to days"
    return seconds / 24 / 3600

def r2d(radians):
    "Convert radians to degrees"
    return radians / pi * 180

fig, ax = plt.subplots()
# 1 degree is 4 minutes if 360 degrees is 24 hours
ax.plot(s2d(t), 4*r2d(pi - mu_alpha_only), label="Effective of Earth's axis tilt")
ax.plot(s2d(t), 4*r2d(pi - mu_e_only), label="Effect of Keplerian orbit")
ax.plot(s2d(t), 4*r2d(pi - mu_eot), label="Effect of both")
ax.grid()
ax.legend()
ax.set_xlabel("Days since perihelion")
ax.set_ylabel("Minutes")
_ = ax.set_title("The Equation of Time")

png

Exact Parametric Expression

While the Equation of Time is often approximated by two terms in a Fourier series (see for example the presentation on wikipedia), the most accurate treatments account for the effect of the moon and other factors by integrating the equations of motion numerically.

There is a middle ground, however - a parametric expression for the Equation of Time that is exact under the assumptions made here, namely a Keplerian orbit for the planet around its star. This exact expression is the one implemented in analemma, and used to create the above plot.

The minutes shown on the y-axis are related to the difference between the apparent hour angle \(\mu(t)\) and mean hour angle \(\mu_m(t)\) of the sun by a scaled factor of \(180/\pi \times 4\) translating from radians to degrees (\(180/\pi\)), then to minutes (\(\times 4\)). In terms of hour angles, the Equation of Time is given by

\[ \Delta_\mu(t) \equiv \mu(t) - \mu_m(t) = \tan^{-1}\left(\frac{\tan(\sigma_t) - \tan(\psi_t)\cos(\alpha))}{\tan(\sigma_t)\tan(\psi_t) + \cos(\alpha)}\right) - \frac{N+1}{N}\frac{2\pi\,N}{T_y}\,(t-t_p)\]

where

  • \(\alpha\) is the tilt of the planet's axis of rotation relative to the plane of its orbit,
  • \(\sigma_t\) and \(\psi_t\) measure the progress of the planet around its orbit and the planet's rotation, respectively,
  • N is the number of days in a mean year, approximately 365.2422 for Earth
  • \(T_y\) is the length of a year, \(365 \times 24 \times 3600\) seconds for Earth
  • \(t_p\) is the time of the most recent perihelion

The angles \(\alpha\), \(\sigma_t\) and \(\psi_t\) (with \(t\)-subscripts added here to emphasize time-dependence) are defined at the beginning of Setup and Definitions, and the expression for \(\mu(t)\) can be obtained by setting the gnomon inclination \(\iota\) equal to the (\(90^\circ\) minus) latitude \(\theta\) in analemma.algebra.result.hour_angle_tan, derived in The Hour Angle. The factor of \(N\) in the numerator and denominator of the second term cancels, but in this form we can relate it to the derivation more easily (see Section 5 of the Reference Article in the project homepage).

The remainder of the Equation of Time lives in the time-dependence of \(\sigma_t\) and \(\psi_t\). The latter is simplest and is given by

\[\psi_t = \rho + \frac{N+1}{N}\,(t-t_p)\]

where \(\rho\) is the difference in \(\sigma_t\) between perihelion and the winter solstice (see Setup and Definitions). In contrast to \(\psi_t\)'s linear time dependence, \(\sigma_t\) is only available in terms of parameterization based on the planetary orbit, defined by \(a\), \(b\) and \(e\), the semi-major axis, semi-minor axis and eccentricity respectively. Given the parameter \(\tau\), we have

\[\sigma_t(\tau) = \frac{b}{a}\frac{\sin(\omega\tau)}{e + \cos(\omega\tau)}\]

and

\[t(\tau) = a\left(\tau + \frac{e}{\omega}\sin(\omega\tau)\right)\]

where \(\omega = \frac{2\pi a}{T_y}\). This is a standard result in elliptical orbits, whose derivation is given in a non-standard but elegant way in Appendix A of the Reference Article.

To evaluate \(\Delta_\mu(t)\) for a given time, say \(t_\ast\) we first need to invert \(t(\tau_\ast) = t_\ast\) to find \(\tau_\ast\), then evaluate \(\sigma_t(\tau_\ast)\). This is performed above by means of analemma.orbit.spinor_time and analemma.orbit.orbital_angle respectively. We can evaluate everything else that goes into \(\Delta_\mu(t_\ast)\) directly.

The inversion of \(t(\tau)\) calls for a numerical method because \(t(\tau) = t_\ast\) is a transcendental equation, but this is simpler than numerical integration, and becomes simpler still when we realize that \(t(\tau)\) is monotonic in \(\tau\), and so can be inverted by interpolation with minimal runtime cost. This is the approach taken in analemma.orbit.spinor_time.

For a more detailed look at the properties of Earth's orbit including the formulae used above, see Orbit Analysis.