docs: def rotation_matrix(axis, theta):
      """The Euler–Rodrigues formula.
  
      Return the rotation matrix associated with counterclockwise rotation about
      the given axis by theta radians.
  
      Parameters
      ----------
      axis: vector to rotate around
      theta: rotation angle, in rad
      """
      axis = np.asarray(axis)
      axis = axis / np.linalg.norm(axis)
      a = np.cos(theta / 2)
      b, c, d = -axis * np.sin(theta / 2)
      aa, bb, cc, dd = a * a, b * b, c * c, d * d
      bc, ad, ac, ab, bd, cd = b * c, a * d, a * c, a * b, b * d, c * d
      return np.array([
          [aa + bb - cc - dd, 2 * (bc + ad), 2 * (bd - ac)],
          [2 * (bc - ad), aa + cc - bb - dd, 2 * (cd + ab)],
          [2 * (bd + ac), 2 * (cd - ab), aa + dd - bb - cc],
      ])