Context Navigation

-              r38591
+              r39920
 % \pdfoutput=1
+% see latex.readme.txt for notes on using the PS1 template
+%\documentclass[12pt,preprint]{aastex}
+%\documentclass[manuscript]{aastex}
+%\documentclass[preprint2]{aastex}
+%\documentclass[preprint2,longabstract]{aastex}
+\RequirePackage{graphicx}
 \RequirePackage{color}
+\RequirePackage{code}
 \input{astro.sty}
 …
 %\def\plotmode{bw}
+% latex uses the .ps regardless of \plotext if passed with ext=.\plotext
+% pdflatex uses the png if it exists, pdf otherwise OK with pdf, png
+% independent of graphicx or graphics
+% pdflatex fails on files with multiple dots, protect with {file.name.long}.pdf
+% use ps2pdf -dEPSCrop input.ps to convert to pdf
+% \includegraphics[width=\hsize,clip]{peaks}
+% \includegraphics[width=\hsize]{temp.pdf}
+% \includegraphics[width=\hsize]{{FWHM.smooth.trend.ps1}.pdf}
 %\def\plotext{pdf}
+\def\plotext{ps}
+%\def\plotext{png}
+\def\plotext{pdf}
 %\def\picdir{/home/eugene/chipresid.20140404}
 \def\picdir{/data/pikake.2/eugene/chipresid.20140404}
+\def\picdir{pics}
 % Pick a terse version of the title here;
 …
 \shortauthors{E.A. Magnier et al}
 \begin{document}
 \title{Pan-STARRS Pixel Analysis : Source Detection \& Characterization}
+\title{Pan-STARRS Pixel Analysis : Source Detection and Characterization}
 % this is a crude trick to get the order of affiliations right.  These
 …
 % list and (2) re-order the list at the bottom (and comment-out as needed)
 \def\IfA{1}
+\def\Princeton{2}
+\def\DUR{3}
 \def\CfA{2}
-\def\MPIA{3}
-\def\Princeton{3}
-\def\USNO{4}
-\def\JHU{1}
 % This example has a first author from UH:
 \author{
 Eugene A. Magnier,\altaffilmark{\IfA}
+IPP Team,
+%PS Builder List
+% R. H. Lupton,\altaffilmark{\Princeton}
+W.~E. Sweeney,\altaffilmark{\IfA}
+K.~C. Chambers,\altaffilmark{\IfA}
+H.~A. Flewelling,\altaffilmark{\IfA}
+M. E. Huber,\altaffilmark{\IfA}
+P.~A. Price,\altaffilmark{\Princeton}
+C. Z. Waters,\altaffilmark{\IfA}
+% PS1 Builders
+L. Denneau,\altaffilmark{\IfA}
+P. Draper,\altaffilmark{\DUR}
+R. Jedicke,\altaffilmark{\IfA}
+K. W. Hodapp,\altaffilmark{\IfA}
+N. Kaiser,\altaffilmark{\IfA}
+R.-P. Kudritzki,\altaffilmark{\IfA}
+N. Metcalfe,\altaffilmark{\DUR}
+C.~W. Stubbs,\altaffilmark{\CfA}
 % W.~S. Burgett,\altaffilmark{\IfA}
-% K.~C. Chambers,\altaffilmark{\IfA}
-% L. Denneau,\altaffilmark{\IfA}
-% P. Draper,\altaffilmark{\DUR}
-% H.~A. Flewelling,\altaffilmark{\IfA}
 % T. Grav,\altaffilmark{\IfA}
 % J. N. Heasley,\altaffilmark{\IfA}
-% K. W. Hodapp,\altaffilmark{\IfA}
-% M. E. Huber,\altaffilmark{\IfA}
-% R. Jedicke,\altaffilmark{\IfA}
-% N. Kaiser,\altaffilmark{\IfA}
-% R.-P. Kudritzki,\altaffilmark{\IfA}
 % G. A. Luppino,\altaffilmark{\IfA}
 % R. H. Lupton,\altaffilmark{\Princeton}
 % E. A. Magnier,\altaffilmark{\IfA}
-% N. Metcalfe,\altaffilmark{\DUH}
 % D. G. Monet,\altaffilmark{\USNO}
 % J.~S. Morgan,\altaffilmark{\IfA}
 % P. M. Onaka,\altaffilmark{\IfA}
-% P.~A. Price,\altaffilmark{\Princeton}
-% C.~W. Stubbs,\altaffilmark{\CfA}
-% W.~E. Sweeney,\altaffilmark{\IfA}
 % J.~L. Tonry, \altaffilmark{\IfA}
+% R. J. Wainscoat,\altaffilmark{\IfA} and
+% C. Z. Waters,\altaffilmark{\IfA}
+R. J. Wainscoat\altaffilmark{\IfA}
 } % this bracket terminates author list
 % The ordering here should be sequential, matching the sequence in the list of authors:
 \altaffiltext{\IfA}{Institute for Astronomy, University of Hawaii, 2680 Woodlawn Drive, Honolulu HI 96822}
+% \altaffiltext{\CfA}{Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138}
+% \altaffiltext{\Princeton}{Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA}
+\altaffiltext{\Princeton}{Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA}
+\altaffiltext{\DUR}{Department of Physics, Durham University, South Road, Durham DH1 3LE, UK}
+\altaffiltext{\CfA}{Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138}
 % \altaffiltext{\USNO}{US Naval Observatory, Flagstaff Station, Flagstaff, AZ 86001, USA}
 % \altaffiltext{\JHU}{Department of Physics and Astronomy, Johns Hopkins University, 3400 North Charles Street, Baltimore, MD 21218, USA}
 …
 \begin{abstract}
+Lorem ipsum dolor sit amet, consectetur adipiscing elit. Vestibulum
+bibendum nisi id tristique posuere. Duis eu mollis nulla. Maecenas est
+turpis, mattis tempor urna vitae, placerat rhoncus sem. Lorem ipsum
+dolor sit amet, consectetur adipiscing elit. Sed quis velit
+nisl. Aliquam erat volutpat. Cras lacinia, nisl tristique auctor
+molestie, dolor nulla rhoncus purus, ac accumsan nunc nunc ac
+nibh. Maecenas vitae mollis mauris. Ut sollicitudin pulvinar purus,
+eget luctus lorem tincidunt vitae. Vestibulum eu mattis neque. Nulla
+in tortor id urna dapibus gravida a vel leo.
+Over 3 billion astronomical objects have been detected in the more
+than 22 million orthogonal transfer CCD images obtained as part of the
+Pan-STARRS\,1 $3\pi$ survey.  Over 85 billion instances of those
+objects have been automatically detected and characterized by the
+Pan-STARRS Image Processing Pipeline photometry software,
+\code{psphot}.  This fast, automatic, and reliable software was
+developed for the Pan-STARRS project, but is easily adaptable to
+images from other telescopes.  We describe the analysis of the
+astronomical objects by \code{psphot} in general as well as for the
+specific case of the 3rd processing version used for the first public
+release of the Pan-STARRS $3\pi$ survey data.
 \end{abstract}
 …
 \keywords{Surveys:\PSONE }
-\section{OUTLINE}
-\begin{verbatim}
-Intro
- Pan-STARRS background
- Scope: Source Detection \& Characterization, Galaxy modeling
- Requirements / Goals
- Comparable programs
- PSPhot
-\end{verbatim}
 \section{INTRODUCTION}\label{sec:intro}
+\note{more PS1 background}
+The Pan-STARRS Image Processing Pipeline is responsible for the basic
+analysis of images from the Pan-STARRS telescopes Gigapixel Camera.
+The photometric and astrometric precision goals for the all-sky
+surveys, as well as the other survey components, are quite stringent:
+\begin{itemize}
+\item relative photometry: 10 millimagnitudes scatter for bright stars
+across the sky in the internal photometric system;
+\item relative astrometry; 10 milliarcseconds scatter for individual
+stars between repeated images.
+\item absolute astrometry: 100 milliarcseconds scatter for all ICRS
+  reference stars (Tycho).
+\end{itemize}
+An additional constraint on the Pan-STARRS system comes from the high
+data rate.  PS1 produces typically $\sim 700$ GB per night of imaging
+data.  These images range from high galactic latitudes to the Galactic
+Bulge, so large numbers of measurable stars can be expected in much of
+the data.  The combination of the high precision goals of the
+astrometric and photometric measurements and the high data rate (and a
+finite computing budget) mean that the process of detecting,
+classifying, and measuring the astronomical objects in the image data
+stream in a timely fashion are a significant challenge.
+This is the fourth in a series of seven papers describing the
+Pan-STARRS1 Surveys, the data reduction techiques and the resulting
+data products.  This paper (Paper IV) describes the details of the
+source detection and photometry, including point-spread-function and
+extended source fitting models, and the techniques for ``forced''
+photometry measurements.
+%Chambers et al. 2017 (Paper I)
+%The Pan-STARRS\,1 Surveys
+\citet[][Paper I]{chambers2017}
+provides an overview of the Pan-STARRS System, the design and
+execution of the Surveys, the resulting image and catalog data
+products, a discussion of the overall data quality and basic
+characteristics, and a brief summary of important results.
+%Magnier et al. 2017 (Paper II)
+%Pan-STARRS Data Processing Stages
+\citet[][Paper II]{magnier2017c}
+describes how the various data processing stages are organised and implemented
+in the Imaging Processing Pipeline (IPP), including details of the
+the processing database which is a critical element in the IPP infrastructure .
+%Waters et al. 2017 (Paper III)
+%Pan-STARRS Pixel Processing : Detrending, Warping, Stacking
+\citet[][Paper III]{waters2017}
+describes the details of the pixel processing algorithms, including detrending, warping, and adding (to create stacked images) and subtracting (to create difference images) and resulting image products and their properties.
+%Magnier et al. 2017 (Paper IV)
+%Pan-STARRS Pixel Analysis : Source Detection
+%\citet[][Paper IV]{magnier2017a}
+%describes the details of the source detection and photometry, including point-spread-function and extended source fitting models, and the techniques for ``forced" photometry measurements.
+%Magnier et al. 2017 (Paper V)
+%Pan-STARRS Photometric and Astrometric Calibration
+\citet[][Paper V]{magnier2017b}
+describes the final calibration process, and the resulting photometric and astrometric quality.
+%Flewelling et al. 2017 (Paper VI)
+%Pan-STARRS 1 Database and Data Products
+\citet[][Paper VI]{flewelling2017}
+describes  the details of the resulting catalog data and its organization in the Pan-STARRS database.
+%
+%
+\citet[][Paper VII]{huber2017}
+%Huber et al. 2017 (Paper VII)
+describes the Medium Deep Survey in detail, including the unique issues and data products specific to that survey. The Medium Deep Survey is not part of Data Release 1. (DR1)
+%
+The Pan-STARRS1 filters and photometric system have already been
+described in detail in \cite{2012ApJ...750...99T}.
+{\color{red} {\em Note: These papers are being placed on arXiv.org to
+    provide crucial support information at the time of the public
+    release of Data Release 1 (DR1). We expect the arXiv versions to
+    be updated prior to submission to the Astrophysical Journal in
+    January 2017. Feedback and suggestions for additional information
+    from early users of the data products are welcome during the
+    submission and refereeing process.}}
+\section{Background}
+The photometric and astrometric precision goals for the Pan-STARRS\,1
+surveys were quite stringent: photmetric accuracy of 10
+millimagnitudes, relative astrometric accuracy of 10 milliarcseconds
+and absolute astrometric accuracy of 100 milliarcseconds with respect
+to the ICRS reference stars.
+An additional constraint on the Pan-STARRS analysis system comes from
+the high data rate.  PS1 produces typically $\sim 500$ exposures per
+night, corresponding to $\sim 750$ billion pixels of imaging data.
+The images range from high galactic latitudes to the Galactic bulge,
+so large numbers of measurable stars can be expected in much of the
+data.  The combination of the high precision goals of the astrometric
+and photometric measurements and the high data rate (and a finite
+computing budget) mean that the process of detecting, classifying, and
+measuring the astronomical objects in the image data stream in a
+timely fashion are a significant challenge.
 In order to achieve these ambitious goals, the object detection,
 …
 astrometry.
