The goal is to calculate sunrise and sunset times given a date and the observation site’s latitude $\phi_{o}$ and longitude $\lambda_o$. The example is to compute the sunrise and sunset time on March 23, 1996 at the site 40N, 0E.

Step 1: Compute reference JD at sunrise and sunset:

$JD_{sunrise}=JD_{0h}+6/24 - \lambda_o/360 = 2450165.75$

$JD_{sunset}=JD_{0h}+18/24 - \lambda_o/360 = 2450166.25$

From Step 2 to Step 7, the numbers are assumed to use $JD_{sunrise}$ as the input. For sunset, we need to compute them using $JD_{sunset}$.

Step 2: Compute Julian centuries:

$T_{UT1}=\frac{JD_{UT1}-2451545}{36525}=-0.037762$ (w.r.t. J2000)

Step 3: Use Julian centuries to compute Sun’s mean ecliptic longitude:

$\bar{\lambda}{ecliptic}=280.4606184^{\circ}+36000.77005361 T{UT1}=1.006488^\circ$

Step 4: Assuming $T_{UT1} \sim T_{TDB}$ (Universal Time 1 vs. Barycentric Dynamical Time), compute Sun’s mean anomaly:

$M=357.5291092^\circ + 35999.05034 T_{TDB}=78.139^\circ$

Step 5: Compute Sun’s true ecliptic longitude. The difference to mean ecliptic longitude is to consider the orbital eccentricity.

$\lambda_{ecliptic}=\bar{\lambda}_{ecliptic} + 1.91466471^\circ SIN(M) + 0.019994643 SIN (2M) = 2.88832^\circ$

Step 6: Compute the obliquity of the Ecliptic (the angle between the Earth’s equatorial plane and the plane of the ecliptic):