import numpy as np
import matplotlib.pyplot as plt
from astropy.modeling import models
X, Y = np.mgrid[:100,:100]
PSF = models.Moffat2D()
p = (50, 50, 50, 5, 2)
Z = PSF.evaluate(X, Y, *p)
plt.imshow(Z, origin='loxer')
plt.show()