docs: def rotate_points_to_axis(points, axis):
      """Rotate all points of a list, such that `axis==[0,1,0]`
  
      This is accomplished by rotating in the x-z-plane by phi into the
      y-z-plane, then rotation in the y-z-plane by theta up to [0,1,0],
      and finally rotating back in the x-z-plane by -phi.
  
      Parameters
      ----------
      points: list-like with elements of length 3
          The Cartesian points. These should be in the same format as
          produced by `sphere_points_from_angles_and_tilt`.
      axis: list-like, length 3
          The reference axis that will be used to determine the
          rotation angle of the points. The points will be rotated
          about the origin such that `axis` matches [0,1,0].
  
      Returns
      -------
      rotated_points: np.ndarray of shape (N,3)
          The rotated points.
      """
      axis = norm_vec(axis)
      u, v, w = axis
      points = np.array(points)
      # Determine the rotational angle in the x-z plane
      phi = np.arctan2(u, w)
      # Determine the tilt angle w.r.t. the y-axis
      theta = np.arccos(v)
  
      # Negative rotation about y-axis
      Rphi = np.array([
          [np.cos(phi), 0, -np.sin(phi)],
          [0, 1,            0],
          [np.sin(phi), 0,  np.cos(phi)],
      ])
  
      # Negative rotation about x-axis
      Rtheta = np.array([
          [1,             0,              0],
          [0, np.cos(theta),  np.sin(theta)],
          [0, -np.sin(theta),  np.cos(theta)],
      ])
  
      DR1 = np.dot(Rtheta, Rphi)
      # Rotate back by -phi such that effective rotation was only
      # towards [0,1,0].
      DR = np.dot(Rphi.T, DR1)
      rotpoints = np.zeros((len(points), 3))
      for ii, pnt in enumerate(points):
          rotpoints[ii] = np.dot(DR, pnt)
  
      # For visualiztaion:
      # import matplotlib.pylab as plt
      # from mpl_toolkits.mplot3d import Axes3D
      # from matplotlib.patches import FancyArrowPatch
      # from mpl_toolkits.mplot3d import proj3d
      #
      # class Arrow3D(FancyArrowPatch):
      #     def __init__(self, xs, ys, zs, *args, **kwargs):
      #         FancyArrowPatch.__init__(self, (0,0), (0,0), *args, **kwargs)
      #         self._verts3d = xs, ys, zs
      #
      #     def draw(self, renderer):
      #         xs3d, ys3d, zs3d = self._verts3d
      #         xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, renderer.M)
      #         self.set_positions((xs[0],ys[0]),(xs[1],ys[1]))
      #         FancyArrowPatch.draw(self, renderer)
      #
      # fig = plt.figure(figsize=(10,10))
      # ax = fig.add_subplot(111, projection='3d')
      # for vec in rotpoints:
      #     u,v,w = vec
      #     a = Arrow3D([0,u],[0,v],[0,w],
      #                 mutation_scale=20, lw=1, arrowstyle="-|>")
      #     ax.add_artist(a)
      #
      # radius=1
      # ax.set_xlabel('X')
      # ax.set_ylabel('Y')
      # ax.set_zlabel('Z')
      # ax.set_xlim(-radius*1.5, radius*1.5)
      # ax.set_ylim(-radius*1.5, radius*1.5)
      # ax.set_zlim(-radius*1.5, radius*1.5)
      # plt.tight_layout()
      # plt.show()
  
      return rotpoints