# Calculation of average FWHM by PSF fitting of all the detected stars in the fits image

I need to calculate the FWHM of all the stars detected via the DaoStarFinder package/any other package by performing PSF fitting in Python. Here is the code for my detection.

hdulist = fits.open('F:/TRGB project/2214868p/ccd39.fits')
data = hdulist.data

pixelscale = np.mean(wcs.utils.proj_plane_pixel_scales(image_wcs) * u.degree).to('arcsec')

# Find the sources using DaoFinder

from photutils import DAOStarFinder
from astropy.stats import sigma_clipped_stats

mean, median, std = sigma_clipped_stats(data, sigma=3.0, maxiters=5)

fwhm=1.5 / float(pixelscale/u.arcsec)
sigma_detection = 5.0
daofind = DAOStarFinder(fwhm=fwhm, threshold=sigma_detection*std)
sources = daofind(hdulist.data - median)


Now I need to calculate the FWHM of the stars, take their average, and then use this value as a parameter in SExtractor to distinguish between stars and galaxies. How do I calculate the FWHM of the stars by using PSF fitting method? Here is what I tried by selecting a particular star.

# Median bkg subtracted image
bkg_subtracted_image = data - median

flux_counts = []
pixel_distance = []

x_cen = int(source['x'])
y_cen = int(source['y'])

# Pixels around detection loop
flux_counts.append(((bkg_subtracted_image[y][x]) / source['peak']))
pixel_distance.append(np.linalg.norm((source['x'] - x, source['y'] - y)))

from astropy.modeling import models, fitting

model = 'moffat'

if model == 'gaussian':
# Fit the data using a Gaussian
fwhm_best_guess = 1.5
model_init = models.Gaussian1D(amplitude=1.0, mean=0, stddev=fwhm_best_guess)
elif model == 'moffat':
# Fit the data using a Moffat
model_init = models.Moffat1D(amplitude=1.0, x_0=0, gamma=2., alpha=3.5)
else:
raise Exception("Unknown model type: %s. Must be gaussian or moffat." % model)

fitter = fitting.SimplexLSQFitter()
fitted_model = fitter(model_init, pixel_distance, flux_counts)

# FHWM conversion
if model == 'gaussian':
iFWHM = 2.355 * fitted_model.stddev * pixelscale.value
elif model == 'moffat':
iFWHM = fitted_model.gamma * 2 * np.sqrt(2 ** (1. / fitted_model.alpha) - 1) * pixelscale.value
else:
raise Exception("Unknown model type: %s. Must be gaussian or moffat." % model)

print ("FWHM estimated: ", iFWHM)


I also need to display the PSF-fitted star.