from numpy import sin,cos,tan,arcsin,arccos,arctan, pi, modf
from datetime import datetime

def Length_of_day(latitude,day_of_year):
    """
    latitude = latitude in degrees (eg. 39.3292  for Athens Ohio)

    Example:
    import datetime
    day_of_year = datetime.datetime.now().timetuple().tm_yday

    Sunrise/Sunset is when the top of the sun is apparently even with horizon
    then p = 0.8333

    Paper: A Model Comparison for Daylength as a Function of Latitude and Day
of Year (1995)

    Author: Edward J Rykiel and Randal Stahl
    """
    J = day_of_year
    L = latitude
    p = 0.8333 # daylength coefficient
    theta = 0.2163108 + 2 * arctan(0.9671396 * tan(.00860 * (J - 186)))
    phi = arcsin(0.39795 * cos(theta))

    D = 24 - (24 / pi) * arccos((sin(p * pi / 180) +\
                                 sin(L * pi / 180) *\
                                 sin(phi)) / (cos(latitude * pi / 180) * cos(phi)))

    # get hours and minutes
    frac, whole = modf(D)
    hour = whole
    mint = frac * 60

    s = 'Length of day is: {:.0f} hours {:.0f} minutes'.format(hour,mint)
    print(s)
    return D

latitude = 39.3292 # Athens, OH
day_of_year = datetime.now().timetuple().tm_yday

Length_of_day(latitude, day_of_year);