def transit_sunrise_sunset(dates, lat, lon, delta_t, numthreads):
if ((dates % <NUM_LIT>) != <NUM_LIT:0.0>).any():<EOL><INDENT>raise ValueError('<STR_LIT>')<EOL><DEDENT>utday = (dates // <NUM_LIT>) * <NUM_LIT><EOL>ttday0 = utday - delta_t<EOL>ttdayn1 = ttday0 - <NUM_LIT><EOL>ttdayp1 = ttday0 + <NUM_LIT><EOL>utday_res = solar_position(utday, <NUM_LIT:0>, <NUM_LIT:0>, <NUM_LIT:0>, <NUM_LIT:0>, <NUM_LIT:0>, delta_t,<EOL><NUM_LIT:0>, numthreads, sst=True)<EOL>v = utday_res[<NUM_LIT:0>]<EOL>ttday0_res = solar_position(ttday0, <NUM_LIT:0>, <NUM_LIT:0>, <NUM_LIT:0>, <NUM_LIT:0>, <NUM_LIT:0>, delta_t,<EOL><NUM_LIT:0>, numthreads, sst=True)<EOL>ttdayn1_res = solar_position(ttdayn1, <NUM_LIT:0>, <NUM_LIT:0>, <NUM_LIT:0>, <NUM_LIT:0>, <NUM_LIT:0>, delta_t,<EOL><NUM_LIT:0>, numthreads, sst=True)<EOL>ttdayp1_res = solar_position(ttdayp1, <NUM_LIT:0>, <NUM_LIT:0>, <NUM_LIT:0>, <NUM_LIT:0>, <NUM_LIT:0>, delta_t,<EOL><NUM_LIT:0>, numthreads, sst=True)<EOL>m0 = (ttday0_res[<NUM_LIT:1>] - lon - v) / <NUM_LIT><EOL>cos_arg = ((np.sin(np.radians(-<NUM_LIT>)) - np.sin(np.radians(lat))<EOL>* np.sin(np.radians(ttday0_res[<NUM_LIT:2>]))) /<EOL>(np.cos(np.radians(lat)) * np.cos(np.radians(ttday0_res[<NUM_LIT:2>]))))<EOL>cos_arg[abs(cos_arg) > <NUM_LIT:1>] = np.nan<EOL>H0 = np.degrees(np.arccos(cos_arg)) % <NUM_LIT><EOL>m = np.empty((<NUM_LIT:3>, len(utday)))<EOL>m[<NUM_LIT:0>] = m0 % <NUM_LIT:1><EOL>m[<NUM_LIT:1>] = (m[<NUM_LIT:0>] - H0 / <NUM_LIT>)<EOL>m[<NUM_LIT:2>] = (m[<NUM_LIT:0>] + H0 / <NUM_LIT>)<EOL>add_a_day = m[<NUM_LIT:2>] >= <NUM_LIT:1><EOL>sub_a_day = m[<NUM_LIT:1>] < <NUM_LIT:0><EOL>m[<NUM_LIT:1>] = m[<NUM_LIT:1>] % <NUM_LIT:1><EOL>m[<NUM_LIT:2>] = m[<NUM_LIT:2>] % <NUM_LIT:1><EOL>vs = v + <NUM_LIT> * m<EOL>n = m + delta_t / <NUM_LIT><EOL>a = ttday0_res[<NUM_LIT:1>] - ttdayn1_res[<NUM_LIT:1>]<EOL>a[abs(a) > <NUM_LIT:2>] = a[abs(a) > <NUM_LIT:2>] % <NUM_LIT:1><EOL>ap = ttday0_res[<NUM_LIT:2>] - ttdayn1_res[<NUM_LIT:2>]<EOL>ap[abs(ap) > <NUM_LIT:2>] = ap[abs(ap) > <NUM_LIT:2>] % <NUM_LIT:1><EOL>b = ttdayp1_res[<NUM_LIT:1>] - ttday0_res[<NUM_LIT:1>]<EOL>b[abs(b) > <NUM_LIT:2>] = b[abs(b) > <NUM_LIT:2>] % <NUM_LIT:1><EOL>bp = ttdayp1_res[<NUM_LIT:2>] - ttday0_res[<NUM_LIT:2>]<EOL>bp[abs(bp) > <NUM_LIT:2>] = bp[abs(bp) > <NUM_LIT:2>] % <NUM_LIT:1><EOL>c = b - a<EOL>cp = bp - ap<EOL>alpha_prime = ttday0_res[<NUM_LIT:1>] + (n * (a + b + c * n)) / <NUM_LIT:2><EOL>delta_prime = ttday0_res[<NUM_LIT:2>] + (n * (ap + bp + cp * n)) / <NUM_LIT:2><EOL>Hp = (vs + lon - alpha_prime) % <NUM_LIT><EOL>Hp[Hp >= <NUM_LIT>] = Hp[Hp >= <NUM_LIT>] - <NUM_LIT><EOL>h = np.degrees(np.arcsin(np.sin(np.radians(lat)) *<EOL>np.sin(np.radians(delta_prime)) +<EOL>np.cos(np.radians(lat)) *<EOL>np.cos(np.radians(delta_prime))<EOL>* np.cos(np.radians(Hp))))<EOL>T = (m[<NUM_LIT:0>] - Hp[<NUM_LIT:0>] / <NUM_LIT>) * <NUM_LIT><EOL>R = (m[<NUM_LIT:1>] + (h[<NUM_LIT:1>] + <NUM_LIT>) / (<NUM_LIT> * np.cos(np.radians(delta_prime[<NUM_LIT:1>])) *<EOL>np.cos(np.radians(lat)) *<EOL>np.sin(np.radians(Hp[<NUM_LIT:1>])))) * <NUM_LIT><EOL>S = (m[<NUM_LIT:2>] + (h[<NUM_LIT:2>] + <NUM_LIT>) / (<NUM_LIT> * np.cos(np.radians(delta_prime[<NUM_LIT:2>])) *<EOL>np.cos(np.radians(lat)) *<EOL>np.sin(np.radians(Hp[<NUM_LIT:2>])))) * <NUM_LIT><EOL>S[add_a_day] += <NUM_LIT><EOL>R[sub_a_day] -= <NUM_LIT><EOL>transit = T + utday<EOL>sunrise = R + utday<EOL>sunset = S + utday<EOL>return transit, sunrise, sunset<EOL>
Calculate the sun transit, sunrise, and sunset
for a set of dates at a given location.

Parameters
----------
dates : array
    Numpy array of ints/floats corresponding to the Unix time
    for the dates of interest, must be midnight UTC (00:00+00:00)
    on the day of interest.
lat : float
    Latitude of location to perform calculation for
lon : float
    Longitude of location
delta_t : float
    Difference between terrestrial time and UT. USNO has tables.
numthreads : int
    Number to threads to use for calculation (if using numba)

Returns
-------
tuple : (transit, sunrise, sunset) localized to UTC