import matplotlib.pyplot as plt
from astropy.visualization import simple_norm
from photutils.datasets import (make_model_image, make_model_params,
                                make_noise_image)
from photutils.psf import GaussianPSF

model = GaussianPSF()
shape = (500, 500)
n_sources = 500
params = make_model_params(shape, n_sources, x_name='x_0',
                           y_name='y_0', min_separation=5,
                           flux=(100, 500), x_fwhm=(1, 3),
                           y_fwhm=(1, 3), theta=(0, 90), seed=123)
model_shape = (25, 25)
data = make_model_image(shape, model, params, model_shape=model_shape,
                        x_name='x_0', y_name='y_0')

noise = make_noise_image(shape, distribution='gaussian', mean=5,
                         stddev=2, seed=123)
data += noise

fig, ax = plt.subplots()
norm = simple_norm(data, 'sqrt', percent=99)
ax.imshow(data, norm=norm, origin='lower')
ax.set_title('Simulated image')