GAU2FIT

Fit a two-component Gausian to a map of an extended source

Description:

This routine finds the parameters of a one- or two-component Gaussian beam by doing a least squares fit to a supplied 2D NDF containing an image of a compact source. The source need not be a point source - the source is assumed to be a sharp-edged circular disc of specified radius (i.e. a simple model of a planet). On each iteration of the fitting process, this source model is convolved with the candidate beam and the residuals with the supplied image are then found. The parameters of the beam are modified on each iteration to minimise the sum of the squared residuals. The beam consists of one or two (see parameter FITTWO) concentric Gaussians added together. The first Gaussian has a fixed amplitude of 1.0. Each candidate beam is normalised to a total data sum of unity before being used. The following parameters are varied during the minimisation process:

Note, if the two components have the same FWHM, they become degenerate (i.e. indistinguishable), which causes problems for the minimisation process. To avoid this, the cost function that is minimised (i.e. the sum of the squared residuals between the smoothed source model and the supplied map) is multiplied by a factor that increases the cost as the FWHM for the two beams becomes more similar. Consequently the second component will always have a larger FWHM than the first beam.

In addition, if parameter FITBACK is TRUE and the FWHM of the second component is similar in size to (or larger than) the data array, it becomes difficult for the minimisation process to distinguish between the second component and a constant background level, again making minimisation difficult. For this reason very wide second components are discouraged by increasing the cost if the second component FWHM is similar to (or larger than) the size of the array.

In addition large amplitude second components are discouraged by increasing the cost for large amplitude second components.

Parameters:

AMP2 = _DOUBLE (Write)
An output parameter holding the amplitude of the second fitted Gaussian relative to the first. For instance, if AMP2 is 0.15 it means that the amplitude of the second component is 15% of the amplitude of the first component.
BACK = _DOUBLE (Write)
An output parameter holding the fitted background level. This will be set to zero if FITBACK is FALSE.
CENTRE = LITERAL (Write)
An output parameter holding the fitted position of the beam centre, in the celestial co-ordinate frame of the NDF (which may or may not be the current co-ordinate system). The string consists of a space separated list of two formatted axis values.
FITBACK = _LOGICAL (Read)
If FALSE, the background level is fixed at zero throughout the minimisation. If TRUE, the background level is included as one of the free parameters in the fit. [TRUE]
FITTWO = _LOGICAL (Read)
If TRUE, the beam consists of two Gaussians. Otherwise only one Gaussian is used. [TRUE]
FWHM1 = _DOUBLE (Write)
An output parameter holding the FWHM of the first fitted Gaussian in arc-sec.
FWHM2 = _DOUBLE (Write)
An output parameter holding the FWHM of the second fitted Gaussian in arc-sec.
IN = NDF (Read)
Input image containing the source to be fitted. Note, if parameter FITBACK is FALSE, the source should be on a zero background. The WCS information must contain a pair of celestial axes, but they need not be the current coordinate system.
INITCENTRE = LITERAL (Read)
An initial guess at the position of the beam centre, in the celestial co-ordinate frame of the NDF (which may or may not be the current co-ordinate system). A space or comma separated list of formatted axis values should be supplied.
LOGFILE = LITERAL (Read)
Name of a log file in which to store a table holding the beam parameters and cost function at each iteration of the minimisation. This table is in TOPCAT " ascii" format. [!]
OUT = NDF (Write)
Output NDF containing the fitted model. This is the source model convolved with the final beam.
RADIUS = _DOUBLE (Read)
The radius of the sharp-edge circular source, in arc-seconds. Must be no smaller than 1.0 arcsec.
RMS = _DOUBLE (Write)
An output parameter holding the RMS residual between the fitted model and the supplied image.