code: def surface_of_section(orbit, plane_ix, interpolate=False):
      
  
      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
docs: 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...