# Examples

## Fitting a PSF

Here is a brief example which shows how to construct a loss function for fitting a PSFModel to some data.

using PSFModels
using PSFModels: fit
using HCIDatasets: BetaPictoris
using Plots
using Statistics

# convenience function for plotting
function imshow(data; kwargs...)
xlim = extrema(axes(data, 1))
ylim = extrema(axes(data, 2))
heatmap(transpose(data); xlim=xlim, ylim=ylim,
aspect_ratio=1, clims=(1e-5, Inf), kwargs...)
end

# get a PSF from HCIDatasets.jl;
psf = BetaPictoris[:psf]
inds = CartesianIndices(psf)

imshow(psf)

We can fit this data with a variety of models, here showcasing the flexible PSFModels.fit function.

### Gaussian

Using gaussian

params = (x=20, y=20, fwhm=5, amp=0.1)
P_gauss, mod_gauss = fit(gaussian, params, psf)
pairs(P_gauss)
pairs(::NamedTuple) with 4 entries:
:x    => 20.0191
:y    => 20.01
:fwhm => 4.80293
:amp  => 0.100216
plot(
imshow(psf, title="Data"),
imshow(mod_gauss.(inds), title="Model"),
cbar=false,
ticks=false,
xlabel="",
ylabel="",
layout=2,
size=(600, 300)
)

and now using a rotated, elliptical Gaussian

params = (x=20, y=20, fwhm=(5, 5), amp=0.1, theta=0)
P_ellip, mod_ellip = fit(gaussian, params, psf)
pairs(P_ellip)
pairs(::NamedTuple) with 5 entries:
:x     => 20.0183
:y     => 20.01
:fwhm  => (4.76786, 4.83557)
:amp   => 0.100245
:theta => 0.0
plot(
imshow(psf, title="Data"),
imshow(mod_ellip.(inds), title="Model"),
cbar=false,
ticks=false,
xlabel="",
ylabel="",
layout=2,
size=(600, 300)
)

### Airy disk

Now with airydisk

params = (x=20, y=20, fwhm=5, amp=0.1, ratio=0.3)
P_airy, mod_airy = fit(airydisk, params, psf)
pairs(P_airy)
pairs(::NamedTuple) with 5 entries:
:x     => 19.9891
:y     => 20.0061
:fwhm  => 5.05739
:amp   => 0.102935
:ratio => 0.375187
plot(
imshow(psf, title="Data"),
imshow(mod_airy.(inds), title="Model"),
cbar=false,
ticks=false,
xlabel="",
ylabel="",
layout=2,
size=(600, 300)
)

### Moffat

And finally, with moffat

params = (x=20, y=20, fwhm=(5, 5), amp=0.1, theta=0, alpha=2)
P_moff, mod_moff = fit(moffat, params, psf)
pairs(P_moff)
pairs(::NamedTuple) with 6 entries:
:x     => 20.0359
:y     => 20.0151
:fwhm  => (4.96809, 4.9635)
:amp   => 0.116352
:theta => 0.0
:alpha => 1.46154
plot(
imshow(psf, title="Data"),
imshow(mod_moff.(inds), title="Model"),
cbar=false,
ticks=false,
xlabel="",
ylabel="",
layout=2,
size=(600, 300)
)

### Changing optimization parameters

Any keyword arguments get passed on to Optim.optimize, and you can change the algorithm used with the alg keyword

# load Optim.jl to use the Newton method
using Optim

params = (x=20, y=20, fwhm=(5, 5), amp=0.1, theta=0, alpha=2)
P_moff, mod_moff = fit(moffat, params, psf; alg=Newton())
pairs(P_moff)
pairs(::NamedTuple) with 6 entries:
:x     => 20.0359
:y     => 20.0151
:fwhm  => (4.96808, 4.96349)
:amp   => 0.116352
:theta => 0.0
:alpha => 1.46154

We can also "freeze" parameters by creating a named tuple and passing it to func_kwargs

params = (;x=10, y=20, fwhm=(5, 5), amp=0.1)
func_kwargs = (;alpha=2)
P_moff2, mod_moff2 = fit(moffat, params, psf; func_kwargs)
pairs(P_moff2)
pairs(::NamedTuple) with 4 entries:
:x    => 20.0349
:y    => 20.0144
:fwhm => (6.51242, 6.51775)
:amp  => 0.110515
plot(
imshow(psf, title="Data"),
imshow(mod_moff2.(inds), title="Model"),
cbar=false,
ticks=false,
xlabel="",
ylabel="",
layout=2,
size=(600, 300)
)