docs: def subtract_psf(data, psf, posflux, subshape=None):
      """
      Subtract PSF/PRFs from an image.
  
      Parameters
      ----------
      data : `~astropy.nddata.NDData` or array (must be 2D)
          Image data.
      psf : `astropy.modeling.Fittable2DModel` instance
          PSF/PRF model to be substracted from the data.
      posflux : Array-like of shape (3, N) or `~astropy.table.Table`
          Positions and fluxes for the objects to subtract.  If an array,
          it is interpreted as ``(x, y, flux)``  If a table, the columns
          'x_fit', 'y_fit', and 'flux_fit' must be present.
      subshape : length-2 or None
          The shape of the region around the center of the location to
          subtract the PSF from.  If None, subtract from the whole image.
  
      Returns
      -------
      subdata : same shape and type as ``data``
          The image with the PSF subtracted
      """
  
      if data.ndim != 2:
          raise ValueError('{0}-d array not supported. Only 2-d arrays can be '
                           'passed to subtract_psf.'.format(data.ndim))
  
      #  translate array input into table
      if hasattr(posflux, 'colnames'):
          if 'x_fit' not in posflux.colnames:
              raise ValueError('Input table does not have x_fit')
          if 'y_fit' not in posflux.colnames:
              raise ValueError('Input table does not have y_fit')
          if 'flux_fit' not in posflux.colnames:
              raise ValueError('Input table does not have flux_fit')
      else:
          posflux = Table(names=['x_fit', 'y_fit', 'flux_fit'], data=posflux)
  
      # Set up contstants across the loop
      psf = psf.copy()
      xname, yname, fluxname = _extract_psf_fitting_names(psf)
      indices = np.indices(data.shape)
      subbeddata = data.copy()
  
      if subshape is None:
          indicies_reversed = indices[::-1]
  
          for row in posflux:
              getattr(psf, xname).value = row['x_fit']
              getattr(psf, yname).value = row['y_fit']
              getattr(psf, fluxname).value = row['flux_fit']
  
              subbeddata -= psf(*indicies_reversed)
      else:
          for row in posflux:
              x_0, y_0 = row['x_fit'], row['y_fit']
  
              y = extract_array(indices[0], subshape, (y_0, x_0))
              x = extract_array(indices[1], subshape, (y_0, x_0))
  
              getattr(psf, xname).value = x_0
              getattr(psf, yname).value = y_0
              getattr(psf, fluxname).value = row['flux_fit']
  
              subbeddata = add_array(subbeddata, -psf(x, y), (y_0, x_0))
  
      return subbeddata