How to calculate the length of a day ? 🕒

During a recent development on energy time series analysis, we needed to determine the length of a day in hours.

We have therefore started using Python libraries such as pvlib and its module solar_position or suntime to obtain sunrise and sunset times.

However, to avoid dependencies, we continued to dig and came across an original method in a post on Stackoverflow , which proposes a simple method to estimate the length of the day based on latitude and day of the year. This method was developed by [Forsythe et al.]¹ in 1995!

This duration is described by the following function: $$ D = 24 - \frac{24}{\pi} \times \arccos \left[ \frac{\sin(0.8333 \times \frac{\pi}{180}) + \sin(L \times \frac{\pi}{180})\times\sin(P)}{\cos(L \times \frac{\pi}{180} \times \cos(P))} \right] $$

with $D$ the length of a day, in hours, $L$ the latitude in degrees, $J$ the day of the year and $P$ :

$$ P = \arcsin \left[ 0.39795 \times \cos\left( 0.2163108 + 2 \times \arctan\left( 0.9671396 \times \tan\left(0.00860\times(J-186)\right)\right)\right)\right] $$

The value 0.8333 can be adjusted to account for different definitions of the length of a day (see article).

The authors provide comparisons and show that their approach has an error of less than 1 minute for latitudes below 40° North/South, and a maximum error of 7 minutes for a latitude of 60°. maximal de 7 minutes pour une latitude de 60°.

Python implementation

def day_length(J, L):
    """
    -----------------------------------------------------------------------------------------
    Based upon : "A model comparison for daylength as a function of latitude and day of year"
    Forsythe et al., 1995, Ecological Modelling 80 (1995) 87-95
    -----------------------------------------------------------------------------------------
    Parameters
    ----------
    J: int / list of int / array 
        day of the year.
    L: float 
        latitude (in °)

    Returns
    -------
    Lenght of the day(s) in hours
    
    To account for various definitions of daylength, modify the "p" value accordingly.
    * Sunrise/Sunset is when the center of the sun is even with the horizon 
    p = 0
    * Sunrise/Sunset is when the top of the sun is even with horizon
    p = 0.26667
    * Sunrise/Sunset is when the top of the sun is apparently even with horizon
    p = 0.8333

    """
    p = 0.8333
    phi = np.arcsin(
            0.39795 * ( np.cos( 0.2163108 + 2 * np.arctan( 0.9671396 * np.tan( 0.00860 * (J-186) ) ) ) )
        )
    D = 24 - (24/np.pi)*np.arccos(
              ( np.sin( p*np.pi/180 ) + np.sin( L*np.pi/180 ) * np.sin( phi ) ) / (np.cos(L*np.pi/180) * np.cos( phi ) )
        )

    return D