Program: FINDGAUSS Purpose: Fitting and subtraction of 2-D gaussians Category: FITTING, MODELS, TABLES File: findgauss.c Author: M. Vogelaar Keywords: INSET= Give set, subsets: Maximum number of subsets is 2048. BOX= Give box in ..... [entire subset] ** CUTOFF= Give cutoff ratio for gaussian [1/100] This is the ratio of the gaussian function values at central point and a point x on the boundary. OUTSET= Give output set (, subsets): Output set and subset(s) for the result. The number of output subsets is the same as the number of input sub- sets. SUBFILE= Use parameter file to subtract gaussians? Y/[N] Instead of calculating gaussians, it is possible to use parameters from a previous session to subtract the fitted gaussians. The parameters are written to an ASCII file and can be edited before the actual sub- traction. If SUBFILE=Y: PARAMFILE= Name of file with stored parameters: [Stop program] File must be created in a previous run of this program. Only a subtraction on INSET= is performed. If SUBFILE=N: PARAMFILE= Name of file to store parameters: [No file] Write gauss parameters to file so that these parameters can be edited and used for subtraction. If a given file name already exists, APPEND= must be specified. APPEND= File exists, ok to append? [Y]/N A specified file name already exists. You can append to this file with APPEND=Y. If APPEND=N you will be prompted for another filename. Pre specification of initial parameters in the fit: AMPLITUDE= Amplitude (map units): [Calculated] BEAM= Beam: major, minor: [Calculated] If indicated, in seconds of arc. X0Y0= Centre in grids: [Calculated] ANGLE= Angle of major axis wrt pos. Y-axis (deg): [Calculated] ZERO= Zero level (map units): [Calculated] MASK= Mask (0=fixed,1=free): [1,1,1,1,1,1,1] FITSIZE= Give size (x*y) of fit box: [Calculated] Give a fixed size for the 'fitbox' instead of determination of suitable sizes by the program. SIGMA= Rms of noise in map: [None] VALRANGE= Range of values in search area: [4*sigma, ->] Or if SIGMA= was not specified: Range of values in search area: [All data] Data values outside given range cannot be candidates for a gauss fit. ** TOLERANCE= Convergence criterion. [0.01] Fitting stops if successive iterations fail to produce a decrement in reduced chi-squared less than TOLERANCE= ** LAB= Mixing parameter: [0.01] Mixing parameter, LAB determines the initial weight of steepest descent method relative to the Taylor method. LAB should be a small value (i.e. 0.01). ** CLIPLOHI= Give clip levels: [Include all data] Examples: 1) CLIPLOHI= include all data 2) CLIPLOHI=3 include all data >= 3 3) CLIPLOHI=3 5 include all data >=3 and <= 5 Subtract only gaussians with amplitudes and beams in range AMPRANGE= and BEAMRANGE=: AMPRANGE= Amplitude range to subtract: [All amplitudes] BEAMRANGE= Beam range to subtract: [All values] TABNAME= Name of the table to store table data: [GAUSS] Description: The program FINDGAUSS fits 2 dimensional gaussians to sources in an image. The fitting is used for example to get rid of point sources in a map or to find positi- on and flux of these point sources. The set that is examined is INSET=. The defined subset must be two dimensional. The number of subsets can be greater than 1 but cannot exceed 2048. The items corresponding with the grid spacings (CDELTn) must be available from the the header of the input set, otherwise the program will abort. If it can find the grid spacing, and the units of the spacing is in degrees (along both subset axes), then input will be done in seconds of arc else it will be done in units as found in the header. With BOX= the area that will be searched for gaussians can be made smaller. The default is the entire (sub)set. OUTSET= defines the name of the set where the result is stored. This result is always the subtracted INSET= and therefore can be used as an inspection set. If you want to write to an already existing set, you will be warned with the message: Will overwrite data, okay ? [Y] and at the OKAY= prompt. The application tries to allocate memory for your image and when it is not able to allocate enough memory, it aborts. Then you can try again with a smaller box. If the result of the subtractions is not exactly what you wanted or expected because a couple of unwanted fits were extracted, there is no need to rerun the application completely. It is possible to write the gauss parameters to an ASCII file with PARAMFILE= This file can be edited. To use it for subtraction, you must give SUBFILE=Y first. Then you are prompted with PARAMFILE= Pressing carriage return will stop the program. If a file name is given and the file can be opened, the parameters in this file are used to define the gaussians. The parameter file is made in a run of FINDGAUSS with SUBFILE=N For each subset a file name can be specified. For each fit, the file is extended with one line with 12 items: 1) An integer 0 or 1 (0: not selected for subtraction) (1: selected for subtraction) 2) Value of fitted amplitude in map units. 3) Major axis in (if possible) seconds of arc. 4) Minor axis in (if possible) seconds of arc. 5) Centre X coordinate in grids. 6) Centre Y coordinate in grids. 7) Angle between major axis and positive y-axis in degrees (In most cases the pos. y-direction is also the direction of the North). 8) Zero level in map units. 9) Lower X value of (fit)box where gauss is fitted. 10) Lower Y value. 11) Upper X value. 12) Upper Y value. The fitting: FINDGAUSS finds the maximum in a subset by examining data with values above a certain level or between two levels (VALRANGE=) If SIGMA= has a value, the default for VALRANGE= is one (lower) clip equal to 4 times the value in SIGMA= If SIGMA= was not specified, all data in a subset is used, but this will result in attempts to fit gaussians with an amplitude comparable to the noise, so it is better to define a threshold. One value, say x will result in a value range starting with x and including all values greater than x. Two values, say x,y indicate a range where the search area is determined by all values between x and y (x, y values included). As mentioned, the fit parameters are: 1) Amplitude, 2) Major axis, 3) Minor axis, 4) Centre X, 5) Centre Y, 6) Position angle of major axis, 7) Zero level. In the default situation, all parameters are free in the fit (i.e. after the fit they can have other values than the pre specified or estimated values). You can fix parameters in the fit with the keyword MASK= The default is MASK=1 1 1 1 1 1 1 (All free). If MASK=1 0 0 1 1 1 1 you fix the major and minor axis etc. Fixed parameters must have an initial value. This can be a value estimated by the program, or it can be a pre specified value. You can give pre specified values for fixed and free parameters (which will in both cases overrule the calculated estimates). Pre specifications are given with the keywords: AMPLITUDE= amplitude in map units, ANGLE= position angle in degrees (angle between major axis and positive Y-axis counted counter clockwise). You will be warned if the pos. Y-axis does not align with the north in your map. BEAM= major, minor axes (in units as prompted). X0Y0= centre X, Y in grids. ZERO= zero level in map units. The fitting can be influenced by two variables given in TOLERANCE= and LAB= TOLERANCE= is the relative tolerance. Fitting of the gaussian function stops when successive iterations fail to produce a decrement in reduced chi-squared less than TOLERANCE=. If its value is less than the minimum tolerance possible, it will be set to this value. This means that maximum accuracy can be obtained by setting TOLERANCE=0.0. LAB= is the mixing parameter, LAB= determines the initial weight of steepest descent method relative to the Taylor method. LAB= should be a small value (i.e. 0.01). The data that is used in a fit can be limited with the keyword CLIPLOHI= If you give no values (default) than all data in the 'fitbox' is used in the fit. If you give one value, then all data greater than or equal to that value will be included and if you give two values, then data with values greater than (or equal to) the first clip and less than (or equal to) the second clip are included in the fit. The decision to subtract or not is limited by the key- words AMPRANGE= and BEAMRANGE= Default is the subtraction of all fitted gaussians. If you give one value for a range keyword, then all amplitudes (or beams) with a value greater than this value will be included in the subtraction. If you give two values, then all amplitudes (or beams) with values greater than (or equal to) the first range value and less than (or equal to) the second value will be included in the subtraction. Finding a 'fitbox': If the maximum in a subset is found, a box is created by appending a 'border' around this pixel. This box has size 3x3. Then the box is extended with a new surrounding border. If the maximum in the new extension exceeds the maximum in the previous extension, the first box becomes the 'fitbox'. This is also true if the mean in the new extension exceeds the mean in the previous extension or if the decrease in the mean is less than SIGMA= (i.e. if SIGMA= was specified). The last criterion is that the mean becomes smaller than the maximum times a cutoff factor given in CUTOFF= The default is 1/100. You have to experiment with the keywords SIGMA= and CUTOFF= to get an idea how the growth of a fitbox is stopped. Fixing parameters: You have the possibility to combine the fixing of parameters with the pre specification of them. Good results can be obtained by pre specifying the beam (BEAM=) and angle (ANGLE=) but with all parameters free in the fit (MASK=1 1 1 1 1 1 1) If, for instance, you fix beam and angle, then more gaussians will be found, but after the subtraction you will also find more residuals because not all sources have the pre specified properties, but just gave a 'good' fit. Use of tables: The results of this application are stored in a table. The default name of the table is GAUSS and the columns are AMPLITUDE, MINOR, MAJOR, X0grid, Y0grid, X0phys, Y0phys (physical counterparts of X0grid, Y0grid), ANGLE, ZEROLEV, AMPerr, MAJerr, MINerr, X0err, Y0err, ANGerr, ZEROerr, XLO, YLO, XHI, YHI (4 box coordinates in grids). The name of the table can be changed, the name of the columns not. If one table name is used twice, the contents of the last table will overwrite the first one. Notes: ....... Example: Example of a default file: AMPLITUDE= / Program calculates estimate AMPRANGE= 0.04 12 / Use in subtraction only fits with ampl. in this range ANGLE= 0 / Pre specify angle, 0 degrees wrt north (to east) BEAM= 19 13 / Pre specify beam, major axis axis 19 arcsec, minor 13 arcsec. BEAMRANGE= 5 40 / Use in subtraction only fits with beams in this range BOX= / Operate on entire subset CUTOFF= 1/100 / Limit 'fitbox' to boundaries with mean > cutoff*maxampl. PARAMFILE=storeparm.txt / Sore fit parameters in ASCII file FITSIZE= / Do not use a fixed size for the fitbox. INSET=AURORA freq 0 / Set and subset MASK= 1 1 1 1 1 1 1 / All parameters free in fit OUTSET=SUBAURORA / Name of subtracted data SIGMA= 0.01 / Rms of noise SUBFILE=N / Do not use an ASCII file with gauss parameters for subtraction TABNAME= / Use default name for table VALRANGE= / Sigma was specified so default is VALRANGE=4*0.01 which implies that only data > 0.04 can be candidates for a fit X0Y0= / Program calculates estimate ZERO= / Program calculates estimate Updates: Sep 22, 1992: VOG, Document created. Oct 4, 1993: VOG, Repaired bug in preset of FWHM Mar,25 1994: VOG, Replaced some 'tofchar's by str2char in routine 'putintable'.