code: def squad(R_in, t_in, t_out):
      if R_in.size == 0 or t_out.size == 0:
          return np.array((), dtype=np.quaternion)
      i_in_for_out = t_in.searchsorted(t_out, side='right')-1
      A = R_in * np.exp((- np.log((~R_in) * np.roll(R_in, -1))
                         + np.log((~np.roll(R_in, 1)) * R_in) * ((np.roll(t_in, -1) - t_in) / (t_in - np.roll(t_in, 1)))
                         ) * 0.25)
      B = np.roll(R_in, -1) * np.exp((np.log((~np.roll(R_in, -1)) * np.roll(R_in, -2))
                                      * ((np.roll(t_in, -1) - t_in) / (np.roll(t_in, -2) - np.roll(t_in, -1)))
                                      - np.log((~R_in) * np.roll(R_in, -1))) * -0.25)
      A[0] = R_in[0]
      A[-1] = R_in[-1]
      B[-2] = R_in[-1]
      B[-1] = R_in[-1] * (~R_in[-2]) * R_in[-1]
      R_ip1 = np.roll(R_in, -1)
      R_ip1[-1] = R_in[-1]*(~R_in[-2])*R_in[-1]
      R_ip1 = np.array(R_ip1[i_in_for_out])
      t_inp1 = np.roll(t_in, -1)
      t_inp1[-1] = t_in[-1] + (t_in[-1] - t_in[-2])
      tau = (t_out - t_in[i_in_for_out]) / ((t_inp1 - t_in)[i_in_for_out])
      R_out = np.squad_vectorized(tau, R_in[i_in_for_out], A[i_in_for_out], B[i_in_for_out], R_ip1)
      return R_out