How to use the photutils.psf.models.EPSFModel function in photutils

shape = np.array((y_shape, x_shape))

        # verify odd sizes of shape
        shape = [(i + 1) if i % 2 == 0 else i for i in shape]

        data = np.zeros(shape, dtype=np.float)

        # ePSF origin should be in the undersampled pixel units, not the
        # oversampled grid units. The middle, fractional (as we wish for
        # the center of the pixel, so the center should be at (v.5, w.5)
        # detector pixels) value is simply the average of the two values
        # at the extremes.
        xcenter = stars._max_shape[0] / 2.
        ycenter = stars._max_shape[1] / 2.

        epsf = EPSFModel(data=data, origin=(xcenter, ycenter),
                         oversampling=oversampling, norm_radius=norm_radius,
                         shift_val=shift_val)

        return epsf

epsf = EPSFModel(data=new_epsf, origin=epsf.origin,
                             oversampling=epsf.oversampling,
                             norm_radius=epsf._norm_radius,
                             shift_val=epsf._shift_val, normalize=False)

            # smooth and recenter the ePSF
            epsf._data = self._smooth_epsf(epsf._data)
            epsf._data = self._recenter_epsf(
                epsf, centroid_func=self.recentering_func)

        # return the new ePSF object, but with undersampled grid pixel
        # coordinates
        xcenter = (epsf._data.shape[1] - 1) / 2. / epsf.oversampling[0]
        ycenter = (epsf._data.shape[0] - 1) / 2. / epsf.oversampling[1]

        return EPSFModel(data=epsf._data, origin=(xcenter, ycenter),
                         oversampling=epsf.oversampling,
                         norm_radius=epsf._norm_radius,
                         shift_val=epsf._shift_val)

should be as close as possible to actual centers.  For stars
            than contain weights, a weighted fit of the ePSF to the star
            will be performed.

        Returns
        -------
        fitted_stars : `EPSFStars` object
            The fitted stars.  The ePSF-fitted center position and flux
            are stored in the ``center`` (and ``cutout_center``) and
            ``flux`` attributes.
        """

        if len(stars) == 0:
            return stars

        if not isinstance(epsf, EPSFModel):
            raise TypeError('The input epsf must be an EPSFModel.')

        # make a copy of the input ePSF
        epsf = epsf.copy()

        # perform the fit
        fitted_stars = []
        for star in stars:
            if isinstance(star, EPSFStar):
                fitted_star = self._fit_star(epsf, star, self.fitter,
                                             self.fitter_kwargs,
                                             self.fitter_has_fit_info,
                                             self.fit_boxsize)

            elif isinstance(star, LinkedEPSFStar):
                fitted_star = []

else:
                    residuals = np.nanmean(residuals, axis=0)

            # interpolate any missing data (np.nan)
            mask = ~np.isfinite(residuals)
            if np.any(mask):
                residuals = _interpolate_missing_data(residuals, mask,
                                                      method='cubic')

                # fill any remaining nans (outer points) with zeros
                residuals[~np.isfinite(residuals)] = 0.

            # add the residuals to the previous ePSF image
            new_epsf = epsf._data + residuals

            epsf = EPSFModel(data=new_epsf, origin=epsf.origin,
                             oversampling=epsf.oversampling,
                             norm_radius=epsf._norm_radius,
                             shift_val=epsf._shift_val, normalize=False)

            # smooth and recenter the ePSF
            epsf._data = self._smooth_epsf(epsf._data)
            epsf._data = self._recenter_epsf(
                epsf, centroid_func=self.recentering_func)

        # return the new ePSF object, but with undersampled grid pixel
        # coordinates
        xcenter = (epsf._data.shape[1] - 1) / 2. / epsf.oversampling[0]
        ycenter = (epsf._data.shape[0] - 1) / 2. / epsf.oversampling[1]

        return EPSFModel(data=epsf._data, origin=(xcenter, ycenter),
                         oversampling=epsf.oversampling,