<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 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 """


    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