-\subsection{Comparable Programs}
 A variety of astronomical software packages perform the basic object
 detection, measurement, and classification tasks needed by the
 …
   pro: well-tested, stable code.  con: limited range of models,
   algorithm converges slowly to a PSF model, limited tests of PSF
   validity, inflexible code base, fortran (P. Schechter)
+  validity, inflexible code base, fortran \citep{1993PASP..105.1342S}.
 \item DAOPhot : Pixel-map PSF model with analytical component.  pro:
   well-tested, high-quality photometry.  con: Difficult to use in an
+  automated fashion, does it handle 2D variations well? (P. Stetson)
+\item Sextractor : pure aperture measurement with rudimentary
+  object subtraction.  pro: fast, widely used, easy to automate.  con:
+  poor object separation in crowded regions, PSF-modeling is only
+  beta (psfex), what models are available? (E. Bertin)
+\item apphot : IRAF-based aperture photometry.  pro: widely used.
+  con: IRAF-based, aperture photometry. (???)
+  automated fashion, does it handle 2D variations well? \citep{1987PASP...99..191S}.
+\item Sextractor : pure aperture measurement with rudimentary object
+  subtraction.  pro: fast, widely used, easy to automate.  con: poor
+  object separation in crowded regions, PSF-modeling was only in beta,
+  not widely used at the time \citep{sextractor}.
 \item galfit : detailed galaxy modeling.  not a multi-object PSF
   analysis tool.  con: does not provide a PSF model, not easily
   automated.  very detailed results in very slow processing.  only a
   galaxy analysis program. (C. Impey)
+  galaxy analysis program \citep{2002AJ....124..266P}.
 \item SDSS phot : con: tightly integrated into the SDSS software
   environment.  (R. Lupton)
+  environment \citep{2001ASPC..238..269L}.
 \end{itemize}
+\note{re-phrase this:} The Pan-STARRS IPP team decided that none of
+the existing packages met all of their needs, particularly given the
+very challenging goals of the project.  We decided to redesign the
+photometry analysis from scratch, using the lessons learned from the
+existing photometry systems.  In the process, the object analysis
+software would be written using the data analysis C-code library
+written for the IPP, \code{psLib}, and the components of the
+photometry code would be integrated into the IPP's mid-level astronomy
+data analysis toolkit called \code{psModules}.  The result is
+'PSPhot', which can be used either as a stand-alone C program, or as
+callable set of functions.
+\note{discuss the psphot program varients}
+\begin{verbatim}
+Other Varients:
+* psphotStack -- 5 filter simultaneous fitting
+* psphotFullForce
+\end{verbatim}
+When the IPP development was starting, the existing photometry
+packages either did not meet the level of accuracy required or were
+required too much human intervention to be considered for the needs of
+PS1.  In the case of the SDSS Photo tool, the software was judged to
+be too tightly integrated to the architecture of SDSS to be easily
+re-integrated into the Pan-STARRS pipelin.  A new photometry analysis
+package was developed using lessons learned from the existing
+photometry systems.  In the process, the object analysis software was
+written using the data analysis C-code library written for the IPP,
+\code{psLib}.  Components of the photometry code were integrated into
+the IPP's mid-level astronomy data analysis toolkit called
+\code{psModules}.  The result software, 'PSPhot', can be used either
+as a stand-alone C program, or as a set of library functions which may
+be integrated into other programs
+The main version of PSPhot is a stand-alone program which is run on a
+single image or a group of related images representing the data read
+from a camera in a single exposure.  The images are expected to have
+already been detrended so that pixel values are linearly related to
+the flux.  The gain may be specific by the configuration system, or a
+variance image may be supplied.  A mask may also be supplied to mark
+good, bad, and suspect pixels.  Several variants of psphot have also
+been used in the PS1 PV3 analysis.
+The version called PSPhotStack accepts a set of images, each
+representing the same patch of sky in a different filter, nominally
+the full $grizy$ filter set for the analysis of the PS1 PV3 stack
+images, though where insufficient data were available in a given
+filter, a subset of these filters was processed as a group.  As
+discussed in detail below, the PSPhotStack analysis includes the
+capability of measuring forced PSF photometry in some filter images
+based on the position of sources detected in the other filters.  It
+also include an option to convolve the set of images to a single,
+common PSF size across the filters for the purpose of fixed aperture
+photometry.
+A second version of PSPhot used in the PV3 analysis is called
+PSPhotFullForce.  In this version, a set of image all representing the
+same pixels are processed together, with the positions of sources to
+be analysed loaded from a supplied file.  In this version the
+analysis, sources are not discovered -- only the supplied sources are
+considered.  PSF models are determined for each exposure and the
+forced PSF photometry is measured for all sources.  A subset of
+sources may also be used to measure forced galaxy shape parameters.
+As described below, a grid of galaxy models are fitted based on the
+supplied guess model.
 \section{PSPhot Design Goals}
 …
   level must be reached for images with 250 mas pixels, implying
   PSPhot must introduce measurement errors less than 1/50th of a
+  pixel. \note{the choice of F32 parameters places a numerical limit
+  of 1e-7 on the accuracy of a pixel relative to the size of a chip
+  (since a single data value is used for X or Y).  For the $4800^2$
+  GPC chips, this yields a limit of about 0.25 milliarcsecond.}
+  pixel. The choice of 32 bit floating point data values for the
+  source centroids places a numerical limit of 1e-7 on the accuracy of
+  a pixel relative to the size of a chip (since a single data value is
+  used for X or Y).  For the $4800^2$ GPC chips, this yields a limit
+  of about 0.25 milliarcsecond.
 \end{itemize}
 …
 \item {\bf Flexible non-PSF models} PSPhot must be able to represent
   PSF-like objects as well as non-PSF sources.  It must be easy to add
   new object models as interesting representations of sources are
   invented.
+  PSF-like objects as well as non-PSF sources (e.g., galaxies).  It
+  must be easy to add new object models as interesting representations
+  of sources are invented.
 \item {\bf Clean code base} PSPhot should incorporate a high-degree of
 …
   provide the user with methods for assessing the different PSF models.
 \item {\bf Careful aperture corrections} PSPhot must carefully measure
   and correct for the photometric and astrometric trends introduced by
   using analytical PSF models.
+\item {\bf Careful systematic corrections} PSPhot must carefully
+  measure and correct for the photometric and astrometric trends
+  introduced by using analytical PSF models.
 \item {\bf User Configurable} PSPhot should allow users to change the
 …
 The PSPhot analysis is divided into several major stages:
 \begin{itemize}
+\begin{enumerate}
 \item {\bf Image preparation} Load data, characterize the image
   background, load or construct variance and mask images.
 \item {\bf Initial object detection} Smooth, find peaks, measure basic
   properties
+  properties.
 \item {\bf PSF determination} Select PSF candidates, perform model
 …
   properties (aperture or PSF)
+\item {\bf Extended Source Analysis} Detailed measurements relevant to
+  galaxies and/or other extended (non-PSF) sources.
 \item {\bf Aperture corrections} Measure the curve-of-growth, spatial
   aperture variations, and background-error corrections.
 …
 \item {\bf Output} Write out objects in selected format, write out
   difference image, variance image, etc, as selected.
 \end{itemize}
+\end{enumerate}
 PSPhot is highly configurable.  Users may choose via the configuration
 system which of the above analyses are performed.  This may be useful
+for testing, but may also allow for specialized use cases.  For
+example, the PSF model may already be available from external
 information, in which case the PSF modeling stage can be skipped.
+system which of the above analyses are performed.  This is useful for
+testing, but also allows for specialized use cases.  For example, the
+PSF model may already be available from external information, in which
+case the PSF modeling stage can be skipped.
 \subsection{Image Preparation}
 …
 circumstance, while a pixel in which persistence ghosts have been
 subtracted might be useful for detection or even analysis of brighter
+sources.  \note{can I identify which functions respect which sets of masks}
+sources.  Table~\ref{tab:mask_values} lists the 16 bit values used for
+PS1 mask images, along with their description \citep[see][for
+  additional information]{waters2017}.
+\begin{table*}
+\caption{\label{tab:mask_values} PSPhot / GPC1 Mask Image Pixel Values}\vspace{-0.5cm}
+\begin{center}
+\begin{tabular}{lcl}
+\hline
+\hline
+{\bf Mask Name} & {\bf Mask Value} & {\bf Description} \\
+\hline
+  DETECTOR & 0x0001 & A detector defect is present. \\
+  FLAT     & 0x0002 & The flat field model does not calibrate the pixel reliably. \\
+  DARK     & 0x0004 & The dark model does not calibrate the pixel reliably. \\
+  BLANK    & 0x0008 & The pixel does not contain valid data. \\
+  CTE      & 0x0010 & The pixel has poor charge transfer efficiency. \\
+  SAT      & 0x0020 & The pixel is saturated. \\
+  LOW      & 0x0040 & The pixel has a lower value than expected. \\
+  SUSPECT  & 0x0080 & The pixel is suspected of being bad. \\
+  BURNTOOL & 0x0080 & The pixel contain an burntool repaired streak. \\
+  CR       & 0x0100 & A cosmic ray is present. \\
+  SPIKE    & 0x0200 & A diffraction spike is present. \\
+  GHOST    & 0x0400 & An optical ghost is present. \\
+  STREAK   & 0x0800 & A streak is present. \\
+  STARCORE & 0x1000 & A bright star core is present. \\
+  CONV.BAD & 0x2000 & The pixel is bad after convolution with a bad pixel. \\
+  CONV.POOR& 0x4000 & The pixel is poor after convolution with a bad pixel. \\
+  MARK     & 0x8000 & An internal flag for temporarily marking a pixel. \\
+\hline
+\end{tabular}
+\end{center}
+\end{table*}
 The variance image, if not supplied is constructed by default from the
 flux image using the configuration supplied values of \code{GAIN} and
 \code{READ\_NOISE} to calculate the appropriate Poisson statistics for
+\code{READ_NOISE} to calculate the appropriate Poisson statistics for
 each pixel.  In this case, the image is assumed to represent the
 readout from a single detector, with well-defined gain and read noise
 …
 Since a typical smoothing or warping operation may introduce
 correlation between 25 - 100 neighboring pixels, the size of such a
+covariance image is prohibitive.  In practice, however, there are two
+extreme cases which generally are relevant.  \note{talk about the
+  covar matrix for a PSF}
+\subsection{Background (Sky) Model}
+covariance image is prohibitive.
+Before sources are detected in the image, a model of the background is
+subtracted.  The image is divided into a grid of background points
+with a spacing of 400 pixels.  Superpixels of size $800\times 800$
+pixels are used to measure the local background for each background
+grid point, thus over-sampling the background spatial variations by a
+factor of 2.  In the interest of speed, 10,000 randomly selected {\em
+  unmasked} pixels in these regions are sampled to determine the
+background.  Bilinear interpolation is used to generate a
+full-resolution image from the grid of background points, and this
+image is then subtracted from the science image.  The background image
+and the background standard deviation image are kept in memory from
+which the values of \code{SKY} and \code{SKY_SIGMA} are calculated for
+each object in the output catalog.  See also the discussion in
+\cite{waters2017}.
 \subsection{Initial Object Detection}
+\subsubsection{Peak Detection}
+\label{sec:peaks}
 The objects are initially detected by finding the location of local
 peaks in the image.  The flux and variance images are smoothed with a
