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