docs: def evaluate(self, x, y, flux, x_0, y_0):
          """
          Discrete PRF model evaluation.
  
          Given a certain position and flux the corresponding image of
          the PSF is chosen and scaled to the flux. If x and y are
          outside the boundaries of the image, zero will be returned.
  
          Parameters
          ----------
          x : float
              x coordinate array in pixel coordinates.
          y : float
              y coordinate array in pixel coordinates.
          flux : float
              Model flux.
          x_0 : float
              x position of the center of the PRF.
          y_0 : float
              y position of the center of the PRF.
          """
  
          # Convert x and y to index arrays
          x = (x - x_0 + 0.5 + self.prf_shape[1] // 2).astype('int')
          y = (y - y_0 + 0.5 + self.prf_shape[0] // 2).astype('int')
  
          # Get subpixel indices
          y_sub, x_sub = subpixel_indices((y_0, x_0), self.subsampling)
  
          # Out of boundary masks
          x_bound = np.logical_or(x < 0, x >= self.prf_shape[1])
          y_bound = np.logical_or(y < 0, y >= self.prf_shape[0])
          out_of_bounds = np.logical_or(x_bound, y_bound)
  
          # Set out of boundary indices to zero
          x[x_bound] = 0
          y[y_bound] = 0
          result = flux * self._prf_array[int(y_sub), int(x_sub)][y, x]
  
          # Set out of boundary values to zero
          result[out_of_bounds] = 0
          return result