def subtract_psfmodel(im, psfmodel, cyx, scaling='peak'):
    """subtract a PSF model from the image

        Parameters
        ----------
        im : array, mandatory
            image
        psfmodel : array, mandatory
            PSF model image
        cyx : tuple, mandatory
            target centroid coordinates `(y, x)`
        scaling : str, optional, default: `peak`
            select scaling method (`peak` fits the height of PSF model to
            the target's image)


        Examples
        --------
        >>> from astropy.io import fits
        >>> from sbpy.imageanalysis import centroid, PSFSubtraction # doctest: +SKIP
        >>> import matplotlib.pyplot as plt # doctest: +SKIP
        >>> hdu = fits.open('test.fits') # doctest: +SKIP
        >>> cyx = centroid(hdu[0].data, (21, 45)) # doctest: +SKIP
        >>> psfmodel = PSFSubtraction.create_psfmodel(hdu[0].data, cyx) # doctest: +SKIP
        >>> diff = PSFSubtraction.subtract_psfmodel(hdu[0].data, psfmodel, cyx) # doctest: +SKIP
        >>> plt.imshow(diff) # doctest: +SKIP

        not yet implemented

        """
