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