推估太陽黃經歷時的數值方法
參考太陽直接輻射的天文模型，太陽黃經和時間的關係有底下一階常微分方程式所控制：
$$ \begin{equation} \frac{d t}{d \lambda} = \frac{K'_{rev}}{(1-\mathbf{e} \cos(\lambda - \lambda_1))^2} \end{equation} $$其中\(K'_{rev}\)為一地球公轉之常數、\(\mathbf{e}\)為地球公轉軌道之離心率（約0.0167）、\(\lambda_1\)為遠日點發生時（約7月初）的太陽黃經。
這裡我們將使用頻譜方法（spectral method）解這道微分方程式。我們把函數\(\frac{1}{(1-e \cos(\lambda - \lambda_1))^2}\)做離散取樣後，做離散傅立葉變換（Discrete Fourier Transform）。利用周波數為零的項的係數\(c_0\)，我們可以找出\(K'_{rev}\)：
$$ \begin{equation} K'_{rev} = \frac{T_{rev}}{2 \pi c_0 \Delta \lambda} \end{equation} $$其中\(T_{rev}\)為太陽年（取兩春分點之歷時，約365.2422天）、\(\Delta\lambda\)為太陽黃經的取樣間隔。如此便可以求出太陽黃經的歷時函數\(t(\lambda)\)。
補充說明：\(K'_{rev}\)求法
連續傅立葉變換當中，\(K'_{rev}=\frac{T_{rev}}{2 \pi F(0)}\) ，其中\(F(0)\)為\(\frac{d t}{d \lambda}\)傅立葉變換後周波數為零之對應值。由於\(F(0) = \int^{2 \pi}_{0} t'(\lambda) d \lambda\)可以近似成\(\sum^{N}_{n=0} t'(\lambda_n) \Delta \lambda = c_0 \Delta \lambda\)，我們便有以上\(K'_{rev}\)的近似式。
補充說明：\(t(\lambda)\)求法
一般而言會定義某一函數\(x(\lambda)\)的離散傅立葉變換為\(X(k)=\sum^{N-1}_{n=0} x(\lambda_n) e^{-j \frac{2 \pi nk}{N}}\)，而逆離散傅立葉變換則定義為\(x(\lambda)=\sum^{N-1}_{n=0} X(k) e^{-j \frac{2 \pi nk}{N}} \Delta \lambda\)（此為傅立葉級數的近似）。注意在我們的例子裡，\(\frac{2 \pi n}{N} = \lambda\)。
當我們以逆離散傅立葉變換\(K_{rev} \left( X(0) + \sum^{N-1}_{k=0} X(k) e^{j k\lambda} \right) \Delta \lambda\)近似\(t'(\lambda)\)時，將此級數各項對\(\lambda\)做積分便可得
$$ \begin{equation} \begin{array}{rcl} \Delta t &=& K_{rev} \Delta \lambda \left( X(0) (\lambda - \lambda_0) + \sum^{N-1}_{k=0}\frac{X(k)}{jk} \left( e^{j k\lambda} - e^{j k\lambda_0} \right) \right) \\ &=& K_{rev} \Delta \lambda \left( X(0) \lambda + \sum^{N-1}_{k=0}\frac{X(k)}{jk} e^{j k\lambda} \right) - t(\lambda_0) \end{array} \end{equation} $$所以\(t(\lambda)=K_{rev} \Delta \lambda \left( X(0) \lambda + \sum^{N-1}_{k=0}\frac{X(k)}{jk} e^{j k\lambda} \right)\)。
推估太陽時和太陽黃經關係的數值方法
由太陽赤經與太陽黃經的關係
$$ \begin{equation} \alpha(t) = \tan^{-1}(\cos(\delta_{max}) \tan(\lambda)) \end{equation} $$以及太陽時的一階常微分方程式
$$ \begin{equation} \dot{h}(t)= \frac{d h}{d \lambda} \frac{d \lambda}{dt} = \Omega_{rot} - \dot{\alpha} = \Omega_{rot} - \frac{d \alpha}{d \lambda} \frac{d \lambda}{dt} \end{equation} $$可以推出太陽時對太陽黃經的微分關係：
$$ \begin{equation} \frac{d h}{d \lambda} = \Omega_{rot} \frac{dt}{d \lambda} - \frac{d \alpha}{d \lambda} = \frac{K'_{rev}\Omega_{rot}}{(1-\mathbf{e} \cos(\lambda - \lambda_1))^2} - \frac{\cos(\delta_{max}) (1 + \tan^2(\lambda))}{1 + \cos^2(\delta_{max}) \tan^2(\lambda)} \end{equation} $$同前一節，利用頻譜方法我們可以求出太陽時對太陽黃經的函數\(h(\lambda)\)。
均時差方程式
有了\(h(\lambda)\)後，我們再定義平均太陽時\(\bar{h}(\lambda)\)：這是假設太陽時的時依變化率為定值，並將一太陽年等分成365.2422個太陽日後的均化函數，如此便可推出均時差\(h(\lambda) - \bar{h}(\lambda)\)，此即均時差方程式（equation of time）。
補充說明：\(\bar{h}(\lambda)\)求法
\(\bar{h}(\lambda)\)可直接從\(t(\lambda)\)求出：
$$ \begin{equation} \bar{h}(\lambda) = \frac{2 \pi t(\lambda)}{\hat{T}_{msd}} \end{equation} $$其中\(\hat{T}_{msd}\)是平均太陽日的周期（約86400秒）。
均時差方程式。負值為太陽時比平均太陽時慢的時段、正值為太陽時比平均太陽時快的時段。
程式碼
# Written in R
### Parameters of Earth Orbit
ecc <- .0167 ## eccentricity
lambda_1 <- (185 * 86400 + 22 * 3600 + 27 * 60) / (365.2422 * 86400) * 2 * pi ## eclipse longitude when aphelion occurs in 2021
lambda_2 <- (78 * 86400 + 9 * 3600 + 37 * 60) / (365.2422 * 86400) * 2 * pi ## spring equinox in 2021
T_rev <- 365.2422 * 86400 ## solar year
Omega_rot <- 2 * pi / 86400 * 366.2422 / 365.2422 ## Earth rotation angular speed
delta_max <- 23.4 / 180 * pi ## Eclipse Plane Declination
### Sampling in the lambda and wave number domain
### Lambda is the eclipse longitude
res <- 365 * 24
lambda <- seq(0, 2 * pi, length.out = res + 1)
lambda <- lambda[-length(lambda)]
k_lambda <- 0:(res / 2 - 1)
k_lambda <- c(k_lambda, -rev(1:(res / 2)))
### Time - Lambda Relationship
diff_t <- 1 / (1 - ecc * cos(lambda - lambda_1))^2
DIFF_t <- fft(diff_t)
K_rev <- T_rev / 2 / pi / (Re(DIFF_t[1]) / res)
Re(DIFF_t[1]) / res
INT_t <- DIFF_t / complex(imaginary = k_lambda)
INT_t[1] <- 0
t <- K_rev * fft(INT_t, inverse = T) / res
t <- Re(t) + K_rev * Re(DIFF_t[1]) / res * lambda
### Alpha - Lambda Relationship
### Alpha is the equatorial longitude
diff_alpha <- cos(delta_max) * (1 + tan(lambda - lambda_2)^2) / (1 + cos(delta_max)^2 * tan(lambda - lambda_2)^2)
DIFF_alpha <- fft(diff_alpha)
INT_alpha <- DIFF_alpha / complex(imaginary = k_lambda)
INT_alpha[1] <- 0
alpha <- fft(INT_alpha, inverse = T) / res
alpha <- Re(alpha) + Re(DIFF_alpha[1]) / res * lambda
### h - Lambda Relationship
### h is the local solar time
diff_h <- Omega_rot * K_rev * diff_t - diff_alpha
DIFF_h <- fft(diff_h)
INT_h <- DIFF_h / complex(imaginary = k_lambda)
INT_h[1] <- 0
h <- fft(INT_h, inverse = T) / res
h <- Re(h) + Re(DIFF_h[1]) / res * lambda
## Mean Solar Time and Error
h_bar <- t / 86400 * 2 * pi
h_error <- h - h_bar
plot(t / T_rev * 365.2422, h_error / 2 / pi * 1440, type = "l", xlab = "Day of Year", ylab = "Equation of Time (Minutes)")
lines(c(0, max(t / T_rev * 365.2422)), c(0, 0), lty = 2)