def surface_of_section(orbit, plane_ix, interpolate=False):
    """
    Generate and return a surface of section from the given orbit.

    .. warning::

        This is an experimental function and the API may change.

    Parameters
    ----------
    orbit : `~gala.dynamics.Orbit`
    plane_ix : int
        Integer that represents the coordinate to record crossings in. For
        example, for a 2D Hamiltonian where you want to make a SoS in
        :math:`y-p_y`, you would specify ``plane_ix=0`` (crossing the
        :math:`x` axis), and this will only record crossings for which
        :math:`p_x>0`.
    interpolate : bool (optional)
        Whether or not to interpolate on to the plane of interest. This
        makes it much slower, but will work for orbits with a coarser
        sampling.

    Returns
    -------

    Examples
    --------
    If your orbit of interest is a tube orbit, it probably conserves (at
    least approximately) some equivalent to angular momentum in the direction
    of the circulation axis. Therefore, a surface of section in R-z should
    be instructive for classifying these orbits. TODO...show how to convert
    an orbit to Cylindrical..etc...

    """

    w = orbit.w()
    if w.ndim == 2:
        w = w[...,None]

    ndim,ntimes,norbits = w.shape
    H_dim = ndim // 2
    p_ix = plane_ix + H_dim

    if interpolate:
        raise NotImplementedError("Not yet implemented, sorry!")

    # record position on specified plane when orbit crosses
    all_sos = np.zeros((ndim,norbits), dtype=object)
    for n in range(norbits):
        cross_ix = argrelmin(w[plane_ix,:,n]**2)[0]
        cross_ix = cross_ix[w[p_ix,cross_ix,n] > 0.]
        sos = w[:,cross_ix,n]

        for j in range(ndim):
            all_sos[j,n] = sos[j,:]

    return all_sos
Generate and return a surface of section from the given orbit.

    .. warning::

        This is an experimental function and the API may change.

    Parameters
    ----------
    orbit : `~gala.dynamics.Orbit`
    plane_ix : int
        Integer that represents the coordinate to record crossings in. For
        example, for a 2D Hamiltonian where you want to make a SoS in
        :math:`y-p_y`, you would specify ``plane_ix=0`` (crossing the
        :math:`x` axis), and this will only record crossings for which
        :math:`p_x>0`.
    interpolate : bool (optional)
        Whether or not to interpolate on to the plane of interest. This
        makes it much slower, but will work for orbits with a coarser
        sampling.

    Returns
    -------

    Examples
    --------
    If your orbit of interest is a tube orbit, it probably conserves (at
    least approximately) some equivalent to angular momentum in the direction
    of the circulation axis. Therefore, a surface of section in R-z should
    be instructive for classifying these orbits. TODO...show how to convert
    an orbit to Cylindrical..etc...