docs: def _find_bounding_points(self, x, y):
          """
          Find the indices of the grid points that bound the input
          ``(x, y)`` position.
  
          Parameters
          ----------
          x, y : float
              The ``(x, y)`` position where the PSF is to be evaluated.
  
          Returns
          -------
          indices : list of int
              A list of indices of the bounding grid points.
          """
  
          if not np.isscalar(x) or not np.isscalar(y):  # pragma: no cover
              raise TypeError('x and y must be scalars')
  
          if (x < self._xgrid_min or x > self._xgrid_max or
                  y < self._ygrid_min or y > self._ygrid_max):  # pragma: no cover
              raise ValueError('(x, y) position is outside of the region '
                               'defined by grid of PSF positions')
  
          x0 = self._find_bounds_1d(self._xgrid, x)
          y0 = self._find_bounds_1d(self._ygrid, y)
          points = list(itertools.product(self._xgrid[x0:x0 + 2],
                                          self._ygrid[y0:y0 + 2]))
  
          indices = []
          for xx, yy in points:
              indices.append(np.argsort(np.hypot(self._grid_xpos - xx,
                                                 self._grid_ypos - yy))[0])
  
          return indices