+small circularly symmetric kernel using a two-pass 1D Gaussian
+(\note{KEYWORD?}).  The smoothed flux and variance images are combined
+to generate a significance image in signal-to-noise units
+\note{including correction for the covariance, if known}. At this
+stage, the goal is only to detect the brighter sources, above a user
+defined S/N limit (configuration keyword: \code{PEAK\_NSIGMA}).  The
+detection efficiency for the brighter sources is not strongly
+dependent on the form of this smoothing function.
+small circularly symmetric kernel using a two-pass 1D Gaussian.  The
+smoothed flux and variance images are combined to generate a
+significance image in signal-to-noise units, including correction for
+the covariance, if known. At this stage, the goal is only to detect
+the brighter sources, above a user defined S/N limit (configuration
+keyword: \code{PEAKS_NSIGMA_LIMIT}).  A maximum of
+\code{PEAKS_NMAX} are found at this stage.  The detection efficiency
+for the brighter sources is not strongly dependent on the form of this
+smoothing function.
 The local peaks in the smoothed image are found by first detecting
 …
 the maximum $X$ and $Y$ corners of the region.
-\subsection{Footprints}
-\note{need to describe the process of generating the source footprints
-  and then culling the insignificant peaks}
-\subsubsection{Moments and related}
-\note{disucss the Kron mags}
-\note{this section is wrong: we no longer use S/N clipping, but a
-  Gaussian window function, chosen based on the measured moment}
-Once a collection of peaks have been identified, basic properties of
-the objects are measured.  First, the local sky flux is measured
-within a square annulus with user-defined dimensions
-(\code{INNER\_RADIUS} and \code{OUTER\_RADIUS}), using the sample
-median.  This local background value is then used to calculate the
-object first and second moments within a small user-defined aperture
-(\code{MOMENT\_RADIUS}).  The first-order moments are a good
-representation of the object position, while the second-order moments
-are a measure of the object shape.  The second-order moments are
-somewhat sensitive to the size of the aperture and the accuracy of the
-background measurement.  The moment calculation is only performed
-using pixels which exceed a S/N of 1.  If, in the process of
-calculating the source moments, the S/N limits reject all but \note{3}
-or fewer of the source pixels, the peak is identified as being
-suspect, and is not used for further analysis.  If the measured
-centroid coordinates differ from the peak coordinates be a large
-amount (\code{MOMENT\_RADIUS}), then the peak is again identified as
-being of poor quality and is rejected.  In both of these cases, it is
-likely that the `peak' was identified in a region of flat flux
-distribution or many saturated or edge pixels.
-\subsubsection{Determination of the Peak Coordinates and Errors}
-\note{this section is wrong: it is a poor estimator of the source
-  position errors.  we gave up a reverted to using the FWHM / (S/N)}
 We use the 9 pixels which include the source peak to fit for the
 position and position errors.  We model the peak of the sources as a
 …
 of only 0 or 1, we can greatly simplify the chi-square equation to a
 square matrix equation with the following values:
+%% fix this:
+\begin{verbatim}
+| 9 0 0 0 6 6 | C_00 | = \sum F_{i,j}
+| 0 6 0 0 0 0 | C_10 | = \sum F_{i,j} x
+| 0 0 6 0 0 0 | C_01 | = \sum F_{i,j} y
+| 0 0 0 6 0 0 | C_11 | = \sum F_{i,j} x y
+| 6 0 0 0 6 4 | C_20 | = \sum F_{i,j} x^2
+| 6 0 0 0 4 6 | C_02 | = \sum F_{i,j} y^2
+\end{verbatim}
+The inverse of the 3x3 matrix terms for $C_{00}$, $C_{20}$, and $C_{02}$ is:
+\begin{verbatim}
+| +5/9 -1/3 -1/3 |
+| -1/3 +1/2    0 |
+| -1/3    0 +1/2 |
+\end{verbatim}
+which can be used to determine the errors on the coefficients:
+\[
+\left( \begin{array}{cccccc}
+& 0 & 0 & 0 & 6 & 6 \\
+& 6 & 0 & 0 & 0 & 0 \\
+& 0 & 6 & 0 & 0 & 0 \\
+& 0 & 0 & 6 & 0 & 0 \\
+& 0 & 0 & 0 & 6 & 4 \\
+& 0 & 0 & 0 & 4 & 6 \\
+\end{array} \right)
+\left( \begin{array}{c}
+C_{00}\\
+C_{10}\\
+C_{01}\\
+C_{11}\\
+C_{20}\\
+C_{02}\\
+\end{array} \right)
+=
+\left( \begin{array}{c}
+\sum F_{i,j}     \\
+\sum F_{i,j} x   \\
+\sum F_{i,j} y   \\
+\sum F_{i,j} x y \\
+\sum F_{i,j} x^2 \\
+\sum F_{i,j} y^2 \\
+\end{array} \right)
+\]
+Inverting the 3x3 matrix terms for $C_{00}$, $C_{20}$, and $C_{02}$,
+the location of the peak is determined from the minimum of the
+bi-quadratic function above, and is given by:
 \begin{eqnarray}
+\sigma^2_{00} & = & \sigma^2 (5/9) \\
+\sigma^2_{10} & = & \sigma^2 (1/6) \\
+\sigma^2_{01} & = & \sigma^2 (1/6) \\
+\sigma^2_{11} & = & \sigma^2 (1/6) \\
+\sigma^2_{20} & = & \sigma^2 (1/2) \\
+\sigma^2_{02} & = & \sigma^2 (1/2) \\
+x_{min} & = & (C_{11} C_{01} - 2 C_{02} C_{10}) D^{-1} \\
+y_{min} & = & (C_{11} C_{10} - 2 C_{20} C_{01}) D^{-1} \\
+D      & = & 4 C_{20} C_{02} - C_{11}^2
 \end{eqnarray}
+The location of the peak is determined from the minimum of the
+bi-quadratic function above, and is given by:
+\begin{figure}[htbp]
+  \begin{center}
+ \includegraphics[width=\hsize,clip]{pics/peaks}
+  \caption{\label{fig:peaks} Illustration of peak finding and culling peaks within a
+    footprint.  Insignificant peaks within the footprint of a brighter
+    peak are ignored in further processing. }
+  \end{center}
+\end{figure}
+\subsubsection{Footprints}
+The peaks detected in the image may correspond to real sources, but
+they may also correspond to noise fluctuations, especially in the
+wings of bright stars.  PSPhot attempts to identify peaks which may be
+formally significant, but are not locally significant.  It first
+generates a set of ``footprints'', contiguous collections of pixels in
+the smoothed significance image above the detection threshold.  These
+regions are grown by a small amount to avoid errors on rough edges --
+an image of the footprints is convolved with a disk of radius 3
+pixels.  Peaks are assigned to the footprints in which they are
+contained (note by definition all peaks must be located in a
+footprint).
+For any peak which is not the brightest peak in that footprint it is
+possible to reach the brightest peak by following the highest valued
+pixels between the two peaks.  The lowest pixel along this path is the
+{\em key col} for this peak (as used in topographic descriptions of a
+mountain).  If the key col for a given peak is less than
+\code{FOOTPRINT_CULL_NSIGMA_DELTA} (4.0) sigmas below the peak of
+interest, the peak is considered to be {\em locally insignificant} and
+removed from the list of possible detections (see
+Figure~\ref{fig:peaks}).  In the vicinity of a saturated star, the
+rule is somewhat more agressive as the flat-topped or structured
+saturated top of a bright star may appear as multiple peaks with
+highly significant cols between them.  However, this is an artifact of
+the proximity to saturation.  In this regime, we require the col to
+also be a fixed fraction (5\%) of the saturation below the peak to
+avoid being marked as locally insignificant.
+\subsubsection{Centroid and higher-order Moments}
+\label{sec:moments}
+\begin{figure}[htbp]
+  \begin{center}
+  \includegraphics[width=\hsize]{{pics/FWHM.smooth.trend.ps1}.\plotext}
+  \caption{\label{fig:moments.window} Example of the biases encountered when measuring the second
+    moments.  A simulated image was generated using the PS1 PSF
+    profile.  Each panel corresponds to a different value of
+    $\sigma_w$, as marked.  The solid red line is the true FWHM of the
+    PSF used to generate the stars.  The blue solid line is the FWHM
+    of the window function ($2.35\sigma_w$).  The gray dots are the
+    FWHM derived from the measured second moments for stars in the
+    image.  The dotted blue line is the target (65\% of the window
+    function).  In this example, we would choose $\sigma_w$ between
+.5 and 0.8 arcseconds so the dotted blue line would match the
+    bright end of the gray dots.}
+  \end{center}
+\end{figure}
+Once a collection of peaks has been identified, a number of basic
+properties of the objects related to the first and second moments are
+measured.  Below, the second moments are used to select candidate
+stellar sources to be used in modeling the PSF.
+In order to measure the moments, it is necessary to define an
+appropriate aperture in which the moments are measured.  We also apply
+a ``window function'', down-weighting the pixels by a Gaussian of size
+$\sigma_W$ which is chosen to be large compared to the PSF size,
+$\sigma_{\rm PSF}$.  This
+window function reduces the noise of the measurement of the first and
+second moments by suppressing the noisy pixels at high radial distance
+as well as by reducing the contaminating effects of neighboring stars.
+The choice of the window function $\sigma_W$ and the aperture is an
+iterative process: for a given value of $\sigma_W$, the PSF stars will
+have a measured value of $\sigma_{\rm PSF}$ which is modified by the effect of
+the window function.  In addition, depending on the size of the window
+function compared to the true PSF size, the measured value of the PSF
+size, $\sigma_{\rm PSF}$, will be biased high or low depending on the
+signal-to-noise of the object.
+These effects are illustrated in Figure~\ref{fig:moments.window} using
+simulated data.  An image was generated with a PSF model matching the
+radial profile of the PS1 PSF model with a FWHM of 1.4 arcseconds.  As
+the window function $\sigma_W$ is increased, the measured FWHM for the
+bright simulated stars rises to meet the truth value.  For small
+values of $\sigma_W$, fainter stars are biased to low measured values
+of the FWHM.  For large values of $\sigma_W$, the faint stars are
+biased to higher values and the scatter increases.
+In a real image, we do not know the true value of the PSF size.  If we
+simply choose a very large window function and rely on bright stars,
+our estimate of the PSF size will be quite noisy.  Compounding this
+problem are the two additional facts that (1) we do not know which are
+the real stars (as opposed to bright galaxies or possible image
+artifacts) and (2) the brighter stars are themselves subject to
+additional biases due to saturation and other non-linear effects
+(c.f., ``the Brighter-Fatter'' effect, REF).  To make a robust
+choice for the window function $sigma_w$, we choose a value
+such that the measured value of $\sigma_{\rm PSF}$ is 65\% of
+$\sigma_w$.  The resulting second moment values are biased somewhat
+low (\approx 75\% of the truth value for the PS1 PSF profile), but are
+relatively unbiased as a function of brightness.
+To choose the value of $\sigma_W$, we try values of (1, 2, 3, 4.5, 6,
+, 12, 18) pixels $\approx$ (0.26, 0.51, 0.77, 1.15, 1.54, 2.3, 3.1,
+.6) arcseconds.  For each of these values, we then select candidate
+PSF stars based on the distribution of the measured $\sigma_{x,x},
+\sigma_{y,y}$ values.  For each test value of $\sigma_w$, we determine
+the ratio $f = \frac{\sigma_{x,x} + \sigma{y,y}}{2 \sigma_w}$, i.e.,
+the ratio of the window size to the observed PSF size.  We interpolate
+to find a value of $\sigma_W$ for which $f$ is expected to be 0.65.
+We call this value the \code{MOMENTS_GAUSS_SIGMA}.  We use an aperture
+with a radius of \code{PSF_MOMENTS_RADIUS} = 4$\times$
+\code{MOMENTS_GAUSS_SIGMA} to select the pixels for the measurement.
+Once \code{PSF_MOMENTS_SIGMA} has been determined, moments are
+measured as defined below.
 \begin{eqnarray}
+Det    & = & 4 C_{20} C_{02} - C_{11}^2 \\
+x_{min} & = & (C_{11} C_{01} - 2 C_{02} C_{10}) / Det \\
+y_{min} & = & (C_{11} C_{10} - 2 C_{20} C_{01}) / Det \\
+x_0      & = & \frac{1}{S} \sum_i (f_i - s_i)x_i w_i \\
+y_0      & = & \frac{1}{S} \sum_i (f_i - s_i)y_i w_i \\
+M_{xx}   & = & \frac{1}{S} \sum_i (f_i - s_i)(x_i - x_0)^2w_i \\
+M_{xy}   & = & \frac{1}{S} \sum_i (f_i - s_i)(x_i - x_0)(y_i - y_0)w_i \\
+M_{yy}   & = & \frac{1}{S} \sum_i (f_i - s_i)(y_i - y_0)^2w_i \\
+M_{xxx}  & = & \frac{1}{S} \sum_i (f_i - s_i)(x_i - x_0)^3w_i / r_i \\
+M_{xxy}  & = & \frac{1}{S} \sum_i (f_i - s_i)(x_i - x_0)^2(y_i - y_0)w_i / r_i \\
+M_{xyy}  & = & \frac{1}{S} \sum_i (f_i - s_i)(x_i - x_0)(y_i - y_0)^2w_i / r_i \\
+M_{yyy}  & = & \frac{1}{S} \sum_i (f_i - s_i)(y_i - y_0)^3w_i / r_i \\
+M_{xxxx} & = & \frac{1}{S} \sum_i (f_i - s_i)(x_i - x_0)^4w_i / r^2_i \\
+M_{xxxy} & = & \frac{1}{S} \sum_i (f_i - s_i)(x_i - x_0)^3(y_i - y_0)w_i / r^2_i \\
+M_{xxyy} & = & \frac{1}{S} \sum_i (f_i - s_i)(x_i - x_0)^2(y_i - y_0)^2w_i / r^2_i \\
+M_{xyyy} & = & \frac{1}{S} \sum_i (f_i - s_i)(y_i - y_0)(y_i - y_0)^3w_i / r^2_i \\
+M_{yyyy} & = & \frac{1}{S} \sum_i (f_i - s_i)(y_i - y_0)^4w_i / r^2_i
 \end{eqnarray}
+Applying error propagation to the above, we find:
+where $f_i$ is the flux in a pixel; $s_i$ is the local sky value for
+that pixel; $w_i$ is the value of the window function for the pixel;
+$S = \sum_i (f_i - s_i) w_i$ is the window-weighted sum of the source
+flux, used to re-normalize the moments; $r_i$ is the radius of a
+pixel, $\sqrt{(x_i - x_0)^2 + (y_i - y_0)^2}$; The sum is performed
+over all pixels in the aperture.  For the centroid calculation ($x_0,
+y_0$), the peak coordinate (see~\ref{sec:peaks}) is used to define the
+aperture and the window function; for higher order moments, the
+centroid is used to center the window function.
+If the measured centroid coordinates ($x_0, y_0$) differs from the
+peak coordinates be a large amount (\code{MOMENT_RADIUS}), then the
+peak is identified as being of poor quality and is rejected.  In
+both of these cases, it is likely that the `peak' was identified in a
+region of flat flux distribution or many saturated or edge pixels.
+In addition to the moments above, a preliminary Kron radius and flux
+are also calculated at this stage.  In this analysis, the 1st and
+half-radial moments are calculated:
 \begin{eqnarray}
