Initial commit

1 files changed, 1296 insertions, 0 deletions

diff --git a/noao/digiphot/apphot/doc/specs/apphot.spc b/noao/digiphot/apphot/doc/specs/apphot.spc
new file mode 100644
index 00000000..ecfdf364
--- /dev/null
+++ b/noao/digiphot/apphot/doc/specs/apphot.spc
@@ -0,0 +1,1296 @@
+.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:\fL
+
+.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\fR
+       wphot	-- perform weighted multi-aperture photometry
+.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> ~=~ I sub {x ~ Gauss} ~+~ I sub {y ~ Gauss}$
+
+	$roundlo ~<=~  round ~=~ {DELTA I} over <I> ~<=~  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
+\fL
+
+                  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
+\fL
+  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
+\fL
+              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
+\fL
+         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
+\fL
+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
+\fL
+              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
+\fL
+               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
+\fL
+                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
+\fL
+
+                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
+\fL
+               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
+.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
+\fL
+
+	... 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
+\fL
+	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
+\fL
+	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
+\fL
+	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
+\fL
+	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 APHOT tasks are attached. Although the working of the
+various tasks may change in detail, the help pages attached here should
+provide a good description of the function of each task.
+.PP
+When the package has stabilized a detailed users guide will be written.