code: def qrotate(vector, axis, angle):
      
      n_axis = axis / vnorm(axis)
      sin_angle = np.expand_dims(np.sin(angle / 2), 0)
      if np.ndim(n_axis) == 1:
          n_axis = np.expand_dims(n_axis, 1)
          p__ = np.dot(n_axis, sin_angle)[:, np.newaxis]
      else:
          p__ = n_axis * sin_angle
  
      q__ = Quaternion(np.cos(angle / 2), p__)
      shape = vector.shape
      return np.einsum("kj, ikj->ij",
                       vector.reshape((3, -1)),
                       q__.rotation_matrix()[:3, :3]).reshape(shape)
docs: Rotate *vector* around *axis* by *angle* (in radians).
  
      *vector* is a matrix of column vectors, as is *axis*.
      This function uses quaternion rotation.