+\sigma_{Det}^2  & = & \sigma_{11}^2 (4 C_{11}^2) + \sigma_{20}^2 (16 C_{02}^2) + \sigma_{02}^2 (16 C_{20}^2) \\
+\sigma_{xn}^2   & = & \sigma_{11}^2 C_{01}^2 + \sigma_{01}^2 C_{11}^2 + \sigma_{02}^2 (4 C_{10}^2) + \sigma_{10}^2 (4 C_{02}^2) \\
+\sigma_{yn}^2   & = & \sigma_{11}^2 C_{10}^2 + \sigma_{10}^2 C_{11}^2 + \sigma_{20}^2 (4 C_{01}^2) + \sigma_{01}^2 (4 C_{20}^2) \\
+\sigma_{x}^2    & = & x^2 (\sigma_{xn}^2 / xn^2 + \sigma_{Det}^2 / Det^2) \\
+\sigma_{y}^2    & = & y^2 (\sigma_{yn}^2 / yn^2 + \sigma_{Det}^2 / Det^2) \\
+M_r & = & \frac{1}{S} \sum_i (f_i - s_i)r_i \\
+M_h & = & \frac{1}{S} \sum_i (f_i - s_i)\sqrt{r_i}
 \end{eqnarray}
+Note that the window function is not applied in the calculation of
+these moments.
+The Kron radius \citep{1980ApJS...43..305K} is defined the be
+.5$\times$ the first radial moment.  The Kron flux is the sum of
+(sky-subtracted) pixel fluxes within the Kron radius.  We also
+calculate the flux in two related annular apertures: the Kron inner
+flux is the sum of pixel values for the annulus $R_1 < r < 2.5 R_1$,
+while the Kron outer flux is the sum of pixel values for $2.5 R_1 < r
+< 4 R_1$.  The first radial moment is limited at the low and high ends
+by $R_{\rm min} < M_r < R_{\rm max}$ where $R_{\rm min}$ is the first
+radial moment of the PSF stars, or 0.75$\times$
+\code{MOMENTS_GAUSS_SIGMA} if that cannot be determined.  $R_{\rm
+  max}$ is set to \code{PSF_MOMENTS_RADIUS}, the size of the moments
+aperture.
 \subsection{PSF Determination}
 …
 for example, a 2-D elliptical Gaussian:
 \begin{eqnarray}
+f(x,y) & = & I_o exp (-z) + S  \\
+    R  & = & \frac{(x - x_o)^2}{2\sigma_x^2} + \frac{(y -
+    y_o)^2}{2\sigma_y^2} + \sigma_{\rm xy}(x - x_o)(y - y_o)
+f(x,y) & = & I_o e^{-z} + S  \\
+    z  & = & \frac{x^2}{2\sigma_x^2} + \frac{y^2}{2\sigma_y^2} + \sigma_{\rm xy} x y \\
+    x  & = & x_{\rm ccd} - x_o \\
+    y  & = & y_{\rm ccd} - y_o
 \end{eqnarray}
 The object model will have a variety of model parameters, in this case
 …
 The point-spread-function (PSF) of an image describes the shape of all
 unresolved objects in the image.  In a typical image, the shape of
+point sources is not well described by a single functional form;
+rather, the shape will vary as a function of position in the image.
+The PSF model therefore must describe the parameter variation as a
+function of the position of the object on the image.  Note that the
+object model consists of a certain number of parameters which are
+defined by the PSF model, and another set of parameters which are
+independent from object to object.  For the case of the elliptical
+Gaussian model, the PSF parameters would be the shape terms
+($\sigma_x, \sigma_y, \sigma_{\rm xy}$) while the independent
+parameters would be the centroid, normalization and local sky values
+($x_o, y_o, I_o, S$).  PSPhot uses a 2-D polynomial to specify the
+variation in the PSF parameters as a function of position in the
+image.  In the case of the elliptical Gaussian, this implies that the
+parameters are each a function of the object centroid coordinates:
+point sources is not well described by a single function.  Instead,
+the shape will vary as a function of position in the image.  The PSF
+model therefore must describe the parameter variation as a function of
+the position of the object on the image.  Note that the object model
+consists of a certain number of parameters which are defined by the
+PSF model, and another set of parameters which are independent from
+object to object.  For the case of the elliptical Gaussian model, the
+PSF parameters would be the shape terms ($\sigma_x, \sigma_y,
+\sigma_{\rm xy}$) while the independent parameters would be the
+centroid, normalization and local sky values ($x_o, y_o, I_o, S$).
+Thus these parameters are each a function of the object centroid
+coordinates:
 \begin{eqnarray}
 \sigma_x    & = & f_1(x,y) \\
 …
 \sigma_{xy} & = & f_3(x,y) \\
 \end{eqnarray}
+PSPhot represents the variation in the PSF parameters as a function of
+position in the image in two possible ways, specified by the
+configuration.  The first option is to use a 2-D polynomial which is
+fitted to the measured parameter values across the image.  The second
+option is to use a grid of values which are measured for objects
+within a subregion of the image.  In the latter case, the value at a
+specific coordinate in the image is determined by interpolation
+between the nearest grid points.  The order of the polynomial or the
+sampling size of the grid is dynamically determined depending on the
+number of available of PSF stars.  In the case of the PV3 analysis,
+the grid of values was used, with a maximum of $6\times 6$ samples per
+GPC1 chip image.  For the earlier PV2 analysis, the maximum grid
+sampling was $3\times 3$ per GPC1 chip image.  For the PV1 analysis,
+the polynomial representation was used, with up to 3rd order terms.
+The higher order representation was used for PV3 in order to follow
+some of the observed PSF variations in the images
+% XXX specify the rule for the polynomial order and grid scale
+% XXX discuss the improvements in the astrometric modeling PV1 - PV3
 PSPhot uses a single structure to represent the object model and
 …
 registered as part of the model function code.  Another function is
 then used to return the appropriate function for a specific model
 type.  For example, the \code{psModelLookup\_GetFunction} will return
+type.  For example, the \code{psModelLookup_GetFunction} will return
 the \code{psModelLookup} function for a given model type.  This
 mechanism makes it very easy to add new model functions into the
 …
 the 4 independent parameters.
-\note{the code may also require that two of the PSF-like parameters
-represent the shape in some way}.
 \subsubsection{PSF Candidate Object Selection}
 …
 their peaks, as well as an approximate signal-to-noise ratio.  All
 objects with a S/N ratio greater than a user-defined parameter
+(\code{PSF\_SHAPE\_NSIGMA} ???) are selected by PSPhot, though objects
+which have more than a certain number of saturated pixels are excluded
+at this stage.  PSPhot then examines the 2-D plane of $\sigma_x,
+\sigma_y$ in search of a concentrated clump of objects.  To do this,
+it constructs an artificial image with pixels representing the value
+of $\sigma_x, \sigma_y$, using a user-defined scale for the size of a
+pixel in this artificial image (note that the units of the $\sigma_x,
+\sigma_y$ plane are the size of the second-moment in pixels in the
+original image).  A typical value for the bin size is approximately
+.1 image pixels.  The binned $\sigma_x, \sigma_y$ plane is then
+examined to find a peak which has a significance greater than XXX.
+Unless the image is extremely sparse, such a peak will be well-defined
+and should represent the objects which are all very similar in shape.
+Other objects in the image will tend to land in very different
+locations, failing to produce a single peak.  To avoid detecting a
+peak from the unresolved cosmic rays, objects which have
+(\code{PSF_SHAPE_NSIGMA} = 20.0) are selected by PSPhot, though
+objects which have more than a certain number of saturated pixels are
+excluded at this stage.  PSPhot then examines the 2-D plane of
+$\sigma_x, \sigma_y$ in search of a concentrated clump of objects (see
+Figure~\ref{fig:moment.class}).  To
+do this, it constructs an artificial image with pixels representing
+the value of $\sigma_x, \sigma_y$, using a user-defined scale for the
+size of a pixel in this artificial image (note that the units of the
+$\sigma_x, \sigma_y$ plane are the size of the second-moment in pixels
+in the original image).  A typical value for the bin size is
+approximately 0.1 image pixels.  The binned $\sigma_x, \sigma_y$ plane
+is then examined to find a peak which has a significance greater than
+XXX.  Unless the image is extremely sparse, such a peak will be
+well-defined and should represent the objects which are all very
+similar in shape.  Other objects in the image will tend to land in
+very different locations, failing to produce a single peak.  To avoid
+detecting a peak from the unresolved cosmic rays, objects which have
 second-moments very close to 0 are ignored.  The only danger is if the
 PSF is very small and too many of these objects are rejected as cosmic
 …
 the image.
