math tool to compute time of sunrise and sunset

Question

What is the mathematical tool that enable to compute time of sunrise and time of sunset at a specific geographical location in the same part of a year, say given month and day exactly?

I have noticed that the time of sunrise and the time of sunset is the same every year in a given month and day, say February 4th for example (at a specific geographical location of course).

Thank you.

Answer 1

The rise and set equations below are from Astronomical Algorithms, and the approximate Sun position is from the Astronomical Almanac. The rise/set algorithm works for any celestial object given it's Right Ascension and Declination. So, for the Sun, you must first compute the RA/Dec of the Sun and use that in the rise/set algorithm. The $h_0$ parameter accounts for the fact that Sun rise/set is defined by the upper limb of the Sun, and a standard atmospheric refraction model. An implementation of the below algorithm, including Julian Date, Sun position, and rise and set is implemented here: Rise and Set Algorithm.

Both algorithms require the Julian Date rather than a Gregorian Calendar Date. Since almost all programming languages can supply the date as the number of seconds or milliseconds elapsed since Jan 1 1970 (Unix Time), the routines below convert between Unix Time and Julian Date. The implementation is in JavaScript but will be trivial to convert to any other language. Since JavaScript returns Unix Time in milliseconds, these functions use that, if you're language uses seconds instead, replace 86400000 below with 86400. Other methods of computing the Julian Date are here: Julian Date Algorithms.

function JulianDateFromUnixTime(t){
    //Not valid for dates before Oct 15, 1582
    return (t / 86400000) + 2440587.5;
}

function UnixTimeFromJulianDate(jd){
    //Not valid for dates before Oct 15, 1582
    return (jd-2440587.5)*86400000;
}

Now use $jd$ to compute the RA/Dec of the Sun:

This algorithm for the Sun position is accurate to 1 degree between 1950 and 2050 according to the Astronomical Almanac. If a higher precision is necessary (and it usually isn't), a method like VSOP87 can be used. $$ \begin{align*} n&=jd-2451545.0 \\ L&=280.460+0.9856474n \\ g&=375.528+.9856003n \\ \lambda &=L+1.915 \sin g+0.020 \sin 2g \\ \beta&=0.0 \\ \epsilon &=23.439-0.0000004n \\ \cos\lambda \tan \alpha &=\cos\epsilon \sin\lambda \\ \sin \delta &=\sin\epsilon \sin \lambda \end{align*} $$

Where $\lambda$ is the ecliptic longitude, $\beta$ is the ecliptic latitude (always 0), $\alpha$ is the Right Ascension, $\delta$ is the Declination. And $ jd $ is the Julian Date.

Now use $\alpha$ and $\delta$ to compute the rise and set times:

The algorithm below is from Meeus' Astronomical Algorithms, that book considers longitudes to be negative to the East, which is opposite of modern definitions, so set $L = -L$ for most applications.

$$ \cos H_0 = \dfrac{\sin h_0 - \sin \phi \sin \delta }{\cos \phi \cos \delta} \ $$

If $\cos H_0$ < -1 or > 1, the point is either always above or below the horizon.

$T=(jd-2451545.0)/36525.0 $

$\Theta_0 = 280.46061837+360.98564736629*(jd-2451545.0)+0.000387933T^2 - T^3/38710000.0$

$$ \begin{cases} transit & \dfrac{\delta + L - \Theta_0 }{360^{\circ}} \\ \\ rise & transit - \dfrac{H_0}{360^{\circ}} \\ \\ set & transit + \dfrac{H_0}{360^{\circ}} \end{cases} $$

$jd$ is the Julian Date for the date in question.

$\delta$ Declination.

$L$ Longitude.

$\phi$ Latitude.

$h_0$ Apparent rise or set angle, -0.8333 for the Sun, +0.125 for the Moon, and -0.5667 for most other objects.

$\Theta_0$ Greenwich sidereal time at 0h for the day in question.

The results $transit$, $rise$, $set$, are in fractions of a day, so multiply by 24 to get the time in hours. If any of them are greater than 1, or less than 0, add or subtract 1 to bring it in the range from 0 to 1. They will be in the UTC time zone, so conversion to the local time zone may be desired.

Final:

Since the Sun moves throughout the day, it may be desired to produce a more accurate result by recomputing everything again, but using the JD for the rise, set or transit time. In the Arctic Circle, it's possible $H_0$ will indicate no rise or set for the first iteration, but the Sun may move to a position where it does rise or set, so it may be desired to retry values throughout the day.