docs: def transit_sunrise_sunset(dates, lat, lon, delta_t, numthreads):
      """
      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
  
      """
      isnumpy = isinstance(dates, ndarray)
      if isnumpy:
          condition = ((dates % 86400) != 0.0).any()
      else:
          condition = (dates % 86400) != 0.0
      if condition:
          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]))))
      if isnumpy:
          cos_arg[abs(cos_arg) > 1] = np.nan
      else:
          if abs(cos_arg) > 1:
              cos_arg = np.nan
      
      H0 = np.degrees(np.arccos(cos_arg)) % 180
  
      if isnumpy:
          m = np.empty((3, len(utday)))
      else:
          m = np.empty((3, 1))
      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]
      
      if isnumpy:
          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
      else:
          if abs(a) > 2:
              a = a %1
          ap = ttday0_res[2] - ttdayn1_res[2]
          if (abs(ap) > 2):
              ap = ap % 1
          b = ttdayp1_res[1] - ttday0_res[1]
          if (abs(b) > 2):
              b = b % 1
          bp = ttdayp1_res[2] - ttday0_res[2]
          if abs(bp) > 2:
              bp = bp % 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