+\begin{figure}[htbp]
+  \begin{center}
+  \includegraphics[width=\hsize]{{pics/moment.class}.\plotext}
+  \caption{\label{fig:moment.class} Illustration of PSF star selection using the FWHM derived
+    from the second moments in $X_{\rm ccd}$ and $Y_{\rm ccd}$
+    directions.  The dominant clump is located in this diagram.
+    Galaxies tend to have a range of sizes and thus spread out above
+    the stars.  Cosmic rays also have a range of sizes, with one
+    dimension smaller than the PSF.  The red circle represents the PSF
+    star candidates. }
+  \end{center}
+\end{figure}
 \subsubsection{PSF Candidate Object Model Fits}
+% \note{link to psLibADD}
 All candidate PSF objects are then fitted with the selected object
 model, allowing all of the parameters (PSF and independent) to vary in
+the fit.  PSPhot uses the Levenberg-Marqardt process for the
+non-linear fitting.  Non-linear fitting can be very computationally
+intensive, particularly for if the starting parameters are far from
+the minimization values.  PSPhot uses the first and second moments to
+make a good guess for the centroid and shape parameters for the PSF
+models.  \note{still true? In order to minimize the impact of close
+  neighbors, the variance values used in the fit are enhanced by a
+  fraction of the deviation of the particular pixel value from the
+  model guess.}  Any objects which fail to converge in the fit are
+flagged as invalid.
+\note{does the variance enhancement introduce too much bias?}
+\note{discuss the convergence criteria, model parameter guesses}
+the fit.  PSPhot uses the Levenberg-Marquardt minimization technique
+for the non-linear fitting.  Non-linear
+fitting can be very computationally intensive, particularly for if the
+starting parameters are far from the minimization values.  PSPhot uses
+the first and second moments to make a good guess for the centroid and
+shape parameters for the PSF models.  Any objects which fail to
+converge in the fit are flagged as invalid.
 For the resulting collection of object model parameters, the
 PSF-dependent parameters of the models are all fitted as a function of
 position to a 2-D polynomial.  The order of this polynomial is (should
+be?) a user-defined parameter.  The fitting process for these
+polynomials is iterative, and rejects the $3-\sigma$ outliers in each
+of three passes.  This fitting technique results in a robust
 measurement of the variation of the PSF model parameters as a function
+of position without being excessively biased by individual objects
+which fail drastically.  Objects whose model parameters are rejected
+by this iterative fitting technique are also marked as invalid and
 ignored in the later PSF model fitting stages.
+position to a 2-D polynomial.  The order of this polynomial is a
+user-defined parameter.  The fitting process for these polynomials is
+iterative, and rejects the $3-\sigma$ outliers in each of three
+passes.  This fitting technique results in a robust measurement of the
+variation of the PSF model parameters as a function of position
+without being excessively biased by individual objects which fail
+drastically.  Objects whose model parameters are rejected by this
+iterative fitting technique are also marked as invalid and ignored in
+the later PSF model fitting stages.
 All of the PSF-candidate objects are then re-fitted using the PSF
 …
 correction is judged to be the best model.
-\subsubsection{Basic Deblending}
-The collection of identified peaks is examined to find peaks which are
-'blended', that is, they are close enough together to make the
-analysis of one of the sources difficult if performed in isolation.
-Saturated stars also result in additional peaks which are likely to be
-invalid; it is useful to restrict a saturated star to a single primary
-position with associated neighboring peaks.
-The deblending process first searches for any peaks which are within
-the image cell of another peak.  All such groups are examined,
-starting with the brightest source.  An isophot is found about the
-primary peak which is at least \code{DEBLEND\_SKY\_NSIGMA} times the sky
-sigma above the local background and which is otherwise
-\code{DEBLEND\_PEAK\_FRACTION} of the primary peak central pixel flux.
-Any secondary sources which are contained within this isophot are
-considered to be blended peaks associated with the primary peak.
 \subsection{Bright Source Analysis}
+After a PSF model has been determined, PSPhot performs the analysis of
+the bright objects in the image.  In this stage, all of the objects
+with an estimated signal to noise (based on the moments analysis)
+greater than a user-set threshold are analysed and subtracted from the
+image.  An optional successive stage then finds fainter sources and
+measures them as well (see Faint Source Analysis,
+Section~\ref{faintsources}).  In the bright source analysis stage, two
+major varients are available.  In the primary version, all objects are
+examined (in decending order of brightness) and an appropriate models
+is determined for each object which is then subtracted; in the
+alternate version, the objects are examined (in decending order of
+brightness) and the PSF-like objects subtracted first, then the
+extended objects are analysed on a second pass.
+%% \subsubsection{Very Bright Stars}
+%%
+%% The PSF modeling code fails to fit the wings of highly saturated stars
+%% if the core of the star is too contaminated by saturated pixels. For
+%% stars with estimated instrumental magnitudes brighter than XXX, we fit
+%% and subtract a radial profile modeled with a spline (?).
 \subsubsection{Fast Ensemble PSF Fitting}
 …
 convenient to subtract off all of the sources, at least as well as
 possible at this stage.  We make the assumption that all sources are
 PSF-like.  We also assume their position can be taken as the peak of a
+D quadratic function fitted to the peak pixel and its surrounding 8
 pixels.  A single linear fit is used to simultaneously measure all
+source fluxes.  Since the local sky has been subtracted, this
+measurement assumes the local sky is zero.
+PSF-like.  If the centroid of the source has been determined, we use
+this value for its position; otherwise, we use the interpolated
+position of the peak. A single linear fit is used to simultaneously
+measure all source fluxes.  Since the local sky has been subtracted,
+this measurement assumes the local sky is zero.  We can write a single
+$\chi^2$ equation for this image:
 \[
 \chi^2 = \sum_{\rm pixels} (F_{x,y} - \sum_{\rm sources} A_i PSF[x,y])^2
+\chi^2 = \sum_{\rm pixels} (F_{x,y} - \sum_{\rm sources} A_i P[x_0,y_0])^2
 \]
+where $F_{x,y}$ is image flux for each pixel, $P[x_0,y_0]$ is the PSF
+model realized at the position of source $i$, and $A_i$ is the
+normalization for the source.
 Minimizing this equation with respect to each of the $A_i$ values
 …
 \[ M_{i,j} \bar{A_i} = \bar{F_j}\]
 where $\bar{A_i}$ is the vector of $A_i$ values, the elements of
 $M_{i,j}$ consist of the dot product of the unit-flux PSF for source
 $i$ and source $i$, and $\bar{F_j}$ is the dot product of the
 unit-flux PSF for source $i$ with the pixels corresponding to source
 $i$.  The dot products are calculated only using pixels within the
+$M_{i,j}$ consist of the dot products of the unit-flux PSF for source
+$i$ and source $j$, and $\bar{F_j}$ is the dot product of the
+unit-flux PSF for source $j$ with the pixels corresponding to source
+$j$.  The dot products are calculated only using pixels within the
 source apertures.  Since most sources have no overlap with most other
 sources, this matrix equation results in a very sparse, mostly
 diagonal square matrix.  The dimension is the number of sources,
+likely to be 1000s or 10,000s.  Such a matrix does not lend itself to
+direct inversion.  However, an interative solution quickly yields a
+result with sufficient accuracy.  In the iterative solution, a guess
+at the solution is made; the guess is multiplied by the matrix, and
+the result compared with the observed vector $\bar{F_j}$.  The
+difference is used to modify the initial guess. This proces is
+repeated several times to achieve a good convergence.
+likely to be 1000s or 10,000s.  Direct inversion of the matrix would
+be computationally very slow.  However, an interative solution quickly
+yields a result with sufficient accuracy.  In the iterative solution,
+a guess at the solution $\bar{A}$ is made assuming $M_{i,j}$ is purely
+diagonal; the guess is multiplied by $M_{i,j}$, and the result
+compared with the observed vector $\bar{F_j}$.  The difference is used
+to modify the initial guess.  This proces is repeated several times to
+achieve a good convergence.
 Once a solution set for $A_i$ is found, all of the objects are
 …
 \subsubsection{PSF Model applied to detected objects}
+% \note{review the discussion below}
 Once a PSF model has been selected for an image, PSPhot attempts to
 fit all of the detected objects, above a user-defined signal-to-noise
 ratio (\note{KEYWORD}) with the PSF model.  For these fits, the
+dependent parameters are fixed by the PSF model and only the 4
+independent object model parameters are allowed to vary in the fit.
+PSPhot again uses the Levenberg-Marqardt process for the non-linear
+fitting.  The objects are fitted in their S/N order, starting with the
 brightest and working down to the user-specified limit.
 Once a solution has been achieved, PSPhot attempts to judge the
+quality of the PSF model as a representation of the object shape.  To
+do this, it calculates the next step of the minimization {\em allowing
+the shape parameters to vary}.  This step, essentially the
 Gauss-Newton minimization distance from the current local minimum,
+ratio with the PSF model.  For these fits, the dependent parameters
+are fixed by the PSF model and only the 4 independent object model
+parameters are allowed to vary in the fit.  PSPhot again uses
+Levenberg-Marquardt minimization for the non-linear fitting.  The
+objects are fitted in their S/N order, starting with the brightest and
+working down to the user-specified limit.
+Once a solution has been achieved for an object, PSPhot attempts to
+judge the quality of the PSF model as a representation of the object
+shape.  To do this, it calculates the next step of the minimization
+{\em allowing the shape parameters to vary}.  This step, essentially
+the Gauss-Newton minimization distance from the current local minimum,
 should be very small if the object is well represented by the PSF, but
 large if the PSF is not a good representation of the object flux.  The
 …
 elliptical Gaussian, these two parameters are $\sigma_x$ and
 $\sigma_y$.  For a generic model, the shape parameters may be defined
+differently, but the should always be two parameters which scale the
+object size in two dimensions (what about a polar-coordinate form?)
+Currently, PSPhot requires the two relevant shape parameters to be the
+first two dependent parameters in the list of model parameters (ie,
+parameters 4 \& 5).
+differently, but there should always be two parameters which scale the
+object size in two dimensions.  Currently, PSPhot requires the two
+relevant shape parameters to be the first two dependent parameters in
+the list of model parameters (ie, parameters 4 \& 5).
 The expected distribution of the variation of the two shape parameters
 will be a function of the signal-to-noise of the object in question
 and the value of the shape parameter itself.  The expected standard
 deviation on the shape parameter is, eg, $\sigma_x / {\rm SN}$.  If
+deviation on the shape parameter is, eg, $\sigma_x / {\rm S/N}$.  If
 the object is well-represented by the PSF, then the shape parameter
 values should be close to their minimization value.  We can thus ask,
 …
 distribution, the remaining flux should be below 1 $\sigma$
 significance.  In practice the cores of bright stars are poorly
+represented and may have larger residual significance. \note{in future
+work, we may choose to enhance the variance to minimize detection of
+objects in the residuals of brighter objects}.
+represented and may have larger residual significance.
 \subsubsection{Blended Sources}
 …
 available non-PSF model or models.
