Program: gaufit2d Purpose: Find position or major-, minor axis and position angle of sources by fitting 2 dim. gauss function parameters to data in a box.. Category: MODELS, FITTING File: gaufit2d.c Author: M.G.R. Vogelaar (GUI: J.P. Terlouw) Description: CONTENTS: -PURPOSE -GRAPHICAL USER INTERFACE -DESCRIPTION OF THE MODEL -FIT PROCEDURE -WARNING GRID CONVERSIONS -LIST OF KEYWORDS/BUTTONS AND INPUT FIELDS -APPENDIX I -EXAMPLE OF OUTPUT IN LOG FILE PURPOSE: Program GAUFIT2d is a tool for finding the exact position of the maximum in a source or to find major & minor axes and position angle of a source by fitting a two dimensional gauss function to the data. Also the volume under the 2d-gauss is returned in units of the amplitude (see appendix 2). Output is a list of fitted parameters on screen or in the GIPSY Log file. The program works either in a coordinate system of grids or in a system of seconds of arc. If you press the button, an example run is shown. GRAPHICAL USER INTERFACE: The graphical interface is closely related to Hermes, i.e. you are still able to run the program from within Hermes and specify keywords on the command line. The GUI is a nicer way to present keywords and an easier way to enter values for these keywords. Keyword values entered on the command line will have an effect on the status of the GUI immediately. As an example start the program with something like this: gaufit2d inset=test f 0 box=-192 337 -182 346 Then the GUI starts with filled input fields for SET and BOX. Note that for input fields the same input rules apply as for keywords in Hermes DESCRIPTION OF THE MODEL: The function that this program uses to find model parameters is: F(x,y) = A * EXP[ -4*Ln(2)*( (xr/HW1)^2+(yr/HW2)^2 ) + Z0 where: xr = -x1 * sin(phi) - y1 * cos(phi) yr = x1 * cos(phi) - y1 * sin(phi) and: x1 = x - x0 y1 = y - y0 The parameters are fitted either in a grid system or an system in which both axis units are arcsec. The parameters are: A, amplitude in units of the image data. x0, center in X in grids. y0, center in Y in grids. HW1, first Full Width at Half Maximum in arcsec or grids HW1, second Full Width at Half Maximum in arcsec or grids phi, angle between X axis and first axis of ellipse, angle always between -45 to 45 degrees. Note that angles in the grid system can differ from angles in the arcsec system if the length of a grid in X is unequal to the length in Y. Z0 Offset of the data (zero level) in units of the image data. FIT PROCEDURE: The method used to fit these parameters is a least squares method. This method needs initial estimates. The fit can be very sensitive to these estimates, therefore you are able to overrule the default estimates set by the program and you are also able to fix one or more initial estimates to decrease the degrees of freedom and to improve the results. Estimates are calculated with first and second moments of the data. The parameters of a 2-dim gauss function can be expressed in these moments. For the calculations of initial estimates, a cutoff value is introduced. This value includes only data if the image value of a pixel is more than 5% of the maximum value in a given box. This maximum is in fact the average of the four pixels with the highest value. The least squares fit stops when successive iterations fail to produce a decrement in reduced chi-squared less than a certain value called the relative tolerance. Maximum accuracy can be obtained by setting this value to 0 which is the default in this program. If a fit needs more iterations than the maximum number, you can set the tolerance to e.g. 0.05 (or increase the maximum number of iterations). GAUFIT2D can do its calculations in either a (square) grid system or a sky system where both units are seconds of arc (from now on called the arcsec system). In the arcsec system a position is calculated by multi- plying the grid position with the grid spacing (which can be different in the x- and y direction) and both estimates- and fit routines do their calculations in the arcsec system. If the header information from a set allows grids converted to arcsec then the calculations are default in the arcsec system. This can be overruled by pressing the button (appears in: EXTRA KEYWORDS). WARNING GRID CONVERSIONS: The fit results (an ellipse with major, minor axis and position angle) are plotted if the button is on. The image data is always displayed in grids and therefore the fit is also in grids. This can result in confusing results because if the grid spacings in both x and y are not equal, the position angle is different in the arcsec and grid system. LIST OF KEYWORDS/BUTTONS AND INPUT FIELDS: Keywords: Each keyword corresponds to a button or input field! EXAMPLE= (Button: ) Run an example. It creates a simple set called "GAUFIT2Dexampleset". It fits parameters and shows the result in GIDS. INSET= (Input field: SET) Input set and optionally subsets. The subset must be two dimensional. Model parameters are converted to arcsec if the input is a spatial map, i.e. the first axis is a spatial longitude axis and the second is a spatial latitude axis. If you want the results in grids in any case use the button under 'Extra keywords'. Maximum number of subsets is 1024. You can change the current subset level by pressing one of the '<' or '>' buttons or specify a new level in the 'subset' input field. Examples: INSET=AURORA FREQ 10:20 INSET=AURORA FREQ INSET=AURORA FREQ 8 10 15:18 BOX= (Input field: BOX) Enter the limits in grid coordinates of the box from which you want to extract data for the fit. Examples: BOX=-20 -10 20 10 BOX=0 0 D 10 10 center at 0,0 size: 10x10 SUBSET= (Input field: SUBSET) NOTE: This field can only be used AFTER you specified subsets in the SET input field. Then you can change the value in this field to force reading data from another subset. Subsets can also be changed with the ">" and "<" buttons to the right of the input field. AMPLITUDE= (Input field: AMPLITUDE) Estimate for the value of the central value in least squares fit. Units are the same as the units of your data. This parameter can be adjusted to be a better estimate in your opinion or can be fixed with the / button. A fixed value will always re-appear in the column with fitted parameters. MAJOR= (Input field: FWHM X) Estimate for the value of the major axis in grids or in arcsec. The value is in arcsec if a conversion from grids to arcsec is possible and the button is off. Otherwise the input must be in grids. MINOR= (Input field: FWHM Y) Estimate for the value of the minor axis in grids or in arcsec. CENTERX0= (Input field: CENTER X0) Estimate in grids for the X-position of the maximum of the source. CENTERY0= (Input field: CENTER Y0) Estimate in grids for the Y-position of the maximum of the source. ANGLE= (Input field: ANGLE) Call the axes of an unrotated frame the X and Y axis, and the main axes of a rotated ellipse X' and Y'. Then ANGLE= is an estimate for the smallest angle (positive or negative) between the X and X' (or Y and Y') axis. The range is always from -45 deg to 45 deg. For the arcsec system this angle is converted to a real position angle, i.e. an angle between the major axis and the direction of the North. ZEROLEVEL= (Input field: ZERO LEVEL) Estimate for the offset of the two dimensional gauss. In some cases (e.g. antenna patterns) this parameter can be set to zero and fixed. All these parameter keywords also have prefixes F_ for the corresponding fixed/free toggle used in the lsq.fit and R_ to restore the previously calculated moments estimate. QUIT= (Button: ) Terminate program. FIT= (Button: ) Start a least squares fit using the estimates that were defined before pressing this button. NEXT= (Button: >) Go to the next subset in the sequence. Works only if more than one subset was entered. PREV= (Button: <) Go to the previous subset in the sequence. LOG= (Button: ) Send all fit information to the Hermes Log and Screen file. The button is only active if there is a fit result. Some buttons/input fields appear when you press the button: TOL= (Input field: ) Fractional tolerance for the chi square (chi2) default set to 0.0. Fitting is stopped when successive iterations fail to produce a decrement in reduced chi-squared less than TOL= The value cannot be less than a certain minimum as set by the system. This means that maximum accuracy can be obtained with TOL=0.0 LAB= (Input field: ) Value for mixing parameter default set to 0.001. Mixing parameter in the least squares fit function. LAB= determines the initial weight of steepest descent method relative to the Taylor method. LAB= should be a small value (e.g. 0.001). CLIP= (Input field: CLIP) Value below which data is excluded in the least squares fit. Default the clip is not set, meaning that all data contributes to the fit. If you image data is very noisy then a better fit could be made if you set a clip near the value of the zero level that you expect. MAXITS= (Input field: MAXITS) Maximum number of iterations allowed in the least squares fit routine. The default is 100. GRIDS= (Button: ) If On: All calculations are in grids. If off: If header information allows conversion from grids to arcsec, then calculations are in arcsec. Display: VIEW= (Button: ) If necessary start GIDS. Display current subset and ellipse representing the major,minor axis with position angle. ERASE= (Button: ) Erase current plotted ellipse in GIDS and restore image. RESIDUE= (Button: ) Calculate difference between data and model and display the result in GIDS. The button is only active if a set is viewed and a fit is made. The residue is stored in a set on disk. The name of the set is "XXXXresidualXXXX". The set is removed after the program has finished. Note that the program only shows the residual in the current box. Notes: APPENDIX 1: A note about chi2. The quantity chi-squared is defined as: chi2 = Sum[ (Y_i-Ymodel_i)^2/sigma_i^2 ] and the reduced chi-squared is: reduced-chi2 = chi2 / (N-n-1) where N is the number of data pints in the fit and n the number of parameters. However, the values of sigma_i (sigma per pixel) is not known and set to 1. Then the reduced chi2 is equal to the variance of the fit (Bevington, Data Reduction and Error Analysis for the Physical Sciences, p 188) and chi2 is not a measure of the goodness of fit anymore. So we prefer to list only the variance of the fit. If you assume the fit is a good approximation of the model then reduced-chi2 ~ 1 and because reduced-chi2 = variance / sigma_av you have an approximate value for sigma_av, the weighted average of the individual variances. APPENDIX 2: F(x,y) = A * EXP[ -4*Ln(2)*( (xr/HW1)^2+(yr/HW2)^2 ) + Z0 The volume V under a 2d-gauss limited by the offset ZO is given by: V = 2.A.HW1.HW2.PI.Erf(Inf).Erf(Inf)/Ln(256) = 1.13309.A.HW1.HW2 This volume is expressed in units of the amplitude, therefore it will be the same number while working in grids or physical coordinates. Example: EXAMPLE OF OUTPUT IN LOG FILE: ===================================== GAUFIT2D ================================ Date of fit : 03-Mar-1998 (17:06:38) Name of set : [test] Fit box : [-192 337 -182 346] Iterations : 8 Variance of fit : 0.192932 Points in box : 110 # excluded by clip : 0 Number of blanks : 0 Volume under gauss : 568.882 in units of the amplitude Coordinate system : Arcsec Size of 1 grid : 4.32 x 5.38978 arcsec Center X: -187.34 (grids) = 179.1631786 (DEGREE) = 11h56m 39.1s Center Y: 341.41 (grids) = 54.1635537 (DEGREE) = 54d 9m48.79s Major axis: 19.54'', Minor axis: 15.66'', Position angle: 155.3 (deg) =============================================================================== Updates: Feb 26, 1998: VOG, Document created. Jul 13, 2005: VOG, Repaired bug in QUIT= keyword handler (avoid use of freefmatrix with null pointer) Apr 14, 2009: VOG, Removed unused macro for NINT() Apr 17, 2010: VOG, Extra fit information in LOG