API/Reference

gaussian([T=Float64], point; x, y, fwhm, amp=1, theta=0, bkg=0)
gaussian([T=Float64], px, py; x, y, fwhm, amp=1, theta=0, bkg=0)

An unnormalized bivariate Gaussian distribution. The position can be specified in (x, y) coordinates as a Tuple, AbstractVector, or as separate arguments. If theta is given, the PSF will be rotated by theta degrees counter-clockwise from the x-axis. If bkg is given it will be added as a scalar to the PSF.

The fwhm can be a scalar (isotropic) or a vector/tuple (diagonal). Keep in mind that theta has no effect for isotropic distributions and is degenerate with the fwhm parameters (i.e., theta=90 is the same as reversing the fwhm tuple)

Functional form

f(x | x̂, FWHM) = exp[-4ln(2) * ||x - x̂|| / FWHM^2]

where x̂ and x are position vectors (indices) ||⋅|| represents the square-distance, and FWHM is the full width at half-maximum. If FWHM is a scalar, the Gaussian distribution will be isotropic. If FWHM is a vector or tuple, the weighting is applied along each axis (diagonal).

normal

An alias for gaussian

airydisk([T=Float64], point; x, y, fwhm, ratio=0, amp=1, theta=0, bkg=0)
airydisk([T=Float64], px, py; x, y, fwhm, ratio=0, amp=1, theta=0, bkg=0)

An unnormalized Airy disk. The position can be specified in (x, y) coordinates as a Tuple, AbstractVector, or as separate arguments. If theta is given, the PSF will be rotated by theta degrees counter-clockwise from the x-axis. If bkg is given it will be added as a scalar to the PSF.

The fwhm can be a scalar (isotropic) or a vector/tuple (diagonal). Keep in mind that theta has no effect for isotropic distributions and is degenerate with the fwhm parameters (i.e., theta=90 is the same as reversing the fwhm tuple)

If ratio is supplied, this will be the Airy pattern for a centrally-obscured aperture (e.g., a Newtonian telescope). This has a slightly expanded functional form, and in general the central Airy disk will be smaller and the first Airy ring will be brighter.

Functional form

The Airy disk is a distribution over the radius r (the square-Euclidean distance)

f(x | x̂, FWHM) = [ 2J₁(q) / q ]^2

where J₁ is the first-order Bessel function of the first kind and

q ≈ π * r * D/ λ ≈ π * r / (0.973 * FWHM)

If user a non-zero central obscuration via ratio, the functional form becomes

f(x | x̂, FWHM, ϵ) = [ 2J₁(q) / q - 2ϵJ₁(ϵq) / q ]^2 / (1 - ϵ^2)^2

where ϵ is the ratio (0 ≤ ϵ < 1).

moffat([T=Float64], point; x, y, fwhm, alpha=1, amp=1, theta=0, bkg=0)
moffat([T=Float64], px, py; x, y, fwhm, alpha=1, amp=1, theta=0, bkg=0)

Two dimensional Moffat model. The position can be specified in (x, y) coordinates as a Tuple, AbstractVector, or as separate arguments. If theta is given, the PSF will be rotated by theta degrees counter-clockwise from the x-axis. If bkg is given it will be added as a scalar to the PSF.

The fwhm can be a scalar (isotropic) or a vector/tuple (diagonal). Keep in mind that theta has no effect for isotropic distributions and is degenerate with the fwhm parameters (i.e., theta=90 is the same as reversing the fwhm tuple)

Functional form

f(x | x̂, FWHM, α) = A / (1 + ||x - x̂|| / (FWHM / 2)^2)^α

where x̂ and x are position vectors (indices) ||⋅|| represents the square-distance, and FWHM is the full width at half-maximum. If FWHM is a vector or tuple, the weighting is applied along each axis.

PSFModels.fit(model, params, image, inds=axes(image); func_kwargs=(;), loss=abs2, maxfwhm=Inf, alg=LBFGS(), kwargs...)

Fit a PSF model (model) defined by the given params as a named tuple of the parameters to fit and their default values. This model is fit to the data in image at the specified inds (by default, the entire array). To pass extra keyword arguments to the model (i.e., to "freeze" a parameter), pass them in a named tuple to func_kwargs. The default loss function is the chi-squared loss, which uses the the square of the difference (i.e., the L2 norm). You can change this to the L1 norm, for example, by passing loss=abs. The maximum FWHM can be set with maxfwhm as a number or tuple.

Additional keyword arguments, as well as the fitting algorithm alg, are passed to Optim.optimize. By default we use forward-mode auto-differentiation (AD) to derive Jacobians for the LBFGS optimization algorithm. Refer to the Optim.jl documentation for more information.

Choosing parameters

The fit function is very powerful because it gives you a great amount of flexibility in the way you fit your models. To demonstrate this, let's start with a simple isotropic gaussian.

model = gaussian
# match params to arguments of PSF
params = (x=20, y=20, fwhm=3, amp=1)

Note that params can follow any order

params = (amp=1, x=20, y=20, fwhm=3)

Now, to extend this interface to the bivariate PSF case, where fwhm is two values, all you need to do is use a tuple or vector

params = (x=20, y=20, fwhm=(3, 3))

and, again, the order does not matter

model = moffat
params = (alpha=1, x=20, y=20, fwhm=3, amp=10)

Fitting a PSF

After selecting your model and parameters, fitting data is easy

P = (x=12, y=13, fwhm=13.2, amp=0.1)
psf = gaussian.(CartesianIndicies(1:25, 1:15); P...)

params, synthpsf = PSFModels.fit(gaussian, P, psf)

here params is a named tuple of the best fitting parameters. It will not include any fixed parameters.

synthpsf is the best-fitting model, for direct comparison with the input data.

psf_fit = synthpsf.(CartesianIndicies(psf))

To alter parameters without fitting them (i.e., "freeze" them) use func_kwargs

P = (x=12, y=13, fwhm=(12.4, 13.2), amp=0.1)
func_kwargs = (alpha=2)
params, synthpsf = PSFModels.fit(moffat, P, psf; func_kwargs)