def transit_sunrise_sunset(dates, lat, lon, delta_t, numthreads):

    if ((dates % 86400) != 0.0).any():
        raise ValueError('Input dates must be at 00:00 UTC')

    utday = (dates // 86400) * 86400
    ttday0 = utday - delta_t
    ttdayn1 = ttday0 - 86400
    ttdayp1 = ttday0 + 86400

    # index 0 is v, 1 is alpha, 2 is delta
    utday_res = solar_position(utday, 0, 0, 0, 0, 0, delta_t,
                               0, numthreads, sst=True)
    v = utday_res[0]

    ttday0_res = solar_position(ttday0, 0, 0, 0, 0, 0, delta_t,
                                0, numthreads, sst=True)
    ttdayn1_res = solar_position(ttdayn1, 0, 0, 0, 0, 0, delta_t,
                                 0, numthreads, sst=True)
    ttdayp1_res = solar_position(ttdayp1, 0, 0, 0, 0, 0, delta_t,
                                 0, numthreads, sst=True)
    m0 = (ttday0_res[1] - lon - v) / 360
    cos_arg = ((np.sin(np.radians(-0.8333)) - np.sin(np.radians(lat))
               * np.sin(np.radians(ttday0_res[2]))) /
               (np.cos(np.radians(lat)) * np.cos(np.radians(ttday0_res[2]))))
    cos_arg[abs(cos_arg) > 1] = np.nan
    H0 = np.degrees(np.arccos(cos_arg)) % 180

    m = np.empty((3, len(utday)))
    m[0] = m0 % 1
    m[1] = (m[0] - H0 / 360)
    m[2] = (m[0] + H0 / 360)

    # need to account for fractions of day that may be the next or previous
    # day in UTC
    add_a_day = m[2] >= 1
    sub_a_day = m[1] < 0
    m[1] = m[1] % 1
    m[2] = m[2] % 1
    vs = v + 360.985647 * m
    n = m + delta_t / 86400

    a = ttday0_res[1] - ttdayn1_res[1]
    a[abs(a) > 2] = a[abs(a) > 2] % 1
    ap = ttday0_res[2] - ttdayn1_res[2]
    ap[abs(ap) > 2] = ap[abs(ap) > 2] % 1
    b = ttdayp1_res[1] - ttday0_res[1]
    b[abs(b) > 2] = b[abs(b) > 2] % 1
    bp = ttdayp1_res[2] - ttday0_res[2]
    bp[abs(bp) > 2] = bp[abs(bp) > 2] % 1
    c = b - a
    cp = bp - ap

    alpha_prime = ttday0_res[1] + (n * (a + b + c * n)) / 2
    delta_prime = ttday0_res[2] + (n * (ap + bp + cp * n)) / 2
    Hp = (vs + lon - alpha_prime) % 360
    Hp[Hp >= 180] = Hp[Hp >= 180] - 360

    h = np.degrees(np.arcsin(np.sin(np.radians(lat)) *
                             np.sin(np.radians(delta_prime)) +
                             np.cos(np.radians(lat)) *
                             np.cos(np.radians(delta_prime))
                             * np.cos(np.radians(Hp))))

    T = (m[0] - Hp[0] / 360) * 86400
    R = (m[1] + (h[1] + 0.8333) / (360 * np.cos(np.radians(delta_prime[1])) *
                                   np.cos(np.radians(lat)) *
                                   np.sin(np.radians(Hp[1])))) * 86400
    S = (m[2] + (h[2] + 0.8333) / (360 * np.cos(np.radians(delta_prime[2])) *
                                   np.cos(np.radians(lat)) *
                                   np.sin(np.radians(Hp[2])))) * 86400

    S[add_a_day] += 86400
    R[sub_a_day] -= 86400

    transit = T + utday
    sunrise = R + utday
    sunset = S + utday

    return transit, sunrise, sunset
Calculate the sun transit, sunrise, and sunset
dataset/ETHPy150Open pvlib/pvlib-python/pvlib/spa.py/transit_sunrise_sunset