docs: def field_subset(f, inds, rank=0):
      """Return the value of a field at a subset of points.
  
      Parameters
      ----------
      f: array, shape (a1, a2, ..., ad, r1, r2, ..., rrank)
          Rank-r field in d dimensions
      inds: integer array, shape (n, d)
          Index vectors
      rank: integer
          The rank of the field (0: scalar field, 1: vector field and so on).
  
      Returns
      -------
      f_sub: array, shape (n, rank)
          The subset of field values.
      """
      f_dim_space = f.ndim - rank
      if inds.ndim > 2:
          raise Exception('Too many dimensions in indices array')
      if inds.ndim == 1:
          if f_dim_space == 1:
              return f[inds]
          else:
              raise Exception('Indices array is 1d but field is not')
      if inds.shape[1] != f_dim_space:
          raise Exception('Indices and field dimensions do not match')
      # It's magic, don't touch it!
      return f[tuple([inds[:, i] for i in range(inds.shape[1])])]