+\note{better description of the acceptance criteria; the FLT model is
+  tried before the DBL is accepted or rejected}.
+\subsubsection{Source Size Assessment}
+\label{sec:source.size}
+After the PSF model has been fitted to all sources, and the Kron flux
+has been measured for all sources, PSPhot uses these two measurements,
+along with some additional pixel-level analysis, to determine the size class
+of the object.  If the object is large compared to a PSF, it is
+considered to be {\em extended} and will be
+fitted with a galaxy model (or possibly another type of extended
+source model in special cases).  If the object is small compared to a
+PSF, it is considered to be a {\em cosmic ray} and masked.
+Extended sources are identified as those for which the Kron magnitude is
+significantly brighter than the PSF magnitude when compared to a PSF
+star.  The value $dMagKP = m_{\rm Kron} - m_{\rm PSF}$, the difference between the PSF
+and Kron magnitudes, is calculated for each object.  The median of
+$dMagKP$ is calculated for the PSF stars.  This median is subtracted
+from $dMagKP$ for each star.  The result is divided by the quadrature
+error of the PSF and Kron magnitudes and called \code{extNsigma}.  If
+\code{extNsigma} is larger than \code{PSPHOT.EXT.NSIGMA.LIMIT} (3.0),
+the object is considered to be extended.
+Cosmic Rays are identified by a combination of the Kron magnitude and
+the second-moment width of the object in the minor axis direction.
+The second-moment in the minor axis direction is calculated from
+$M_{xx}, M_{xy}, M_{yy}$ as follows:
+\[
+M_{\rm minor} = \frac{1}{2}(M_{xx} + M_{yy}) - \frac{1}{2}\sqrt{(M_{xx} - M_{yy})^2 + 4 M_{xy}^2}
+\]
+If $M_{\rm minor} < 1.2$ pixels$^2$ and the instrumental Kron
+magnitude is $< -5.5$, then the object is identified as a cosmic ray
+and the associated pixels are masked.
 \subsubsection{Non-PSF Objects}
 …
 PSPhot will use the user-selected galaxy model to attempt the galaxy
 model fits.  In the configuration system, the keyword \code{GAL\_MODEL}
+model fits.  In the configuration system, the keyword \code{GAL_MODEL}
 is set to the model of interest.  All suspected extended objects are
 fitted with the model, allowing all of the parameters to float.  The
 …
 not modified.
-\note{we have no code yet to select the best of several models for a
-  given objects.  The relative value of the Chi-Square is the obvious
-  test in this case}.
 \subsection{Faint Sources}
-\note{this is not done : should use the ensemble PSF fitting to fit
-  just the new significant peaks}
 After a first pass through the image, in which the brighter sources
 …
 The objects which are measured in this faint-object stage are clearly
