code: def circulation(self):
          
          L = self.angular_momentum()
  
          # if only 2D, add another empty axis
          if L.ndim == 2:
              single_orbit = True
              L = L[...,None]
          else:
              single_orbit = False
  
          ndim,ntimes,norbits = L.shape
  
          # initial angular momentum
          L0 = L[:,0]
  
          # see if at any timestep the sign has changed
          circ = np.ones((ndim,norbits))
          for ii in range(ndim):
              cnd = (np.sign(L0[ii]) != np.sign(L[ii,1:])) | \
                    (np.abs(L[ii,1:]).value < 1E-13)
              ix = np.atleast_1d(np.any(cnd, axis=0))
              circ[ii,ix] = 0
  
          circ = circ.astype(int)
          if single_orbit:
              return circ.reshape((ndim,))
          else:
              return circ
docs: Determine which axes the Orbit circulates around by checking
          whether there is a change of sign of the angular momentum
          about an axis. Returns a 2D array with ``ndim`` integers per orbit
          point. If a box orbit, all integers will be 0. A 1 indicates
          circulation about the corresponding axis.
  
          TODO: clockwise / counterclockwise?
  
          For example, for a single 3D orbit:
  
          - Box and boxlet = [0,0,0]
          - z-axis (short-axis) tube = [0,0,1]
          - x-axis (long-axis) tube = [1,0,0]
  
          Returns
          -------
          circulation : :class:`numpy.ndarray`
              An array that specifies whether there is circulation about any of
              the axes of the input orbit. For a single orbit, will return a
              1D array, but for multiple orbits, the shape will be
              ``(3, norbits)``.