import numpy as np
import matplotlib.pyplot as plt
from astropy.modeling import models
from spectractor.tools import fit_gauss2d_outlier_removal
X, Y = np.mgrid[:50,:50]
PSF = models.Gaussian2D()
p = (50, 25, 25, 5, 5, 0)
Z = PSF.evaluate(X, Y, *p)
guess = (45, 20, 20, 7, 7, 0)
bounds = ((1, 10, 10, 1, 1, -90), (100, 40, 40, 10, 10, 90))
fit = fit_gauss2d_outlier_removal(X, Y, Z, guess=guess, bounds=bounds, circular=True)
plt.imshow(Z-fit(X, Y), origin='lower')
plt.show()