c gaufit fits gaussians to a profile and can write the output to
a miriad dataset, a logfile or the terminal.
The fitting is done using an adapted version of fitting routines in numerical recipes.
Obligatory parameters are: a) either or both of 'in=' and 'parinp=' b) 'rmsest=' c) and either 'estim=' or 'options=findestim'
The ease of fitting depends strongly on the initial estimates. These can be given using the estim= keyword, or they can be automatically found (options=findestim). The latter is usually preferable, though the former may be necessary in pathological cases.
If multiple gaussians are fit but the fit is bad, another try is made with one less gaussian. This is repeated until the fit works. A fit is considered bad if the rms of the residual is higher than 1.8 times the rms estimate, or if the parameters lie outside the range given by cutoff=, crange= and wrange=, or if the fitting routine does not converge. If after all possible ways of retrying a fit the rms is between 1.8 and 5.4 times the rms estimate, accept the fit after all; maybe a low-level component increased the rms, or the profile is not perfectly gaussian, but the fit is somewhat reasonable.
Automatic initial estimates are made by first finding the velocity and amplitude of the peak. An estimated width is found by looking to both sides of the peak for the nearest zero of the 2nd derivative and at where the half-maximum lies. The most consistent combination gives the width estimate. For low S/N profiles the width is found from the integral out to the nearest zero. This estimated component is then subtracted from the profile and the process is repeated until the maximum amplitude is too low, or until the maximum number of gaussians has been found.
If parinp= is used gaussian parameters are taken from this dataset for all pixels outside the specified region. Inside the region new fits are made. All results are written to the params= dataset. The number of fitted gaussians does not have to be the same between the parinp= and params= datasets. This creation of an extra dataset is needed because MIRIAD very-deep-down disallows opening an existing dataset for writing. So it is not possible to add new fits to an existing parameter set.
If parinp= is present, but in= is not, then no new fits are made, but the gaussians in parinp are sorted as specified by the cmpsort= keyword. They can also be transformed as specified by options=integral,fwhm,dispersion.
images(z1,z2)Select image planes z1 to z2 inclusive. z2 defaults to z1.
quarter(z1,z2)Select the inner quarter of the image planes z1 to z2 inclusive. If both z1 and z2 are missing, then all planes are selected. If only z2 is omitted, z2 defaults to z1.
boxes(xmin,ymin,xmax,ymax)(z1,z2)Select the pixels within a box with corners xmin,ymin,xmax,ymax. z1 and z2 are the same as in the "image" subcommand. If z1 and z2 are omitted, all planes are selected.
polygon(x0,y0,x1,y1,x2,y2,...)(z1,z2)Select the pixels within the polygon defined by the list of vertices. z1 and z2 are the same as in the "image" subcommand. If z1 and z2 are missing, all planes are selected. If only z2 is omitted, it defaults to z1.
mask(file)Select pixels according to the mask given in the file.
The units of the numbers given in the above commands are, in general, absolute pixels. But this can be changed (and rechanged) by using one of the following subcommands.
abspixelCoordinates are interpreted as absolute pixel values, the default.
relpixelCoordinates are relative to the reference pixel of the map.
relcenterCoordinates are relative to the central pixel of the map, (defined as (naxis1/2+1,naxis2/2+1)).
arcsecCoordinates are in arcseconds, relative to the reference pixel.
kmsCoordinates in the third dimension are in km/s. Note: the region=mask option is not implemented. The mask of the input dataset is used however. If there is one, the profile value at masked datapoints is set to zero before doing the fit.
nofit: output the initial estimates, don't make fits findestim: let gaufit determine the initial estimates
noprint: do not print the fit results on the terminal supbad: suppress results for fits outside ranges given by cutoff, crange and wrange and results for bad fits. estimout: print initial estimates intermout: print some intermediate results for multi-component fits abspix: x, y coordinates on output are relative to lower left, rather than relative to crpix abscoo: x, y coordinates on output are absolute coordinates
wrprof: write out a file with the data and the fit so that it at least is possible to use plotting programs to compare them; a kludge until gaufit itself can plot. filenames will be 'profile_at_$x_$y' (or given by prof=)
integral: write out integral of gaussian instead of amplitude dispersion: write out dispersion of gaussian instead of fwhm pixel: write center and width in pixels, not in axis units (for these three: also interpret input for cutoff, vrange, wrange and estim keywords as int/disp/pix)
average: first make an average profile of the selected region and then fit one single gaussian to this profile summed: first make a summed profile of the selected region and then fit one single gaussian to this profile
negative amplitudes may be both positive and negative, instead if just positive fixvelo: fix the velocities to the initial estimate during fit fixwidth: fix the width to the initial estimate while fitting (fixvelo and fixwidth can be combined)
velocity, amplitude, integral, width, vdiff, vrangeOption 'velocity' and 'fwhm' result in components sorted on increasing velocity or width. Option 'amplitude' and 'integral' result in components sorted on decreasing amplitude or integral. If vdiff is used, then a second parameter gives a center velocity; components are sorted based on the difference between the fitted velocity and this center velocity. If vrange is used, the second and third parameter give a velocity range. If one component is within this range, it becomes the first. If none or more than one is within this range, they are sorted on velocity. Usually, cmpsort is applied for every pixel of the dataset. This is wanted when originally fitting (in= used). It is also generally wanted when refitting part of the dataset (in= and parinp= used), especially when more gaussians are to be added in selected regions. However, when only parinp= is present, the sorting is done only in the selected region and everything outside is left alone.