import numpy as np
import matplotlib.pyplot as plt
from astropy.modeling import models
from spectractor.tools import fit_moffat2d_outlier_removal
X, Y = np.mgrid[:100,:100]
PSF = models.Moffat2D()
p = (50, 50, 50, 5, 2)
Z = PSF.evaluate(X, Y, *p)
guess = (45, 48, 52, 4, 2)
bounds = ((1, 10, 10, 1, 1), (100, 90, 90, 10, 10))
fit = fit_moffat2d_outlier_removal(X, Y, Z, guess=guess, bounds=bounds, niter=3)
plt.imshow(Z-fit(X, Y), origin='loxer')
plt.show()