+low significance detections.  A typical threshold for the bright
+object analysis would S/N of 5 - 10.  The lower limit cutoff for the
+faint object analysis would typically be S/N of 2 - 4.  In this stage,
+PSPhot can perform one of three types of analysis.  The difference
+between these options is one of speed vs detail.
+In the first option, PSPhot can repeat the analysis described above in
+sections XXX and XXX, performing a PSF fit followed by a non-PSF fit
+to the objects which are not PSF-like, and subtracting them.  The
+advantage of this option is that the faint objects are treated
+identically to the bright objects, and all potential parameters are
+measured, even for marginally extended sources.  The disadvantage of
+this option is that the low signal-to-noise of the objects in this
+stage limits the amount of information which can be extracted from
+them.  The marginal gain may not be worth the large expense of
+processing time.
+In the second option, PSPhot can perform only the PSF model fit to the
+remaining peaks, but ignore any further questions of the shape of the
+objects.  The advantage of this option is that it is substantially
+faster than performing the more complex (and less stable)
+multi-parameter non-linear fits on all faint objects.  On the
+downside, less information is learned about the objects.
+Finally, PSPhot can simply ignore the fitting process and instead
+glean information about the fainter sources on the basis of the peak
+characteristics.  In this option, the image is smoothed with the PSF
+model, and the peak for each object is measured.  The peak flux and
+the local peak curvature theoretically give sufficient information to
+recover the object flux, the centroid coordinates, and the centroid
+errors.  The advantage of the stage is speed, especially for the very
+faintest levels: if the lower limit is not sufficiently faint, the
+time spent in performing the smoothing (3 FFTs) cannot make up for the
+time which would have been spent applying the PSF model to the peaks.
+The downside of this method is an increased sensitivity to the local
+sky model (the local sky must be correctly subtracted) and fewer
+constraints on the quality of the detection (no Chi-Square is
+measured, for example).
+\note{In the ideal case, if we were only interested in detecting PSFs,
+and we had a good model for the PSF, we could optimally find the
+sources by smoothing the image and the variance image with the PSF model.
+\em write out the description of Nick's optimal PSF finding}.
+PSPhot allows the user to select between these three options for the
+analysis of the faint sources.  Three separate user-controlled
+signal-to-noise ratio limits are defined.  One specifies the depth to
+which the PSF / non-PSF analysis is performed.  A second (which must
+be smaller) specifies the depth to which only the PSF is fitted.  A
+third specifies the depth to which the analysis is performed using on
+the peak statistics.  If two of these are identical, then certain of
+these options are simply skipped.  For example, if the peak analysis
+threshold is set to the same value as the PSF-only threshold, no peak
+analysis is performed.
+low significance detections.  The PV3 threshold for the bright object
+analysis is a signal-to-noise of 20.  The lower limit cutoff for the
+faint object analysis in PV3 is a signal-to-noise of 5.0.  Objects
+detected in the faint object stage are fitted with the PSF model using
+the linear, ensemble fitting process.
 \subsection{Aperture Correction Measurement}
 …
 number of very bright stars is limited in any image, and of course the
 brighter stars are more likely to suffer from non-linearity or
+saturation.
+\note{this discussion sucks: put in some more details of my point:
+  amplitude of systematic vs random sky errors}
+How important is this effect?  Consider a typical bright object with a
+flux of (say) 40,000 counts in an image of background 1000 counts per
+pixel, with FWHM of 4 pixels.  In principle, the flux of this object
+should be measurable with an accuracy of roughly 0.57\%
+($\frac{\sqrt{40000 + 1000 \times 12}}{40000}$).  However, the
+measurement of the sky is limited at some finite level by Poisson
+statistics.  If we are required to use an aperture of (say) 25 pixels
+in radius (eg, 5 arcseconds for an 0.2 arcsec / pixel detector), and
+we have an annulus of twice this radius to measure the local sky, then
+we will have an error of XXX.
+\note{outline the variation of {\em ApResid} as a function of
+magnitude}.
+PSPhot measures the aperture correction ({\em ApResid}) for every PSF
+candidate object, then calculates the trend of this correction as a
+function of the magnitude.  This trend is fitted with a line.  The
+resulting function can be used to determine the effective aperture
+correction for an infinite flux object and the average bias inherent
+in the sky measurement for the image.  The scatter of the
+PSF-candidate object measurements about this trend is a measure of how
+well we can measure photometry from the image by applying the specific
+PSF model.  The slope of this trend is a measure of the bias in the
+local sky measurment for each object.  In principal, the measured sky
+levels could be modified by this bias.  More generally, the measured
+bias in a collection of images could be used to improve the model
+fitting or sky fitting portion of the software the remove the bias
+term.
+saturation.  PSPhot measures the aperture correction ({\em ApResid})
+for every PSF candidate object and applies this correction to the PSF
+model photometry.
+% How important is this effect?  Consider a typical bright object with a
+% flux of (say) 40,000 counts in an image of background 1000 counts per
+% pixel, with FWHM of 4 pixels.  In principle, the flux of this object
+% should be measurable with an accuracy of roughly 0.57\%
+% ($\frac{\sqrt{40000 + 1000 \times 12}}{40000}$).  However, the
+% measurement of the sky is limited at some finite level by Poisson
+% statistics.  If we are required to use an aperture of (say) 25 pixels
+% in radius (eg, 5 arcseconds for an 0.2 arcsec / pixel detector), and
+% we have an annulus of twice this radius to measure the local sky, then
+% we will have an error of XXX.
+%
+% \note{outline the variation of {\em ApResid} as a function of
+% magnitude}.
+%%% PSPhot measures the aperture correction ({\em ApResid}) for every PSF
+%%% candidate object, then calculates the trend of this correction as a
+%%% function of the magnitude.  This trend is fitted with a line.  The
+%%% resulting function can be used to determine the effective aperture
+%%% correction for an infinite flux object and the average bias inherent
+%%% in the sky measurement for the image.  The scatter of the
+%%% PSF-candidate object measurements about this trend is a measure of how
+%%% well we can measure photometry from the image by applying the specific
+%%% PSF model.  The slope of this trend is a measure of the bias in the
+%%% local sky measurment for each object.  In principal, the measured sky
+%%% levels could be modified by this bias.  More generally, the measured
+%%% bias in a collection of images could be used to improve the model
+%%% fitting or sky fitting portion of the software the remove the bias
+%%% term.
 PSPhot allows a collection of PSF model functions to be tried on all
 …
 value for the ApResid scatter is then used by PSPhot as the best PSF
 model for this image.  The number of models to be tested is specified
 by the configuration keyword \code{PSF\_MODEL\_N}.  The configuration
 variables \code{PSF\_MODEL\_0}, \code{PSF\_MODEL\_1}, through
 \code{PSF\_MODEL\_N - 1} specify the names of the models which should be
+by the configuration keyword \code{PSF_MODEL_N}.  The configuration
+variables \code{PSF_MODEL_0}, \code{PSF_MODEL_1}, through
+\code{PSF_MODEL_N - 1} specify the names of the models which should be
 tested.
+\subsubsection{Types of Object / PSF models currently available}
+\note{the discussion of the model types needs to be extended}
+Several likely PSF model classes are available within \code{psphot}:
 \begin{itemize}
+\item GAUSS  : Pure elliptical Gaussian
+\item PGAUSS : polynomial approximation to a Gaussian (PGAUSS)
+\item QGAUSS : power law with variable exponential term
+\item SGAUSS :
+\item Gaussian : $f = I_0 e^{-z}$
+\item Pseudo-Gaussian : $f = I_0 (1 + z + \frac{1}{2} z^2 + \frac{1}{6} z^3)^{-1}$ \code{[PGAUSS]}
+\item Variable Power-Law : $f = I_0 (1 + z + z^{\alpha})^{-1}$ \code{[RGAUSS]}
+\item Steep Power-Law : $f = I_0 (1 + \kappa z + z^{2.25})^{-1}$ \code{[QGAUSS]}
+\item PS1 Power-Law : $f = I_0 (1 + \kappa z + z^{1.67})^{-1}$ \code{[PS1_V1]}
 \end{itemize}
+\note{discuss the stability issues with the galaxy model(s)}
+where $z \propto r^2$ ($z = \frac{x^2}{2\sigma_x^2} +
+\frac{y^2}{2\sigma_y^2} + \sigma_{\rm xy} x y $).  The Pseudo-Gaussian
+is a Taylor expansion of the Gaussian and is used by Dophot
+\citep{1993PASP..105.1342S}.  The latter profiles are similar to the
+Moffat profile form \citep{1969AA.....3..455M,1983AA...126..278B},
+with small differences.  For the PS1 GPC1 analysis, we used the
+\code{PS1_V1} model, which we found by experimentation to match well
+to the observed profiles generated by PS1.
+Figure~\ref{fig:radial.profiles} shows example radial profiles for
+moderately bright stars in fairly good (0.9 arcsec) and poor (2.2
+arcsec) seeing.  Using a fixed power-law exponent results in somewhat
+faster profile fitting compared to the variable power-law exponent
+model.
+% moffat : 1969A&A.....3..455M
+% buonanno : 1983A&AS...51...83B
+\begin{figure}[htbp]
+  \begin{center}
+  \includegraphics[width=\hsize]{{pics/radial.profiles}.\plotext}
+  \caption{\label{fig:radial.profiles} Radial profiles of stellar images from PS1.  These two
+    profiles illustrate the radial trend of the PS1 PSFs for a star
+    with FWHM 0.9 arcsec (red) and 2.2 arcsec (blue).  The black line
+    shows the PSF model with radial trend of the form $(1 + \kappa r^2 + r^{3.33})^{-1}$.}
+  \end{center}
+\end{figure}
+\subsection{Radial Profiles}
+Galaxies with regular profiles, such as elliptical galaxies and
+regular spiral galaxies, may be described as primarily a radial
+surface brightness profile, with additional structure acting as a
+perturbation on that profile.  For many galaxies, the azimuthal shape
+at a given isophotal level may be described as an elliptical contour.
+To first order, a galaxy may be well decribed with a single elliptical
+contour and radial profile.
+In order to facilitate the Petrosian photometry analysis below, PSPhot
+generates a radial profile for each suspected galaxy.  This analysis
+starts by generating a radial profile in 24 azimuthal segments.  Near
+the center of the galaxy, the profile is defined for radial
+coordinates in steps of 1 pixel, with the closest pixel values
+interpolated to that radial position.  Further from the center,
+profile is defined using the median of the pixels landing in an
+annular segment of size $\delta R = r \sin \theta$, rounded up to the
+nearest integer pixel value.  The median of all pixels within a
+rectangular approximation to the radial wedge is used.
+The resulting 24 radial profiles are subject to contamination from
+neighboring sources or to NAN values from masked pixels.  To clean the
+profiles, pairs of radial profiles from opposite sides of the source
+are compared.  Any masked values are replaced by the corresponding
+value in the other profile.  The minimum of both profiles is the kept
+for both profiles.  The result of this analysis is a set of profiles
+of the form $f_i(r_i)$.  In this case, $f_i$ is effectively the
+surface brightness for each radius in instrumental counts per pixel.
+The surface brightness profiles are then used to define the radial
+contour at a specific isophotal level.  This contour will be used to
+rescale the radial profiles into a single set of profiles normalized
+by the elliptical contour.  This contour is defined by determining the
+median radius for profile bins with surface brightness in the range
+$F_{\rm min} + 0.1 F_{\rm range}$ to $F_{\rm min} + 0.5 F_{\rm
+  range}$.  The result of this analysis is a value for the radius as a
+function of the angle for a well-defined surface brightness regime.
+We then determine the elliptical shape parameters for this elliptical
+contour: $R_{\rm major}, R_{\rm minor}, \theta$.  This ellipse is then
+used to redefine a single radial profile normalized by the elliptical
+contour:
+\[
+\rho = \sqrt{\frac{x^2}{S^2_{xx}} + \frac{y^2}{S^2_{yy}} + x y S_{xy}} \\
+\]
+The surface brightness values are sampled at a number of radial
+annuli, with the radii defined in the configuration
+(\code{RADIAL.ANNULAR.BINS.LOWER} \&
+\code{RADIAL.ANNULAR.BINS.UPPER}).  For each source, the resulting
+surface brightness profile is saved in the output cmf-file as an
+N-element value in the FITS table (\code{PROF_SB}).  The flux at each
+radial position and the fill-factor (fraction of pixels used to the
+total possible) as also saved as equal-length vectors in the FITS
+table (\code{PROF_FLUX} and \code{PROF_FILL}).  The values of the
+radial bins are saved in the cmf header (\code{RMIN_NN},
+\code{RMAX_NN}).
+% \note{these profiles are not saved in PSPS}
+\subsection{Petrosian Radii and Magnitudes}
+\cite{1976ApJ...209L...1P} defined an adaptive aperture based on a
+ratio of surface brightnesses.  The motivation is to define an
+aperture which can be determined for galaxies without significant
+biases as a function of distance.  Since surface brightness in a
+resolved object is conserved, using a ratio of surface brightness to
+define a spatial scale results in a spatial scale which is constant
+regardless of galaxy distance.
+To measure the Petrosian radius and flux, we start by defining a
+series of radial apertures with power-law spacing: $r_{i + 1} = 1.25
+r_i$.  We calculate the surface brightness for the annulus from $r_i -
+r_{i+1}$ by calculating the median of the values in the range $r_i /
+\sqrt{1.25}$ to $r_{i+1} \sqrt{1.25}$ and dividing the the effective
+area of the annulus corresponding to $r_i - r_{i+1}$.
+For any annulus $i$ spanning the radii $r_{\rm min}$ to $r_{\rm max} =
+\beta r_{\rm min}$, the
+Petrosian Ratio for that annulus is defined as the ratio of the
+surface brightness in the annulus to the average surface brigthness
+within $r_{\rm max}$.  The Petrosian Radius is defined to be $r_{\rm
+  max}$ for the annulus for which the Petrosian Ratio = 0.2, i.e., the
+point on the galaxy radial profile at which the surface brightness is
+\% of the average surface brightness at that point.
+We determine the Petrosian Radius for the galaxy by quadratic
+interpolation between the last two of the fixed annuli with Petrosian
+Ratio $> 0.2$ and the first annulus with Petrosian Ratio $< 0.2$.  In
+general, the Petrosian Ratio for a galaxy with a regular morphology
+(spiral or elliptical) is falling monotonically, so this interpolation
+is unambiguous.  However, irregular galaxy morphologies, noise, and/or
+significant masking can cause the Petrosian Ratio to have rises as
+well as drops.  We track the Petrosian Ratio until the value is no
+longer significant ($\sigma_{\rm Ratio} < 2 {\rm Ratio}$).  If the
+Petrosian Ratio drops below 0.2 for more than one radius, we choose
+the largest such radius.
+Once the Petrosian Radius has been determined, we can now measure the
+Petrosian Flux : this is defined to be the total flux within an
+aperture corresponding to 2 $\times$ the Petrosian Radius.  Using the
+Petrosian Flux, we can calculate two other interesting radii: $R_{50}$
+and $R_{90}$, the radii inside which 50\% and 90\% of the total
+Petrosian flux is contained.
+\subsection{Radial Profile Wings}
+We attempt to measure the radial profile of sources in order to find
+the radius at which the flux of the object is matches the sky.  In
+this analysis, a series of up to 25 radial bins with power-law spacing
+are defined and the flux of the object in each annulus is measured.
+The ``sky radius'' is defined to be the radius at which the (robust
+median) flux in the annulus is within 1 $\sigma$ of the local sky
+level.  If this limit is not reached before the slope of the flux from
+one annulus to the next is less than a user-defined limit, then the
+annulus at which the slope reaches this limit is used to define the
+sky radius.  These values are saved in the output smf / cmf files, but
+not sent to the PSPS.  The sky radius value is used below in the
+calculation of the kron magnitude.
+\subsection{Kron Magnitudes}
+Preliminary Kron radius and flux values \citep{1980ApJS...43..305K}
+are calculated soon after sources are detected (Section~\ref{sec:moments}).
+However, these preliminary values are not accurate due to the
+window-functions applied.  After sources have been characterized and
+the PSF model is well-determined, the Kron parameters are
+re-calculated more carefully.  In this version of the calculation, the
+image is first smoothed by Gaussian kernel with $\sigma = 1.7$ pixels,
+corresponding to a FWHM of 1.0\arcsec\ in the PS1 stack images.  Next,
+the Kron radius is determined in an iterative process: the first
+radial moment is measured using the pixels in an aperture 6$\times$
+the first radial moment from the previous iteration.  On the first
+iteration, the sky radius is used in place of the first radial moment.
+By default, 2 iterations are performed.  The Kron radius is defined
+the be 2.5$\times$ the first radial moment.  The Kron flux is the sum
+of pixel fluxes within the Kron radius.  We also calculate the flux in
+two related annular apertures: the Kron inner flux is the sum of pixel
+values for the annulus $R_1 < r < 2.5 R_1$, while the Kron outer flux
+is the sum of pixel values for $2.5 R_1 < r < 4 R_1$.
+Two details in the calculation above should be noted.  First, for
+faint sources, noise in the measurement of the 1st radial moment may
+result in an excessively small aperture for the successive
+calculations.  The window used for the calculations is constrained to
+be at least the size of the aperture based on the PSF stars
+(Section~\ref{sec:moments}).  At the other extreme, noise may make the radius
+grow excessively, resulting in an unrealistically low effective
+surface brightness.  The aperture is constrained to be less than a
+maximum value defined such that the minimum surface brightness is
+/2$times$ the effective surface brightness of a source detected at the
+$5\sigma$ limit.
+Second, the measurement of the 1st radial moment includes a filter to
+reduce contamination from outlier pixels.  Pairs of pixels on
+opposites sides of the central pixel are considered together.  The
+geometric mean of the two fluxes is used to replace the flux values.
+If the object has 180\degree\ symmetry, this operation has no impact.
+However, if one of the two pixels is unusually high, the value will be
+surpressed by the matched pixel on the other side.  This trick has the
+effect of reducing the impact of pixels which include flux from near
+neighbors.
+\subsection{Convolved Galaxy Model Fits}
+In the galaxy model fittting stage, sources which meet certain
+criteria are fitted with analytical models for galaxies.  Three
+traditional analytical galaxy models are implemented in \code{psphot}
+and used in the PV3 analysis:
+\begin{itemize}
+\item Exponential profile : $f = I_0 e^{-\rho}$
+\item DeVaucouleur profile \citep{1948AnAp...11..247D}: $f = I_0 e^{-\rho^{1/4}}$
+\item Sersic \citep{1963BAAA....6...41S} : $f = I_0 e^{-\rho^{1/n}}$
+\end{itemize}
+where $\rho$ is a normalized radial term: $\rho =
+\sqrt{\frac{x^2}{R^2_{xx}} + \frac{y^2}{R^2_{yy}} + x y R_{xy}}$.  The
+terms ($R_{xx}$, $R_{yy}$ , $R_{xy}$) describe the elliptical contour
+and the profile scale in all three models and the coordinates $x$ \&
+$y$ are determined relative to the centroids ($x,y = X_{\rm chip} -
+x_0, Y_{\rm chip} - y_0$).  Including the normalization ($I_0$) and a
+local sky value, the Exponential and DeVaucouleur profiles have 7 free
+parameters and the Sersic profile has the additional free parameter of
+the Sersic index $n$.
+In this stage, the galaxy model is convolved with an approximation to
+our best guess for the PSF model at the location of the galaxy.  For
+the PV3 analysis, all sources detected in the 'bright source' analysis
+step ($S/N > 20$) were fitted with all three galaxy models, unless (a)
+the morphological test identified the source as a likely cosmic ray
+(Section~\ref{sec:source.size}) or (b) the peak of the PSF profile was
+above the saturation limit for the chip \citep[see the discussion in
+][ regarding the masking of saturated pixels]{waters2017}.  Sources in
+the denser portions of the Galactic plane and bulge were not included
+in the analysis.  This restriction limited the total time spent on the
+galaxy modeling analysis at the expense of galaxy photometry in the
+plane (though Kron photometry is available for those objects).  The
+Galactic Plane region was defined by $|b| > b_{\rm min}$ where $b_{\rm
+  min} = b_0 + r_b e^{\frac{-l^2}{2 \sigma_b^2}}$.  For the PV3
+analysis, $b_0 = $20\degree, $r_b = $15\degree, $\sigma_b = $50\degree.
+%  \note{need a discussion of the detector saturation behavior
+% \note{more detail below?}
+Before the non-linear fitting may be performed, it is necessary to
+determine initial values for the parameters to be fitted.  For each of
+the three model types, the position determined from the PSF fitting
+analysis is used as the initial centroid $x_0,y_0$.  A guess for the
+terms ($R_{xx}$, $R_{yy}$ , $R_{xy}$) is generated based on the second
+moments.  The guess does not attempt to use the PSF model to adjust the
+($R_{xx}$, $R_{yy}$ , $R_{xy}$) values; it was found that such a guess
+tended to be too small and resulted in more iterations rather than
+fewer. The 1st radial moment (see
+\ref{sec:moments}) is used to estimate the effective radius of the
+model based on the results of Graham \& Driver (2005, Table 1).  They
+quantive the relationships between the first radial moment used to
+calculated a Kron Magnitude and the effective radius for different
+Sersic index values, $n$.  Since the Exponential and DeVaucouleur
+models are equivalent to Sersic models with $n$ = 1 and 4,
+respectively, this work can be used to generate the initial effective
+radius values for all 3 model types.  Once the effective radius is
+chosen, the second moments are used to define the aspect ratio and
+position angle of the elliptical contour.  The Kron flux is used to
+generate a guess for the normalization, applying an appropriate scale
+factor based on the ($R_{xx}$, $R_{yy}$ , $R_{xy}$) values, generated
+by integrating normalized Sersic models and determining the
+relationship between the central intensity and the integrated flux as
+a function of the Sersic index.
+The PSF-convolved galaxy model fitting analysys uses the
+Levenberg-Marquardt minimization method to determine the best fit.  In this
+process, the $\chi^2$ value to be minimized is:
+\[
+\chi^2 (\bar{a}) = \sum_p \frac{1}{\sigma_p^2} \left[I_p - M_p(\bar{a}) \otimes \mbox{PSF} \right]^2
+\]
+where $I_p$ represents the pixel values in the image (within some
+aperture) and $M_p(\bar{a})$ represents the unconvolved galaxy model, a
+function of a number of parameters $\bar{a}$, which is then convolved
+with the PSF model.
+We simplify this by defining:
+\begin{eqnarray}
+f_p (a_m)         & = & \frac{1}{\sigma_p} (I_p - M_p \otimes \mbox{PSF}) \\
+\end{eqnarray}
+To determine the minimization, we need the gradient and laplacian of
+$\chi^2$ with respect to the model parameters, $a_m$:
+\begin{eqnarray}
+\chi^2 (\bar{a})  & = & \sum_p f_p^2  \\
+\nabla   \chi^2  & = & \sum_p f_p \frac{\partial f_p}{\partial a_m} \\
+\nabla^2 \chi^2  & \approx & H_{m,n} \\
+H_{m,n}  & = & \sum_p \frac{\partial f_p}{\partial a_m} \frac{\partial f_p}{\partial a_n}
+\end{eqnarray}
+where we have approximated the Laplacian with the Hessian matrix,
+$H_{m,n}$ by dropping the second-derivatives (which are assumed to be
+a small perturbation).  Since
+\[
+\frac{\partial f_p}{\partial a_m} = -\frac{1}{\sigma_p}\frac{\partial M_p \otimes \mbox{PSF}}{\partial a_m}
+\]
+and since the order of the derivative and convolution may be
+exchanged, we can write these in terms of the convolved image of the
+model and the convolved images of the derivatives of the model $M_p$ with respect to the model parameters, $a_m$:
+\begin{eqnarray}
+\mathcal{M}_{p}   & = & M_p \otimes \mbox{PSF} \\
+\mathcal{M}^\prime_{p,m} & = & \frac{\partial M_p}{\partial a_m} \otimes \mbox{PSF} \\
+\nabla \chi^2    & = & -\sum_p \frac{I_p - \mathcal{M}_p}{\sigma_p} \mathcal{M}^\prime_{p,m} \\
+H_{m,n}  & = &  \sum_p \frac{1}{\sigma_p^2} \mathcal{M}^\prime_{p,m} \mathcal{M}^\prime_{p,n}
+\end{eqnarray}
+The gradient vector and Hessian matrix are used in the
+Levenberg-Marquardt minimization analysis using the standard
+techinique of determining a step from the current set of model
+parameters to a new set by solving the matrix equation:
+\[
+(1 + \lambda_{m,n}) H_{m,n} = \delta \nabla \chi^2
+\]
+where $\lambda_{m,n}$ is zero for $m \neq n$ and for $m = n$ set to be
+large when the last iteration produced a large change in the
+parameters compared to the local-linear expectation and small when the
+last change was small.  The iteration ends when the change in the
+parameters is small and/or the change in the $\chi^2$ value is small.
+In the analysis, convolved galaxy fit, the galaxy model image and the
+model derivative images must be convolved with the PSF at each
+iteration step.  To save computation time, this convolution is
+performed using a circularly symmetric approximation of the PSF model,
+with the PSF model scale size set to the average of the major and
+minor axis direction scale size of the full PSF model, with the same
+radial profile term as the PSF model.  The convolution is performed
+directly using the circular symmetry to reduce the number of
+multiplications performed: all points in the 2D circularly symmetric
+PSF model which have the same radial pixel coordinate can be evaluated
+in the convolution by summing up the corresponding pixels in the
+(galaxy model) image to be convolved before multiplying by the PSF
+model profile at that radial coordinate.  This approximation reduces
+the number of multiplications by a factor of near 8 for larger radii.
+For the small size of the PSF model used to convolve the galaxy model
+images, it was found that this direct convolution was faster than
+using an FFT-based convolution.
+% \note{(examples?)}
+For the Exponential and DeVaucouleur fits, all parameters are fitted
+in the non-linear minimization stage.  For the Sersic model, we do not
+fit the index within the Levenberg-Marquardt analysis.  Instead, we
+start with a coarse grid search over a range of possible index values,
+($n = 0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 6.0$) and a range of possible
+values for $R_{\rm eff}$ based on the value of $R_1$, the first radial
+moment.  For a given value of the Sersic index, the $R_{\rm eff}$ is
+related to the 1st radial moment by the scale factor specificy by
+Graham \& Driver.  We use the observed value of the 1st radial moment
+and try $R_{\rm eff}$ values of a factor of (0.8, 0.9, 1.0, 1.12,
+.25) times the value predicted by the Graham and Driver equation.
+For each of these steps, the aspect ratio and position angle are held
+constant and the normalization is determined to minimize the $\chi^2$.
+We next perform 3 Levenberg-Marquardt minimization fits allowing the
+shape parameters ($R_{xx}$, $R_{yy}$ , $R_{xy}$) and the normalization
+to be fitted, holding the centroid ($x_0, y_0$), Sersic index $n$, and
+sky constant.  In these fits, the index $n$ is set to the minimum
+value previously calculated as well as values halfway to the next, and
+previous, values in the grid above.  E.g., if the minimum fitted index
+value is 3.0, then the LMM fits are performed using $n$ = 2.5, 3.0, 3.5.
+The resulting $\chi^2$ values are then used to perform quadratid
+interpolation to find the index $n$ which produces the locally minium
+$\chi^2$ value.  Finally, this best-fit index value is held constant
+while Levenberg-Marquardt minimization is used to find the best fit
+values of all other parameters.
+% Graham & Driver : Graham A. W., Driver S. P.  2005, PASA 22, 118
+% DOI: https://doi.org/10.1071/AS05001
+\subsection{Convolved Radial Aperture Photometry}
+For some science goals, a well-measured color of a galaxy is more
+important than an accurate total magnitude.  In the case of PS1, the image
+quality variations for stacks of different filters presents a serious
+challenge for the determination of precise colors.  PSPhot determines
+a set of PSF-matched radial aperture flux measurements in order to
+minimize the impact of the stack image quality variations.
+In PSPhotStack, the stack analysis version of PSPhot, the 5 filter
+images are processed together.  After the PSF models have been fitted
+and a best set of galaxy models have been determined, three sets of
+radial apertures are measured.  In the first set, the fluxes in the
+radial apertures are measured using the raw stack images.  The centers
+of the apertures for each object across the 5 filters are fixed so
+that the pixels represent the equivalent portions of the same galaxy
+for all 5 filters.  In this analysis, the best model for each object
+is subtracted from the image pixels for all objects excluding the
+object in consideration.  The 'best model' is determined based on the
+minimum $\chi^2$ value for the model fits.
+% \note{more discussion of the selection of the best model}.
+In addition to the raw radial apertures, the stack images are each
+convolved with a circular Gaussian with $\sigma$ chosen to yield an
+image with a typical FWHM of 6\arcsec.  The full set of radial
+apertures are again measured on these convolved images.  Again, the
+best object models are subtracted from the image for objects not being
+measured.  This subtraction includes the convolution to smooth the
+model to the effective FWHM of the convolved image.  The entire
+procedure is then repeated with a target FWHM of 8\arcsec.
+% \note{is the first convolution done with the Alard-Lupton technique?}
+\acknowledgments
+The Pan-STARRS1 Surveys (PS1) have been made possible through
+contributions of the Institute for Astronomy, the University of
+Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its
+participating institutes, the Max Planck Institute for Astronomy,
+Heidelberg and the Max Planck Institute for Extraterrestrial Physics,
+Garching, The Johns Hopkins University, Durham University, the
+University of Edinburgh, Queen's University Belfast, the
+Harvard-Smithsonian Center for Astrophysics, the Las Cumbres
+Observatory Global Telescope Network Incorporated, the National
+Central University of Taiwan, the Space Telescope Science Institute,
+the National Aeronautics and Space Administration under Grant
+No. NNX08AR22G issued through the Planetary Science Division of the
+NASA Science Mission Directorate, the National Science Foundation
+under Grant No. AST-1238877, the University of Maryland, and Eotvos
+Lorand University (ELTE) and the Los Alamos National Laboratory.
+\bibliographystyle{apj}
+% \bibliography{lib}{}
+\input{analysis.bbl}
+\end{document}
+\subsection{Forced Photometry : PSFs}
+\subsection{Forced Photometry : galaxies}
 \subsection{Output Options}
 \note{need to discuss tests}
 \note{need to discuss failings and holes}
+% \note{need to discuss tests}
+% \note{need to discuss failings and holes}
 \section{Alternative Scenarios}
 …
 \subsection{Trailed Sources}
-\subsection{Fixed / Known-position Sources}
 \subsection{Difference Images}
 The variance map for a difference image must be generated from the two
 images use to construct the difference.  Otherwise, the low sky level
+images used to construct the difference.  Otherwise, the low sky level
 will automatically result in inconsistent interpretation of the variance.
 …
 \section{Sample Tests}
+\end{document}
+\begin{verbatim}
+Configuration variables affecting the peak detection process:
+PEAKS_SMOOTH_SIGMA [2.5]   : Gaussian sigma of smoothing kernel, in pixels.
+PEAKS_SMOOTH_NSIGMA [2.0]  : Gaussian smoothing kernel window size in sigmas.
+PEAKS_NSIGMA_LIMIT [20.0]  : Detection limit on first pass (sigmas).
+PEAKS_NSIGMA_LIMIT_2 [5.0] : Detection limit on faint detection pass (sigmas).
+\end{verbatim}
+Recipe parameters which affect the PSF-convolved galaxy model fitting
+process:
+\begin{verbatim}
+EXT_FIT_NSIGMA_CONV [9] : number of sigma
+EXT_FIT_ITER
+EXT_FIT_MIN_TOL
+EXT_FIT_MAX_TOL
+LMM_FIT_CHISQ_CONVERGENCE
+LMM_FIT_GAIN_FACTOR_MODE
+\end{verbatim}
+Figures Needed for this document:
+* aperture - PSF model example
+* map of the sky with galaxy region illustrated?
+* plots showing the quality of the data?
+Tables needed:
+* table of models?
+Work still needed:
+* Tables
+* refereces for other programs
+* authors
+* PSF residual map
+* section 3.5.3 Model applied to detected objects needs to be reviewed
+* read for english & phrasing
+* convert XXX
+* background model description (see waters)
+* reduce / remove design goal wording?
+* too many bullets?
+* reduce coding description?
+* put engineering docs (psLib, psModules) on public website

Note: See TracChangeset for help on using the changeset viewer.

Context Navigation

Changeset 39920 for branches/czw_branch/20160809/doc/release.2015/ps1.analysis/analysis.tex

Legend:

branches/czw_branch/20160809/doc/release.2015/ps1.analysis

branches/czw_branch/20160809/doc/release.2015/ps1.analysis/analysis.tex

Download in other formats: