.EQ delim $$ .EN .RP .TL Specifications for the Aperture Photometry Package .AU Lindsey Davis .AI .K2 "" "" "*" Revised October 1987 .AB The APPHOT package is a set of tasks for performing aperture photometry on uncrowded or moderately crowded fields, in either interactive or batch mode. The photometric technique employed is fractional pixel aperture integration. Point spread function fitting techniques are not used and no knowledge of the point spread function is required for the computation of magnitudes. Separate tasks are provided for creating and modifying object lists, computing accurate centers for a list of objects, computing sky values for a list of objects and performing aperture photometry inside circular or polygonal apertures. .AE .NH Introduction .PP The APPHOT package is a set of tasks for performing aperture photometry on uncrowded or moderately crowded fields. The photometric technique employed is fractional pixel integration. Point spread function techniques are not used and no knowledge of the point spread function is required for the computation of magnitudes. .PP The APPHOT package performs multi-aperture photometry on digitized starfields maintained as IRAF image files. Input to the package consists of an image file, a list of object coordinates, numerous parameters controlling the analysis algorithms and, optionally, the graphics terminal and/or display. APPHOT output consists of successive lines of text where each line summarizes the results of the analysis for a particular object. .PP Given starting coordinates and an IRAF image, the principal APPHOT task computes accurate centers, sky values and magnitudes for a list of objects. Separate IRAF callable tasks in APPHOT exist to create and modify object lists, to determine image characteristics such as the full width half maximum of the point spread function or standard deviation of the sky pixels, to compute accurate centers for a list of objects, to compute accurate local sky values for a list of objects, and to compute magnitudes inside a polygonal aperture. .NH APPHOT Package Requirements .NH 2 APPHOT Package Input .NH 3 The IRAF Image .PP APPHOT assumes that images exist on disk in IRAF format. Facilities for reading and writing images exist elsewhere in IRAF such as in package DATAIO or MTLOCAL. .PP APPHOT assumes that the image data is linear. The input image is assumed to have been corrected for those instrumental defects which affect the intensity of a pixel. These include pixel to pixel variations in the bias values, pixel to pixel gain variations, cosmic rays, cosmetic defects, geometric distortion, and detector non-linearities. .PP APPHOT assumes that the local sky background is approximately flat in the vicinity of the object being measured. This assumption is equivalent to assuming that the local sky region has a unique mode. Therefore any variations in the sky background which occur at the same scale as the sky region should be removed prior to entering APPHOT. .PP The point spread function of the image is assumed to be constant for all regions of the image. This is particularly critical in the case of faint objects for which small apertures which minimize the effects of crowding and sky noise in the aperture are used. The wings of the object will almost certainly extend beyond the aperture and good results will only be obtained if objects are consistently well centered and the shape and diameter of an object is constant throughout the image and invariant to magnitude. .NH 3 The Coordinate Lists .PP All APPHOT tasks operate on lists of object coordinates or interactive cursor input. Lists are maintained as text files, one object per line with the x and y coordinates in columns one and two respectively. List files may be created interactively with either the graphics or image cursor, by a previously executed APPHOT task, by a previously executed IRAF task or by a user task. .PP The x and y coordinates contained in the list file can be either the starting centers for the objects of interest or the true centers. In either case the aperture centered around the object position must not extend beyond the boundary of the image. To obtain reliable results from the centering algorithm the starting centers should not be more than 1 fwhmpsf pixels from the true centers. .NH 3 Algorithm Parameters .PP Many tasks in the APPHOT package have an associated set of algorithm parameters which control the operation of the analysis routines. For example the centering and sky fitting algorithms each have an associated group of parameters. The APPHOT tasks supply reasonable defaults for these parameters. However the user may alter them at any time in interactive mode or by editing the appropriate parameter file before running the task in batch. .NH 3 Terminal Graphics and the Image Display .PP Most APPHOT programs may be run in either interactive or batch mode. In interactive mode, the user can mark the positions of objects by interactively, positioning the image cursor on the display device. Simple cursor commands can interactively alter the algorithm parameters. In batch mode the algorithm parameters are fixed and object coordinates are taken from the user supplied coordinate list. The display device is not required in batch mode. .PP At present there is no high level IRAF display interface. Therefore the present version of APPHOT replaces the display device with terminal graphics. For example it is possible to load a contour plot of an image or image section, and run the APPHOT tasks, interactively marking positions on the plot. This is a temporary measure to tide thing over until the display interface is available. Furthermore this option is only really suitable for those terminal which have both text and graphics windows displayed at the same time. .PP Certain APPHOT tasks such as FITSKY occasionally require an interactive graphics capability. For example it is sometimes desirable to plot the histogram of the sky pixels and mark the peak of the histogram interactively rather than using the automatic sky fitting routines. Graphics capablity has been added to APPHOT tasks as deemed useful. .NH 2 APPHOT Package Functions .NH 3 Creating Coordinate Lists .PP APPHOT task(s) shall exist to create coordinate lists interactively within APPHOT by loading the IRAF image into the image display and successively marking the objects to be measured with the image cursor. It shall be possible to mark the positions of objects on the display, and draw circles of arbitrary size around the objects of interest. It shall be possible to verify existing coordinate lists by marking the object coordinates on the image display. .PP At present the full image display cababilities do not exist in IRAF. Therefore as a temporary measure most tasks will run with a contour plot and the graphics cursor as a substitute for the image display and image cursor by setting the image cursor to stdgraph and the display parameter for each task to stdgraph. The output coordinates will be written to a text file and may be used as starting coordinates for the APPHOT routines. .PP Those sites which have SUN workstations can use the IMTOOL facilites to create cursor lists. .PP APPHOT tasks shall exist to automatically detect objects by specifying the IRAF image to be searched, plus a few image characteristics such as a detection threshold and the full width half maximum of the point spread function. The output coordinates plus a few object statistics will be written to a text file. .NH 3 Coordinate List Operations .PP General list processing tasks are or will be available in the IRAF lists package. Examples of useful currently existing tasks are list sorting by for example magnitude and performing a linear transformations on coordinate lists. .PP It may eventually be necessary to add some specialized list processing tasks to APPHOT. One example is a list matching facility in which objects in common to two lists are combined and output to a third list. .NH 3 Determining the Image Characteristics .PP APPHOT tasks shall exist to estimate certain image characteristics such as the full width half maximum of the image point spread function and the standard deviation of the background. In order to obtain meaningful error estimates it is also necessary to specify the noise charactersitics of the detector and a noise model. In this version of APPHOT two noise functions will be supported a constant sky background noise, and constant background noise plus Poisson statistics in the object. .PP The reason for this capability is that the majority of the APPHOT algorithm parameters scale with these two image characteristics. For example all the pixel scale dependent parameters such as the centering aperture, object apertures, the two sky annuli and the maximum shift parameter all scale with the full width half maximum of the point spread function. Similarly most of the pixel intensity dependent parameters scale with the standard deviation of the sky pixels. .NH 3 Centering .PP An APPHOT task shall exist to determine accurate centers for a list of objects, using the approximate object coordinates as a starting point. The centering aperture may extend beyond the boundaries of the image but a warning message will be generated. .PP A choice of centering algorithms with a range of efficiency and accuracy will be available. The center determination must be resistant to the affects of nearby contaminating objects. .PP The user may opt to use the starting center as the actual center in all cases, to use the starting center as the actual center if the object is very faint, or tweak up the center if the signal to noise is above a certain threshold. .PP The output of the centering algorithm will the set of computed coordinates and their errors. .NH 3 Fitting the Sky .PP An APPHOT task shall exist to compute a constant background value by analysis of an annular region surrounding the object. The background is assumed to be flat in the region of the object, but may contain contaminating objects or defects which shall be detected and eliminated by the fitting algorithm. It shall be permissible for the background region to extend beyond the boundaries of the image; the out of bounds region of the background shall be excluded from the fit. .PP The user may supply a constant background, with zero-valued noise and skew value or read appropriate values from a file instead of computing the background. .PP The output of the sky fitting algorithm shall be the mode, standard deviation and skew of the sky as well as the number of pixels left in the background region after pixel rejection and the number rejected. .NH 3 Multi-aperture Photometry .PP APPHOT routines shall exist to compute the integral of object minus background within one or more circular apertures centered upon the object. The integration shall be performed using partial pixel techniques, to minimize the effects of sampling. If the aperture contains any indefinite pixels, or if the aperture extends beyond the boundary of the image, an indefinite result shall be returned. Both normal and weighted photometry routines shall be available. .PP It shall be possible to normalize the output magnitudes to a user specified integration time using values stored in the IRAF image header. .PP The output shall consist of magnitudes and errors for each of the specified apertures. .NH 3 Polygonal Aperture Photometry .PP Determine the integral of the object within a specified polygon. Aperture integration shall be by fractional pixel summation of successive image lines. .NH Specifications .NH 2 Apphot CL Callable Tasks .PP The CL callable part of the APPHOT package consists of the following routines:\f(CW .nf .na mark -- create/verify coordinate lists interactively polymark -- create/verify polygon lists interactively daofind -- automatic image finder from DAOPHOT lintran -- linearly transform a list of coordinates (LISTS package) radprof -- compute the radial profile of an object. fitpsf -- model the PSF center -- compute accurate centers for a list of objects fitsky -- compute accurate sky values for a list of objects phot -- perform multi-aperture photometry on a set of objects polyphot -- compute the flux inside a set of irregular polygons wphot -- perform weighted multi-aperture photometry\fR .ad .fi .NH 2 Standard Analysis Procedures .PP A standard reduction procedure would be something as follows. .IP 1. Estimate the image and data characteristics in the following manner. .RS .IP 1. Use the interactive setup option (i key) in the PHOT task to specify the fwhm of the psf, the centering apperture, the sky annuli, the photometry apertures as well as estimate the sky sigma. Examine the results (:show command) and optionally store the new parameters in the parameter files (w key). .IP 2. Use the EPAR task to edit the parameter files directly. For example the detector gain and readout noise must be known before running PHOT in order to get reliable error estimates. .RE .IP 2. Prepare a coodinate list in one of the following ways. .RS .IP 1. Running MARK or POLYMARK. .IP 2. Run DAOFIND or some other automatic image finder. .IP 3. Run another APPHOT task such as CENTER . .IP 4. Run any other IRAF or user program which generates the appropriate format. .IP 5. Transform an existing list using the LISTS package facilities .RE .IP 3. Run PHOT or WHOT in non-interactive mode to perform the multi-aperture photometry. The user should be familiar with the algorithms used to measure the object center, to fit the background, and compute the aperture integral, magnitude, and errors. The values of all visible and hidden PHOT parameters should be inspected before doing any serious processing. Note that CENTER and FITSKY can be run independently of PHOT and their output used as input to PHOT. .EQ delim $$ .EN .NH 2 The APPHOT Algorithms .PP The principal APPHOT algorithms are described below. .NH 3 The RADPROF Algorithm .PP The function of the RADPROF task is to compute the radial profile of selected, isolated, high signal to noise objects in an IRAF image. .PP Given the IRAF image, an initial center, an aperture and a step size, RADPROF computes the radial profile of an object in the following manner. .IP 1. The APPHOT centering routines are used to compute an accurate center for the object. The new center is used to compute the radial coordinate for each pixel in the subraster. .IP 2. The APPHOT sky fitting routines are used to compute an accurate sky value for the object. The new sky value is used to compute the object intensity for each pixel in the subraster. .IP 3. The CURFIT package routines are used to produce a smooth intensity versus radius curve by fitting a cubic least squares spline function to the data using linear least squares techniques. .IP 4. The fitted curve is evaluated at equally spaced intervals to produce the final radial profile. The profile is normalised using the r = 0 value. .IP 5. The smoothed curve is integrated using the 1D integration routines IMINTERP to evalute the fraction of the total integral at each interval. .IP 6. The three quantities $I sub r$ the raw intensity, $Inorm sub r$ the normalized intensity, and $Ifraction sub r$ the fraction of the total intensity inside radius at r, are plotted on the terminal and output as a function of r to a text file. .PP RADPROF is used principally to compute the stellar image and setup the parameters for doing photometry. Its default output is a radial profile plot with the phot parameters marked on it. .NH 3 The FITPSF Algorithm .PP The function of the FITPSF task is to fit an analytic model to an object in order to derive information about the point spread function. Given an IRAF image, an initial center and an aperture, FITPSF computes the model parameters in the following manner. .IP 1. The fitting function is chosen. The present version of FITPSF offers two function choices a 2D radial Gaussian function and an elliptical Gaussian function. A third option, moment analysis, has been added to the FITPSF package. The image characteristics are evaluted by using the 0th, 1st and 2nd order moments of the image. .IP 2 Initial guesses for the parameters are made from the subraster data. .IP 3 The parameters and their errors are derived using standard non-linear least squares techniques and the NLFIT package. .NH 3 The DAOFIND Algorithm .PP The function of the DAOFIND task is to automatically detect objects in an IRAF image above a certain intensity threshold. The user specifies an intensiyt threshold, the estimated full width half maximum of the point spread function, and limits on the sharpness and roundness of the objects to be detected. DAOFIND produces a list of x and y coordinates and a sharpness and roundness statistic for each detected object. .PP The DAOFIND algorithm works in the following way. .IP 1. The user specifies an intensity threshold above local sky and estimates the full width half maximum of the point spread function. .IP 2. DAOFIND computes a convolution kernel which is a linear function of the brightness values in an approximately circular array of pixels. The original data is convolved with the computed kernel. The equivalent mathematical operation is a convolution of the original data with truncated, lowered circular Gaussian function with the specified FWHM, computed in such a way as to be the mathematical equivalent of fitting a Gaussian stellar profile to the object data by least squares. The coefficients of the linear function sum to zero so that the overall bias level (local sky) cancels out exactly. Since the function is symmetric about a single pixel, a smooth gradient in the sky brightness also cancels exactly. Therefore the user does not have to specify an absolute brightness threshold but only a threshold above local sky. .IP 3 The convolved data is searched for local maxima. Local maxima which are less than the specified threshold are rejected. A pixel is defined to be a local maximum when it is the brightest pixel within $n ~*~ sigma sub Gauss$. .IP 4. Using the original data, DAOFIND computes an object sharpness statistic which can be used to reject hot or cold pixels. The sharp parameter is defined below. Detected objects with a sharpness parameter outside the range specifed by sharplo and sharphi are rejected. Typical values for sharplo and sharphi are 0.2 and 1.0 respectively. .nf $I sub delta$ = height of best fitting $delta$ function = pixel value $I sub Gauss$ = height of best fitting Gaussian function = convolved value sharplo <= sharp = ${I sub delta} over {I sub Gauss}$ <= sharphi .fi If the brightness enhancement detected in the convolved data is due to a single bright pixel, then the best fitting delta function will have an amplitude equal to the height of this pixel above the local background, while the amplitude of the best fitting Gaussian will be pulled down by the surrounding low valued pixels. Sharp will be a number significantly greater than one. A single cold pixel will produce brightness enhancements approximately 0.5 * FWHM away from the pixel in question in all directions. In this case the height of the delta function which best fits these spurious maxima is close to zero, while the height of the best fitting Gaussian is some small positive number. In this case sharp will be close to zero. .IP 5. The roundess statistic is computed from the original picture data by fitting 1D Gaussians to the marginal sums in x and y. If the height of either 1D Gaussian is negative the object is rejected. If both heights are positive the roundness parameter is computed. .nf $DELTA I ~=~ I sub {x ~ Gauss} ~-~ I sub {y ~ Gauss}$ $ ~=~ I sub {x ~ Gauss} ~+~ I sub {y ~ Gauss}$ $roundlo ~<=~ round ~=~ {DELTA I} over ~<=~ roundhi$ .fi Objects which are elongated in the x direction have roundness values less than zero and those elongated in the y direction have roundness greater than zero. Round discriminates only against objects which are elongated along either rows or columns. Objects elongated at a 45 degree angle will not be rejected. .IP 6. Finally the approximate x and y centroids of the detected objects, rough magnitudes relative to threshold are computed, and object is added to the output file. .NH 3 The CENTER Task .NH 4 Centering Package Routines .PP The following are the principal programmer callable routines routines in the centering package. .nf \f(CW apcinit (ap, cfunction, cbox, fwhmpsf, noise) ier = apfitcenter (ap, im, xinit, yinit) ier = aprefitcenter (ap) value = apstat[ir] (ap, param) apstats (ap, param, str, maxch) apset[sir] (ap, param, value) apcfree (ap) \fR .fi .PP The following quantities can be examined or set by apstat/apset calls .nf \f(CW data 1. fwhmpsf full width half maximum of psf in pixels parameters 2. threshold minimum intensity for centering 3. noise noise model 4. sigma standard deviation of background 5. readnoise readout noise of CCD in electrons 6. epadu electrons per adu centering 1. calgorithm centering algorithm parameters 2. positive emission / absorption features 3. cbox centering radius in fwhmpsf 4. cmaxiter maximum number of fitting iterations 5. maxshift maximum permitted shift in fwhmpsf 6. minsnratio minimum permitted signal to noise ratio 7. clean clean subraster before centering 8. rclean cleaning radius in fwhmpsf 9. rclip clipping radius in fwhmpsf 10.kclean k-sigma rejection for cleaning in sigma \fR .fi .PP The following computed quantities can be examined by apstat calls. .nf \f(CW 1. xcenter computed x coordinate 2. ycenter computed y coordinate 3. xerr computed x coordinate error 4. yerr computed y coordinate error \fR .fi .PP See the manual pages for CENTER and CENTERPARS for a detailed description of the centering parameters. .NH 4 The Centering Algorithm .PP The function of the CENTER task is to compute accurate centers and errors for objects in an IRAF image given a list of initial positions and a centering aperture. .IP 1. If the centering algorithm is disabled, the computed center is set to the initial center and the uncertainty estimate is set to zero. .IP 2. If the cleaning algorithm is enabled, clean the subraster before centering. .IP 3. Estimate is made of the signal to noise of the object. If this quantity is less than a certain minimum value the computed center is kept but a warning message is issued. .IP 4. The center and errors are computed using one of several algorithms. .IP 5. If the computed center is greater than a user specified distance from the initial center then the computed center is returned with an error message. .NH 4 Symmetry Clean Algorithm .PP The symmetry-clean algorithm attempts to remove defects in the centering subraster by assuming that the object has radial symmetry and comparing pixels on the opposite sides of the center of symmetry. The algorithm works in the following way. .IP 1. The center of symmetry is computed. Normally the center of symmetry is assumed to coincide with the position of the brightest pixel in the subarray. However if the maximum pixel is more than a user specified distance away from the intial center, the initial center is used as the center of symmetry. .IP 2. Pixels inside a specified cleaning radius are unaltered. .IP 3. Pairs of pixels diametrically opposed about the center in the cleaning region between the cleaning and clipping radii are tested for equality. If the difference between the pixels is greater than a specified k-sigma rejection limit, the larger value is replaced by the smaller. In this region sigma is computed from Poisson statistics. .IP 4. Pairs of pixels in the clippping region are compared in the same manner as those in the cleaning region except that sigma is the standard deviation of the sky pixels. .PP The effect of the symmetry-clean algorithm is to edit the raster, removing any contaminating objects in the vicinity of the primary object. This simplifies the fitting algorithm and increases its reliability, since it does not have to deal with multipeak marginal distributions. .NH 4 Signal to Noise Estimate .PP The signal to noise of the object is estimated from the data values in the subraster in the following way. .nf $SNR ~=~ N sub object over {sqrt {n sub pix ~*~sigma sub sky sup 2}}$ or $SNR ~=~ N sub object over {sqrt {N sub object ~+~ n sub pix~*~sigma sub sky sup 2}}$ .fi where $N sub *$ is the number of counts in the object above threshold, $sigma sub sky$ is the standard deviation of the pixels in the sky region and $n sub pix$ is the number of pixels in the object aperture. The first approximation corresponds to constant sky noise only, the second includes Poisson noise in the object. .NH 4 Centroid .PP The centroid centering algorithm is similar to that used in MPC and can be briefly described as follows. For more detailed description see (Stellar Magnitudes From Digital Pictures, 1980, Adams et al). .IP 1. The marginal distributions in x and y are accumulated. .IP 2. The intensity weighted centroid positions for the marginals is computed using only data pixels which are above the threshold intensity. If the threshold parameter is 0 the mean intensity of the marginal is used in place of the threshold .nf $I sub i ~=~ I sub i ~-~ threshold ~~~~~ I sub i ~>~ 0.0$ $x sub c ~=~ {sum I sub i ~ x sub i} over {sum I sub i}$ .fi .IP 3. The errors are estimated in the following way .nf $sigma sup 2 ~=~ {sum I sub i ~ x sub i sup 2} over {sum I sub i} ~-~ x sub c sup 2$ $err sub xc ~=~ sqrt {{sigma sup 2} ~/~ sum I sub i}$ .fi .NH 4 Gaussian Fit to the Marginals .PP The fit is performed in the following way. .IP 1. The marginal distributions in x and y are accumulated. .IP 2. Initial guesses for the parameters of the 1D Gaussians $I sub Gauss$, $x sub c$ and $I sub sky$ are derived from the marginal distributions themselves. The width $sigma$ is held constant. .IP 3. The best fit parameters and their errors are derived using non-linear least-squares techniques and the NLFIT package. .NH 4 Optimal Filtering of Marginals .PP The fit is performed is the following way. .IP 1. The marginal distributions in x and y are accumulated. .IP 2. The centroid of the observed distribution is computed by solving the following equation. .nf $PSI (x sub c ) ~=~ sum omega sub i (x sub c ) ~ PHI sub i ~=~ 0$ $omega sub i ~=~ {{partial phi sub i} ~/~ {partial x sub c}} over {phi sub i~+~b}$ .fi The assumptions are that the observed distribution $PHI sub i$ is correctly modeled by the profile function $phi sub i$. The usual choise of profile model for centering is a Gaussian. Howver in the interests of speed a triangle function has been substituted. This causes at most a 7% increase in the random centering error. .IP 3. The startup procedure is based on the following fact. .nf $PSI (x sub n )~>~0~~for~~ x sub n ~<~ x sub c$ $PSI (x sub n ) ~<~0~~for~~ x sub n ~>~ x sub c$ .fi The iteration is initialized by assuming that $x sub 1$ = $x sub c$ and computing $PSI (x sub 1 )$. The initial x value is incremented by $+- sigma sub Gauss$ depending on the sign of $PSI$. The search is repeated until $PSI (x sub {n-1})$ and $PSI (x sub n)$ have opposite signs. At this point the true center is bracketed by the estimated positions $(x sub n, ~ x sub {n-1})$ and we have a table of at least 2 values of $PSI (x sub n )$. .IP 4. The computation proceeds by interpolating in the table of values for the estimated position $x sub {n+1}$ where $PSI$ = 0. If the table contains only two values as may be the case for the inititial interpolation, linear interpolation is used. In all other cases a quadratic fit to the three most recent $PSI$ values is used. The computation is complete when two successive estimates differ by less than some tolerance typically 0.001 pixel. .IP 5. The errors are estimated as follows. .nf $sigma sup 2 ~=~ ( {int {( partial phi ~/~ partial x sub c )} sup 2 over { phi + b}} ) sup -1$ $err sub xc ~=~ sqrt {sigma sup 2}$ .fi .NH 4 Other Centering Methods .PP The code is constructed in such a way that new algorithms may be easily added at a future date, such as the more sophisticated techniques required for accurate astrometry. .EQ delim $$ .EN .NH 3 The FITSKY Task .NH 4 Sky Fitting Routines .PP The following are the main entry points in the sky fitting package. .nf \f(CW apsinit (ap, function, annulus, dannulus, fwhmpsf, noise) ier = apfitsky (ap, im, xpos, ypos, sd, gd) ier = aprefitsky (ap, sd, gd) value = apstat[ir] (ap, param) apstats (ap, param, str, maxch) apset[irs] (ap, param, value) apsfree (ap) \fR .fi .PP The following quantities can be examined or set with apset/apstat calls. .nf \f(CW data 1. fwhmpsf full width half maximum of the psf parameters 2. sigma standard deviation of sky pixels sky fitting 1. annulus inner radius of sky annulus in fwhmpsf parameters 2. dannulus outer radius of sky annulus in fwhmpsf 3. sfunction sky fitting algorithm 4. smaxiter maximum number of fitting iterations 5. kreject k-sigma rejection limit sigma 6. nreject maximum number of rejection cycles 7. rgrow region growing radius in fwhmpsf 8. khist half-wdth of histogram in sigma 9. binsize binsize of histogram in sigma 10.smooth Lucy smooth the histogram 11.skyfile name of text file containing sky values 12.skyvaluser supplied constant value for sky \fR .fi .PP The following computed quantities can only be examined with apstat calls. .nf \f(CW 1. skymode computed sky value 2. skysig computed sky sigma 3. skyskew computed sky skew 4. nsky number of sky pixels 5. nsky_reject number of rejected sky pixels \fR .fi .NH 4 The Sky Fitting Algorithms .PP A good background fit is essential to aperture photometry. Fitting the background is trivial in a sparse field, but difficult in a crowded field. In general the background region will contain contaminating objects which must be detected and excluded if a good fit is to be obtained. .PP The main algorithm used in APPHOT is the following. .IP 1. If the skyfitting switch is disabled either read the sky values from a text file or accept a user supplied constant for the sky. .IP 2. Perform the initial sky fit using one of the specified algorithms. The sky fitting algorithms fall into three general categories, those that use the actual sky pixel array itself, those that operate on a histogram of sky values and those that rely on user interaction. .IP 3. If the pixel rejection flags are set perform pixels rejection with optional region growing. .NH 4 Sky Pixel Array Techniques .NH 5 Median .IP 1. Sort the array of sky pixels. This is necessary to avoid quantization effects. .IP 2. Compute the median, and the standard deviation and skew with respect to the mean. .IP 3. If the k-sigma rejection limit is greater than zero and the maximum number of rejection cycles is greater than one, perform pixel rejection. Pixels greater than k-sigma from the median are rejected. Region growing is optional. .IP 4. Stop the rejection cycle on any given iteration if the maximum number of rejection cycles ix exceeded, no more sky pixels are left or no more pixels are rejected. .NH 5 Mode .IP 1. Sort the array of sky pixels. This is necessary to avoid quantization effects. .IP 2. Compute the mode, and the standard deviation and skew with respect to the mean. .nf $I sub mode ~=~ 3.0 ~*~ I sub median ~-~ 2.0 ~*~ I sub mean$ .fi .IP 3. If the k-sigma rejection limit is greater than zero and the maximum number of rejection cycles is greater than one, perform pixel rejection. Pixels greater than k-sigma from the mode are rejected. Region growing is optional. .IP 4. Stop the rejection cycle on any given iteration if the maximum number of rejection cycles ix exceeded, no more sky pixels are left or no more pixels are rejected. .NH 4 Histogram Techniques .PP The following three techniques all operate on the histogram of the sky pixels. The routines all construct the histogram in the following identical manner. .IP 1. The mean of the sky distribution is computed. .IP 2. If the user specified standard deviation of the sky pixels is INDEF the algorithm computes the standard deviation of the sky pixels with respect to the mean. .IP 3. All pixels within plus or minus sigma standard deviations are accumulated into a histogram. The user specifies the bin size. .IP 4. The histogram may optionally be Lucy smoothed before any operation is performed on it. .NH 5 Centroid .PP The mode, sigma and skew of the sky pixels are computed in the following manner. .IP 1. The histogram is compiled as above. .IP 2. The mode, standard deviation and skew of the sky pixels are computed in the following manner. .nf $I sub 0 ~=~ sum I sub i$ $I sub 1 ~=~ sum I sub i ~ x sub i$ $I sub 2 ~=~ sum I sub i ~ x sub i sup 2$ $I sub 3 ~=~ sum I sub i ~ x sub i sup 3$ $I sub mode ~=~ {I sub 1} ~/~ {I sub 0}$ $sigma ~=~ {( I sub 2 ~/~ I sub 0 ~-~ I sub mode sup 2 )} sup {1/2}$ $skew ~=~ ( {I sub 3 ~/~ I sub 0 ~-~ I sub mode ~*~ sigma sup 2 ~-~ I sub mode sup 3} ) sup {1/3}$ .fi .IP 3 If pixel rejection is enabled sky pixels within a user supplied limit of the mode are rejected with optional region growing. .NH 5 Gaussian Fit .PP The mode, standard deviation and skew of the sky pixels are derived from a model fit in the following way. .IP 1. The histogram of the sky pixels is compiled as above. .IP 2. Initial guesses to the model parameters, $N sub max$, $I sub mode$, $sigma$, and $skew$ are made from the histogram itself. .IP 3. Final parameters and their errors are derived using non-linear least squares techniques and the NLFIT package. .IP 4. If pixel rejection is enabled sky pixels within a user supplied limit of the computed mode are rejected with optional region growing. .NH 5 Optimal Filtering .PP The method is as follows. .IP 1. The histogram is compiled as above. .IP 2. Using the mean of the sky pixels as the intital value of the sky mode, a new mode is computed using the optimal filtering technique described for centering. .IP 4 If pixel rejection is enabled sky pixels within a user supplied limit of the computed mode are rejected with optional region growing. .NH 5 Cross Correlation .PP The method is as follows. .IP 1. The histogram is compiled as above. .IP 2. The noise function is estimated using the standard deviation of the sky pixels and the cross-correlation function is computed. .IP 3. The mode is computed using quadratic interpolation around the peak of the distribution. .IP 4 If pixel rejection is enabled sky pixels within a user supplied limit of the mode are rejected with optional region growing. .NH 4 Interactive Techniques .NH 5 Histogram Plot .PP The histogram is compiled as described above and the user marks the peak on the histogram plot with the graphics cursor. The sigma and skew of the sky distribution with respect to the mean is also computed. .NH 5 Radial Distribution .PP A radial profile plot of the sky region is plotted and the user marks the sky value on the plot with the graphics cursor. The sigma and skew of the sky distribution with respect to the mean is computed. .NH 4 Pixel Rejection and Region Growing .PP All the sky fitting algorithms permit pixel rejection and optional region growing. Pixel rejection and region growing are performed by locating all pixels more than k * sigma from the mode, and blindly rejecting all pixels within a certain radius of each deviant pixel. This simple algorithm works well because the sample is large, and therefore there is little penalty for discarding pixels that might not be deviant. Region growing also tends to accelerate convergence significantly. .PP Very faint contaminating objects are difficult to detect and reject. If there are enough such objects, they should not be rejected, because there are probably a few in the object aperture as well. A higher sky sigma will be calculated and the computed uncertainty in the magnitude will increase. The best solution to this problem may be to increase the size of the annulus to minimize the bias. .NH 4 The Principal PHOT Routines .PP The main entries in the photometry routine are the following. .nf \f(CW apinit (ap, cfunction, cbox, sfunction, annulus, dannulus, aperts, napert, fwhmpsf, noise) ier = apfitcenter (ap, im, xinit, yinit) ier = aprefitcenter (ap) ier = apfitsky (ap, im, xcenter, ycenter, sd, gd) ier = aprefitsky (ap, sd, gd) ier = apmag (ap, im, xcenter, ycenter, skyval, skysig, nsky) ier = apwmag (ap, im, xcenter, ycenter, positive, skyval, skysig, nsky) ier = apremag (ap, positive, skyval, skysig, nsky) ier = apwremag (ap, positive, skyval, skysig, nsky) value = apstat[ir] (ap, param) apstats (ap, param, str, maxch) apset[sir] (ap, param, value) apfree (ap) \fR .fi .PP The following parameters can be examined or altered by apset/apstat calls. .nf \f(CW 1. weighting weighting scheme for wphot 2. aperts list of apertues 3. naperts number of apertures 4. zmag zero point of magnitude scale 5. itime effective integration time \fR .fi .PP The following quantities can be examined with apstat calls. .nf \f(CW 1. sums array of aperture sums 2. areas array of areas 3. mags array of magnitudes 4. magerrs array of magnitude errors \fR .fi .NH 4 The PHOT Aperture Integration Algorithm .PP The integral of the flux within a circular aperture is computed by fractional pixel techniques. Pixels are assumed to be square apertures arranged in a rectangular grid. The fraction of a pixel which lies within the circular APPHOT aperture is computed by an approximation, and all such contributions are summed to produce the total integral. .PP The inclusion of a partial pixel inside the aperture is done as follows. .IP 1. If the distance of the current pixel from the center of the star, r, is exactly equal to the radius of the aperture R then one-half the counts in the pixel are included. .IP 2. If r < R - 0.5 the entire pixel is included while if r > R + 0.5 the pixel is wholly excluded. .IP 3. In between the fraction of the counts varies linearly. A circular aperture is approximated by an irregular polygon. .PP The simplicity of aperture photometry limits the amount of information available for error analysis. The following three sources of error are considered. .IP 1. The error due to sky noise in the aperture. .nf $error sub 1 ~=~ sigma sub sky ~*~ {A sub apert} sup {1/2}$ .fi .IP 2. The error in the aperture sum. .nf $error sub 2 ~=~ ( {A sub "sum" ~/~ phpadu} ) sup {1/2}$ .fi .IP 3. The mean error of the sky. .nf $error sub 3 ~=~ sigma sub sky ~*~ A sub apert ~/~ nsky sup {1/2}$ .fi where $sigma sub sky$ is either computed by the background fitting algorithm or set by the user, and $A sub apert$ is the fractional pixel area of the aperture. .NH 4 The WPHOT Algorithm .PP The WPHOT algorithm computes a weighted aperture sum in an attempt to minimize noise in the sky. The algorithm is the following where w is the weight for each pixel, p is the noise free profile value and $sigma$ is the noise per pixel from all sources. (See the paper by Stover and Allen 1987 for details) .nf $A sub sum ~=~ sum {w sub i ~*~ (I sub i ~-~ sky)}$ $w sub i ~=~ C ~*~ p sub i ~/~ sigma sup 2 sub i$ $C ~=~ sum {p sub j} / sum {p sub j ~*~ w sub j}$ .fi .NH 4 The POLYPHOT ROUTINES .PP The principal polyphot routines are the following. .nf .na \f(CW apyinit (ap, sfunction, annulus, dannulus, noise) ier = apfitcenter (ap, im, wx, wy) ier = aprefitcenter (ap) ier = apfitsky (ap, im, xcenter, ycenter, sd, gd) ier = aprefitsky (ap, sd, gd) ier = polyfit (ap, im, xver, yver, nver) value = apstat[ir] (ap, param) apstats (ap, param, str, maxch) apset[sir] (ap, param, value) apfree (ap) \fR .ad .fi .PP .NH 4 The POLYPHOT Algorithm .PP The function of the POLYPHOT task is to compute the flux inside an irregular polygon given a list of the coordinates of the vertices of a polygon. The polygon must be entirely inside the image and the vertices of the polygon must be specified in clockwise or counterclockwise order. The actual algorithm used is as follows. .IP 1. The range of image lines which intersect the polygon are computed. .IP 2. For each image line in the specified range the intersection points with the polygon are computed. .IP 3. The flux between pairs of limits is summed using a fractional pixel approximation for the endpoints. .IP 4. The sky is fitted using any of the methods previously discussed and a user specified annular region. .IP 5. The errors are computed as specified in the PHOT specifications. .PP .EQ delim $$ .EN .NH Example .PP A brief example may help illustrate the use of the package. Suppose we want to process a few hundred stars on image "blue". .PP The first step is to prepare a list of objects to be measured. The simplest way to do this is to interactively mark the objects with the image cursor using the display (graphics) device and the RIMCURSOR (RGCURSOR) task. .nf \f(CW ... load image on the display ... ap> rimcursor > starlist ... move cursor and mark stars ... ... load contour plot on graphics terminal ... ap> rgcursor > starlist ... move cursor and mark stars ... \fR .fi Alternatively one can run DAOFIND to compute a list of candidate objects in the frame. The name of the coordinate file is stored in the PHOT parameter set. .nf \f(CW ap> phot.coords=starlist \fR .fi .PP The next step is to set up the PHOT parameters interactively. First we load the image (contour plot) blue on the display (graphics terminal). Next we call up PHOT in interactive mode. .nf \f(CW ap> phot blue ... cursor appears ... \fR .fi .PP PHOT takes input by reading the image (graphics) display (terminal) cursor. In order to display the available commands we tap the ? key and the following text appears on the screen. .nf \f(CW Interactive Phot Commands ? Print options : Colon command see below i Setup PHOT parameters interactively w Write PHOT parameters to the parameter files l Process the remainder of the coordinate list r Rewind the coordinate list c Fit center around the current cursor position t Fit sky around the current cursor position s Fit sky around the current center position p Compute magnitudes around the cursor position f Fit center, sky and compute magnitudes sp Fit center, sky, compute magnitudes, and save q Exit program Phot parameters are listed or set with the following commands. :m [n] Move cursor to the [nth] object in the coordinate list :n [n] Measure the [nth] object in the coordinate list :show [center/sky/phot/all] List the aphot parameters :fwhmpsf [value] Full width half maximum of the PSF :noise [string] Noise model :threshold [value] Threshold value for centering :sigma [value] Standard deviation of the background :ccdread CCD readout noise keyword :readnoise Readout noise in electrons :gain Gain keyword :epadu Electrons per adu :calgorithm [string] Centering function :positive [y/n] Emission or absorption feature :cbox [value] Width of the centering box in fwhmpsf :cmaxiter [value] Maximum number of centering iterations :maxshift [value] Maximum shift in fwhmpf :minsnratio [value] Minimum signal to noise ratio of pixels :clean [y/n] Clean subraster before centering :rclip [value] Clipping radius in fwhmpsf :rclean [value] Cleaning radius in fwhmpsf :kclean [value] Sigma for clean algorithm :mkcenter [y/n] Mark the centers on the display :salgorithm [string] Sky fitting algorithm :annulus [value] Inner radius of sky annulus in fwhmpsf :dannulus [value] Width of sky annulus in fwhmpsf :skyvalue [value] User supplied sky :smaxiter [value] Maximum number of rejection cycles :skreject [value] +/- Pixel rejection limits in sky sigma :snreject [value] Maximum number of rejection interations :khist [value] +/- Sky histogram size in sky sigma :binsize [value] Resolution of sky histogram in sky sigma :smooth [y/n] Lucy smooth the sky histogram :rgrow [value] Region growing radius in fwhmpsf :marksky [y/n] Mark the sky annuli on the display :weighting Weighting for wphot :aperts [string] Aperture radii in fwhmpsf :zmag [value] Zero point of magnitude scale :exposure [string] Exposure time keyword :itime [value] Integration time \fR .fi .PP We select the interactive setup option, move the image cursor to a high signal-to-noise, isolated star and tap the i key. PHOT responds by plotting the radial profile of the star on the screen and requesting the user to mark the fwhm of the psf, the centering aperture, the inner and outer sky annuli, the sky background and sigma and the set of circular apertures. The parameters so set can be examined and/or reset with the : commands as shown above. Sample measurements can be made of several stars by moving the cursor and typing the f command. Finally when we are happy with the parameter set we type w to store the parameters and q to exit the program. .PP Now we are ready to do photometry. We enter the PHOT program in batch mode. .nf \f(CW ap> phot blue inter- & \fR .fi The batch job is now running, appending output lines to the file "blue.mag.#". We can proceed to set up the job for the red image, in much the same way that we set up the job for the blue image. When both jobs finish, we can use the list processing tools to filter out the good objects and calculate colors. .NH The APPHOT Tasks .PP Manual pages for the APPHOT tasks are available in the IRAF on line help database.