\documentclass[iop,floatfix]{emulateapj}
% \documentclass[10pt,preprint]{aastex}
% \pdfoutput=1

\RequirePackage{graphicx}
\RequirePackage{color}
\RequirePackage{code}
\input{astro.sty}

\usepackage[T1]{fontenc}% (2) specify encoding

% online version may use color, but print version needs b/w
\def\plotmode{col}
%\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{png}
\def\plotext{pdf}

%\def\picdir{/home/eugene/chipresid.20140404}
\def\picdir{pics}
%\def\picdir{.}

% Pick a terse version of the title here;
\shorttitle{Pixel Analysis in PS1}
\shortauthors{E.A. Magnier et al}
\begin{document}
\title{Pan-STARRS Pixel Analysis : Source Detection and Characterization}

% this is a crude trick to get the order of affiliations right.  These
% names are used in the affiliations below.  The user needs to (1) set
% the order and numbers to have the correct sequence in the author
% 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\MPE{4}
\def\CfA{5}

% This example has a first author from UH:
\author{
Eugene A. Magnier,\altaffilmark{\IfA}
% 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.~W. Draper,\altaffilmark{\DUR}
D. Farrow,\altaffilmark{\DUR,\MPE}
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}
% T. Grav,\altaffilmark{\IfA}
% J. N. Heasley,\altaffilmark{\IfA}
% G. A. Luppino,\altaffilmark{\IfA}
% R. H. Lupton,\altaffilmark{\Princeton}
% E. A. Magnier,\altaffilmark{\IfA}
% D. G. Monet,\altaffilmark{\USNO}
% J.~S. Morgan,\altaffilmark{\IfA}
% P. M. Onaka,\altaffilmark{\IfA}
% J.~L. Tonry, \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{\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}
% \altaffiltext{\MPIA}{Max Planck Institute for Astronomy, K\"onigstuhl 17, D-69117 Heidelberg, Germany}
\altaffiltext{\MPE}{Max-Planck-Institut f\"ur extraterrestrische Physik, Giessenbachstrasse 1, 85748 Garching, Germany}
\begin{abstract}

Over 3 billion astronomical sources 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
sources have been automatically detected and characterized by the
Pan-STARRS Image Processing Pipeline photometry software,
\ippprog{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 sources by \ippprog{psphot} in general as well as for the
specific case of the 3rd processing version used for the first
\textmod{two public releases} of the Pan-STARRS $3\pi$ survey data.
\end{abstract}

% insert additional keywords as appropriate:
\keywords{methods: data analysis -- Surveys:\PSONE -- techniques: image processing -- techniques: photometric}

\section{Introduction}
\label{sec:intro}

% \begin{verbatim}
% here is a list of things to do:
%  * explain use of covariance
%  * add example for sky model
%  * define more tests and generate examples
%   * simulation example of background subtraction
%     at different densities
%   * real example of oversubtracted galaxy
% \end{verbatim}

The 1.8m Pan-STARRS\,1 telescope (PS1) is located on the summit of Haleakala
on the Hawaiian island of Maui.  The wide-field optical design of the
telescope \citep{2004SPIE.5489..667H} produces a 3.3 degree field of view with
low distortion and minimal vignetting even at the edges of the
illuminated region.  The optics and natural seeing combine to yield
good image quality: 75\% of the images have full-width half-max values
less than (1.51, 1.39, 1.34, 1.27, 1.21) arcseconds for (\grizy), with
a floor of $\sim 0.7$ arcseconds.

The \PSONE\ camera \citep{2009amos.confE..40T}, known as GPC1, consists of a
mosaic of 60 back-illuminated CCDs manufactured by Lincoln Laboratory.
The CCDs each consist of an $8\times8$ grid of $590 \times 598$
pixel readout regions, yielding an effective $4846 \times 4868$
detector.  Initial performance assessments are presented in
\cite{2008SPIE.7014E..0DO}.  Routine observations are conducted remotely from the
Advanced Technology Research Center in Kula, the main facility of the
University of Hawaii's Institute for Astronomy (IfA) operations on Maui.
The Pan-STARRS1 filters and photometric system have already been
described in detail in \cite{2012ApJ...750...99T}.

For nearly 4 years, from 2010 May through 2014 March, this telescope
was used to perform a collection of astronomical surveys under the
aegis of the Pan-STARRS Science Consortium.  The majority of the time
(56\%) was spent on surveying the $\frac{3}{4}$ of the sky north of
$-30$ Declination with \grizy\ filters in the so-called $3\pi$ Survey.
Another $\sim 25\%$ of the time was concentrated on repeated deep
observations of 10 specific fields in the Medium-Deep Survey.  The
rest of the time was used for several other surveys, including a
search for potentially hazardous asteroids in our solar system.  The
details of the telescope, surveys, and resulting science publications
are described by \cite{chambers2017}.

Since 2014 March, PS1 has been re-dedicated to a mission of searching
for hazardous asteroids, funded by the NASA NEO Program. Additional
partners collaborate with the Pan-STARRS team to harvest the transient
sources such supernovae and graviational wave counterparts.  A second
Pan-STARRS telescope \citep[PS2][Chambers et al 2020 in prep]{chambers2017},
generally matching the PS1 design \citep{2012SPIE.8444E..0HM} has since been
constructed and has been producing science results since early 2018.

Pan-STARRS produced its first large-scale public data release, Data
Release 1 (DR1) on 16 December 2016.  DR1 contains the results of the
third full reduction of the Pan-STARRS $3\pi$ Surveyo archival data,
identified as PV3.  Previous reductions \citep[PV0, PV1, PV2;
  see][]{magnier2017.datasystem} were used internally for pipeline
optimization and the development of the initial photometric and
astrometric reference catalog \citep{magnier2017.calibration}.  The
products from these reductions were not publicly released, but have
been used to produce a wide range of scientific papers from the
Pan-STARRS 1 Science Consortium members \citep{chambers2017}.  DR1
contained only average information resulting from the many individual
images obtained by the $3\pi$ Survey observations.  A second data
release, DR2, was made available 28 January 2019.  DR2 provides
measurements from all of the individual exposures, and includes an
improved \textmod{astrometric calibration as well as improvements to the
  photometric calibration of the stack and `forced warp' measurements
from} the PV3 processing of that dataset.

\textadd{The Pan-STARRS public data releases are hosted by the {\em
    Barbara A. Mikulski Archive for Space Telescopes} (MAST) at the
  Space Telescope Science Institute (STScI).  MAST provides access to
  the image data products and a hierachical database of measurements
  using a system developed specifically for the Pan-STARRS dataset.
  Development of this database systems was the product of a
  collaboration between the Pan-STARRS Project and Alex Szalay's
  database development group at The Johns Hopkins University (JHU)
  \citep{2008AIPC.1082..352H}.  The resulting system, called the
        {\em Published Science Products Subsystem}, or PSPS
        \citep{Heasley2006}, was initially used within the Pan-STARRS
        Science Consortium for large-scale data access.  A duplicate
        PSPS installation was created at MAST for the DR1 and DR2
        public releases.}

This is the fourth in a series of seven papers describing the
Pan-STARRS1 Surveys, the data reduction techniques 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 model fitting, and the techniques for ``forced''
photometry measurements.  \textadd{The same analysis software, called \ippprog{psphot}, is used
  for individual images, image stacks, and difference images.}
The software described here was used with a
single consistent set of parameters for the complete PV3 analysis,
used for both DR1 and DR2.  \textadd{The software was also used for the
analysis of the Medium Deep Survey data, though with a different
software version and some modifications of
the analysis parameters to better suite the longer exposures.  This
program as well as the rest of the Pan-STARRS Image Processing
Pipeline (IPP) software suite is available for download from \url{http://ipp.ifa.hawaii.edu}}.

% \note{Generate a tarball of just the programs (skip certain directories)}

%Chambers et al. 2017 (Paper I)
%The Pan-STARRS\,1 Surveys
\citet[][Paper I]{chambers2017}
provide 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]{magnier2017.datasystem}
describe how the various data processing stages are organized and implemented
in the \textmod{Image 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} describe 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]{magnier2017.analysis}
%describe 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]{magnier2017.calibration}
describe 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}
describe  the details of the resulting catalog data and its
organization in the Pan-STARRS database system, PSPS. 

%Huber et al. 2017 (Paper VII)
\citet[][Paper VII]{huber2017} describe 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 Releases 1 or 2 and
will be made available in a future data release.

%%{\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.}}

\textadd{In this article, we use the following type-faces to distinguish
different concepts:}
\begin{itemize}
\item \ippstage{Small caps} for the analysis stages.
\item \ippdbtable{Italics} for database tables and columns.
\item \ippprog{Fixed-width} font for program names, variables, and
  miscellaneous constants.
\end{itemize}

\textadd{
The latter catagory includes a number of configuration parameters used
to define the \ippprog{psphot} analysis.  In those cases, unless the
values used for the PV3 analysis are explicitly discussed, we include
the PV3 value immediately after the configuration variable name in parenthesis.}

\section{Background}

The photometric and astrometric precision goals for the Pan-STARRS\,1
surveys were quite stringent.  The astrometric goals were relative
astrometric accuracy of 10 milliarcseconds and absolute astrometric
accuracy of 100 milliarcseconds with respect to the ICRS reference
stars.  For photometry, the goal was 10 millimagnitudes accuracy
within the internal photometric system across the sky, though the tie
to an absolute standard was not required to meet this standard.

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 sources in the image data stream in a
timely fashion are a significant challenge.

In order to achieve these ambitious goals, the source detection,
classification, and measurement process must be both precise and
efficient.  Not only is it necessary to make a careful measurement of
the flux of individual sources, it is also critical to characterize
the image point spread function (PSF), and its variations across the
field and from image to image.  Since comparisons between images must
be reliable, the measurements must be stable for both photometry and
astrometry.  

A variety of astronomical software packages perform the basic source
detection, measurement, and classification tasks needed by the
Pan-STARRS IPP.  Each of these programs have their own advantages and
disadvantages.  Below we discuss some of the most widely used of these
other packages, highlighting the features of the programs which are
particularly desirable, and noting aspects of the programs which are
problematic for the IPP.

\begin{itemize}

\item DoPhot : analytical fitted model with aperture corrections.
  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 \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? \citep{1987PASP...99..191S}.

\item Sextractor : pure aperture measurement with rudimentary source
  subtraction.  pro: fast, widely used, easy to automate.  con: poor
  source 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-source 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 \citep{2002AJ....124..266P}.

\item SDSS phot : con: tightly integrated into the SDSS software
  environment \citep{2001ASPC..238..269L}.

\end{itemize}

When the IPP development was starting, the existing photometry
packages either did not meet the accuracy requirements or 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 pipeline.  A new photometry analysis package was
developed using lessons learned from the existing photometry systems.
In the process, the source analysis software was written using the
data analysis C-code library written for the IPP, \code{psLib}
\citep{magnier2017.datasystem}.  Components of the photometry code
were integrated into the IPP's mid-level astronomy data analysis
toolkit called \code{psModules} \citep{magnier2017.datasystem}.  The
resulting software, `\ippprog{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

Several variants of \ippprog{psphot} have been used in the PS1 PV3
analysis.  The main variant of \ippprog{psphot} operates on a single
image, or a group of related images representing the data read from
\textmod{the multiple chips of a mosaic 
camera from} 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 specified by the configuration system, or a variance
image may be supplied.  A mask may also be supplied to mark good, bad,
and suspect pixels.  This variant of \ippprog{psphot} can be called as a
stand-alone program, also called \ippprog{psphot}.  In standard IPP
operations, this variant is used as a library call within the analysis
program \ippprog{ppImage} during the \ippstage{chip} analysis stage.

\textadd{In the standard IPP analysis, the initial stage of processing
  is performed in parallel on each of the individual CCDs in the
  camera.  This so-called \ippstage{chip} stage analysis includes the
  detrending of CCD image as well as the detection and analysis of
  sources in the image using the basic version of \ippprog{psphot}.
  The next stage of the analysis, the \ippstage{camera} stage,
  consists of photometric and astrometric calibration.}

  \textadd{After the calibrations are available, the detrended CCD images from
  an entire exposure are geometrically transformed to a common pixel
  grid in the \ippstage{warp} stage of the pipeline.  The resulting
  warped images are generated on a pre-define tessellation of the sky
  which starts with projection centers spaced roughly
  4\degrees\ across the sky.  Around each of these projection centers,
  a large regular pixel grid is defined, and then subdivided along
  pixel boundaries into smaller units which are well-matched to the
  memory footprint of our processing computesr.  These smaller images,
  called `skycells' are defined with 1 arcminute of overlap with their
  neighbors to that any modest-sized object can be analysed entirely
  on a single pixel grid.  Note that the term skycell is used to
  describe the particular subdivision of the sky.  A typical exposure
  from the GPC1 camera generates warp images on roughly 70 skycells.
  We refer to the specific warped images from an exposure as `warps'.}

  \textadd{Multiple warps for the same skycell are combined together in the
  \ippstage{stack} stage of the IPP by co-adding the flux to generate
  a deep `stack' image.  Alternatively, one warp may be subtracted
  from another warp of the same skycell, or a stack image may be
  subtracted from a warp image, or indeed from another stack.  These
  subtraction operations are used to detect moving and transient
  objects within the IPP.  Different variants of \ippprog{psphot} are
  used for the source detection and analysis for each of these
  different analysis stages.}

The variant called \ippprog{psphotStack} accepts a set of images, each
representing the same patch of sky \textadd{(with pixels aligned)} in
a different \textmod{filter.  This version was used in the IPP for the
  analysis of the deep ``stacks'' produced by the IPP \ippstage{stack}
  stage.  Nominally, the full $grizy$ filter set was used 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
\ippprog{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 includes an option to
convolve the set of images to a single, common PSF size across the
filters for the purpose of fixed aperture photometry.

Another variant of \ippprog{psphot} used in the PV3 analysis is called
\ippprog{psphotFullForce}.  In this variant, a set of images all representing the
same \textadd{co-aligned} pixels are processed together, with the positions of sources to
be analyzed loaded from a supplied file.  In this variant of 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{\nocode{psphot} Design Goals}

% \subsection{Astronomy Measurement Goals}

\textadd{The top-level design goals of \ippprog{psphot} are to detect and
determine the instrumental positions and fluxes of astronomical
sources in the images.  For extended sources, the goals also include
the measurement of a variety of morphological information, including
galaxy model parameters and non-parametric measurements of the sizes
and profiles of the galaxies to aid in classification and for
weak-lensing analysis.  For trailed asteroids, the goal also includes
the measurement of the length and direction of the trail.}

\textmod{Beyond these basic elements, \ippprog{psphot} has a number of
  design goals} which we believe will help to make it usable in a wide
range of circumstances.  The critical astronomy-driven measurement
goals of the Pan-STARRS project, which drive the design of
\ippprog{psphot}, are the photometric accuracy goal (10
millimagnitudes) and the \textadd{relative} astrometric accuracy goal
(10 milliarcseconds) \textadd{for bright stars for which the photon
  shot-noise is small compared to the systematic errors.}

For \ippprog{psphot}, the photometry accuracy goal implies that the
measured photometry of stellar sources must be substantially better
than this 10 mmag goal since the photometry error per image is
combined with an error in the flat-field calibration and an error in
measuring the atmospheric effects.  We have set a goal for
\ippprog{psphot} of 3 mmag photometric consistency for bright stars
between pairs of images obtained in photometric conditions at the same
pointing, ie to remove sensitivity to flat-field errors.  This goal
splits the difference between the three main contributors and still
allows some leeway.  This goal must be met for well-sampled
images and images with only modest undersampling.

The relative astrometric calibration depends on the consistency of the
individual measurements.  The measurements from \ippprog{psphot} must
be sufficiently representative of the true source position to enable
astrometric calibration at the 10 mas level.  The error in the
individual measurements will be folded together with the errors
introduced by the optical system, the effects of seeing, and by the
available reference catalogs.  We have set a goal for \ippprog{psphot}
of 5 mas consistency between the true source postion and the measured
position given reasonable PSF variations under simulations.  This
level must be reached for images with 250 mas pixels, implying
\ippprog{psphot} must introduce measurement errors less than 1/50th of
a 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 mas.

% \subsection{Software System Goals}

The design goals for \ippprog{psphot} are chosen to make the program flexible,
general, and able to meet the unknown usage cases future projects may
require:

\begin{itemize}
\item {\bf Flexible PSF model} Different image sources require
  different ways of representing the PSF.  Ideally, both analytical
  and pixel-based versions should be possible.

\item {\bf PSF spatial variation} Most images result in some spatial
  PSF variations at a certain level.  The PSF representation should
  naturally incorporate 2-D variations.

\item {\bf Flexible non-PSF models} \ippprog{psphot} must be able to represent
  PSF-like sources as well as non-PSF sources (e.g., galaxies).  It
  must be easy to add new source models as interesting representations
  of sources are invented.

\item {\bf Clean code base} \ippprog{psphot} should incorporate a high-degree of
  abstraction and encapsulation so that changes to the code structure
  can be performed without pulling the code apart and starting from scratch.

\item {\bf PSF validity tests} \ippprog{psphot} should include the ability to
  choose different types of PSF models for different situations, or to
  provide the user with methods for assessing the different PSF models.

\item {\bf Careful systematic corrections} \ippprog{psphot} must carefully
  measure and correct for the photometric and astrometric trends
  introduced by using analytical PSF models.

\item {\bf User Configurable} \ippprog{psphot} should allow users to change the
  options easily and to allow different approaches to the analysis.

\end{itemize}

\textadd{The success of the \ippprog{psphot} implementation in meeting
  the photometry and astrometry design requirements is demonstrated by
  the achieved accuracy for the Pan-STARRS $3\pi$ Survey data.  For a
  survey like the Pan-STARRS\,1 $3\pi$ survey to achieve photometry
  and astrometry accuracy at the level of our goals, not only must the
  measurement of the astronomical detections be precise, but it is
  necessary for the detrending (instrumental signature remove) and
  calibration processes to correct for a wide variety of systematic
  effects and it is also necessary for the observations to be
  performed in such a way that the data can be calibrated well.  These
  others aspects of the process are discussed in detail elsewhere
  (Papers I, III, V).  In the end, the goals were largely achieved for
  the Pan-STARRS\,1 $3\pi$ survey. As reported in Paper V, the
  resulting photometric system is consistent across the sky to between
  7 and 12.4 millimagnitudes depending on the filter.  The systematic
  error floor for individual photometry measurements is $(\sigma_g,
  \sigma_r, \sigma_i, \sigma_z, \sigma_y) = (14, 14, 15, 15, 18)$
  millimagnitudes.  The bright-star systematic error floor for
  individual astrometric measurements is 16 milliarcseconds and the
  Pan-STARRS Data Release 2 (DR2) astrometric system is tied to the
  Gaia DR1 \citep{2016AA...595A...4L} coordinate frame with a
  systematic uncertainty of $\sim 5$ milliarcseconds. }

\section{Basic Analysis}

\subsection{Overview}

The basic \ippprog{psphot} analysis is divided into several major
stages, as listed below.

\begin{enumerate}
\item {\bf Image Preparation} Load data, characterize the image
  background, load or construct variance and mask images.

\item {\bf Initial Source Detection} Smooth, find peaks, measure basic
  properties with focus on the point sources to measure the PSF.

\item {\bf PSF Determination} Select PSF candidates, perform model
  fits, build PSF model from fits, select best PSF model class.

\item {\bf Bright Source Analysis} Fit sources with PSFs, determine
  PSF validity, subtract PSF-like sources, fit non-PSF model(s),
  select best model class, subtract model.

\item {\bf Faint Source Analysis} Detect low-level sources, measure
  properties (aperture or PSF)

\item {\bf Aperture corrections} Measure the curve-of-growth, spatial
  aperture variations, and background-error corrections.  

\item {\bf Output} Write out sources in selected format, write out
  difference image, variance image, etc, as selected.  
\end{enumerate}

In addition to this basic sequence, additional analysis steps may be
performed.  An ``extended source'' analysis mode is available to
measure photometry and morphology of galaxies and other resolved
sources.  Forced photometry may be performed for both point-like and
extended sources.  A special mode is available for the photometry of
sources detected in difference images.  These different modes are
discussed in their own sections below.

Table~\ref{tab:measurements} lists the types of
analyses performed by \ippprog{psphot}, specifying which of the
\ippprog{psphot} usage cases performs the given analysis.  The table
also provides a reference to the section of this paper in which the
analysis is described.  Not all analyses are relevant to all sources
in all images.  The table identifies thoses cases where the analyses
are applied to only a subset of all sources.  

\ippprog{psphot} is highly configurable.  Users may choose via the configuration
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.

\textadd{Ultimately, all measurements of individual astronomical
  sources from \ippprog{psphot} are reported in one of the tables in
  the PSPS database.  As discussed in detail in Paper VI, measurements
  from individual exposures are available from the
  \ippdbtable{Detection} table.  Measurements of objects in the
  stacked images are stored in one of several \ippdbtable{Stack...}
  tables while the `forced' measurements from individual warp images
  are stored in tables beginning with \ippdbtable{ForcedWarp...}.}

\begin{table*}
\caption{\label{tab:measurements} Measurements performed by
  \nocode{psphot}, and whether performed in each of the 4 IPP analysis
  stages.  The analysis is described in this article in the listed Sections. } % \vspace{-0.5cm}
\begin{center}
\footnotesize
\begin{tabular}{lccccll}
\hline
\hline
{\bf Measurement} & {\sc \bf CHIP} & {\sc \bf STACK} & {\sc \bf FORCED
  WARP} & {\sc \bf DIFF} & {\bf Section} & {\bf Details} \\
\hline
  Background Subtraction     & Y & Y & Y & N$^1$ & \ref{sec:image.preparation}      & N/A \\
  Peaks                      & Y & Y & N & Y     & \ref{sec:peaks}                  & All detections \\
  Footprints                 & Y & Y & N & Y     & \ref{sec:footprints}             & All detections \\
  Moments                    & Y & Y & Y & Y     & \ref{sec:moments}                & All detections \\
  PSF Model                  & Y & Y & Y & N$^2$ & \ref{sec:PSF.Model}              & Selected bright stars \\
  Bright Star Profile        & Y & Y & N & Y     & \ref{sec:very.bright.star}       & Saturated Stars \\
  Radial Profiles v1         & Y & Y & N & Y     & \ref{sec:radial.profile}         & All detections \\
  Kron Fluxes                & Y & Y & Y & Y     & \ref{sec:kron.mags}              & All detections \\
  Source-Size Tests          & Y & Y & N & Y     & \ref{sec:source.size}            & All detections \\
  Non-Linear PSF Fits        & Y & Y & N & N     & \ref{sec:nonlinear.psf.model}    & $S/N > 20$ \\
  Unconvolved Galaxy Model   & Y & Y & N & N     & \ref{sec:nonlinear.galaxy.model} & $S/N > 20$, extended \\
  Unconvolved Streak Model   & N & N & N & Y     & \ref{sec:nonlinear.galaxy.model} & $S/N > 20$, extended \\
  Linear PSF Fits            & Y & Y & Y & Y     & \ref{sec:faint.psf.model}        & All detections \\
  Radial Profiles v2         & Y & Y & N & Y     & \ref{sec:radial.profile.v2}      & Gal. Latitude Cut \\
  Petrosian Fluxes           & N & Y & Y & N     & \ref{sec:petrosian}              & Gal. Latitude Cut \\
  Convolved Galaxy Models    & N & Y & N & N     & \ref{sec:galaxy.conv.fit}        & Gal. Latitude Cut, mag cut \\
  Fixed Aperture Photometry  & N & Y & Y & N     & \ref{sec:fixed.aperture.photom}  & All detections \\
  Convolved, Fixed Apertures & N & Y & N & N     & \ref{sec:fixed.aperture.photom}  & All detections \\
  Aperture Corrections       & Y & Y & Y & N     & \ref{sec:aperture.correction}    & All detections \\
  Forced PSF Fluxes          & N & N & Y & N     & \ref{sec:psf.forced.fit}         & All detections \\
  Forced Galaxy Models       & N & N & Y & N     & \ref{sec:galaxy.forced.fit}      & Requires stack galaxy models \\
  Lensing Parameters         & N & Y & Y & N     & \ref{sec:lensing.params}         & All detections \\
\hline
\multicolumn{5}{l}{$^1$ Background subtraction is performed by {\tt ppSub} before calling {\tt psphot}} \\
\multicolumn{5}{l}{$^2$ PSF modeling is performed by {\tt ppSub} on the input warps before calling {\tt psphot}} \\
\end{tabular}
\end{center}
\end{table*}

\subsection{Informational and Warning Bit Flags}

During the \ippprog{psphot} analysis, there are a wide variety of
conditions which are identified by the analysis software.  As part of
the output data for each detected source, two fields are provided
which encode these conditions as bit values in the two 32-bit
integers.  The following two tables list the individual bit values in
these two fields.  These informational and warning bits are described
in more detail later in this article.

%
Table~\ref{tab:det_flag_values} lists the flags recorded in the output
field \ippmisc{FLAGS}.  When data from \ippprog{psphot} is loaded into
a DVO database \citep{magnier2017.calibration}, these values are
stored in the field \ippdbtable{Measure.photFlags} and exposed in the public
database \citep[PSPS][]{flewelling2017} in the fields
\ippdbtable{Detection.infoFlag}, \ippdbtable{StackObjectThin.XinfoFlag} (where
\ippdbtable{X} is one of {$grizy$}), and
\ippdbtable{ForcedWarpMeasurement.FinfoFlag}.
%
Table~\ref{tab:det_flag2_values} lists the flags recorded in the
output field \ippmisc{FLAGS2}.  When data from \ippprog{psphot} is
loaded into a DVO database \citep{magnier2017.calibration}, these
values are stored in the field \ippdbtable{Measure.photFlags2}, and they are exposed in PSPS in the fields
\ippdbtable{Detection.infoFlag2}, \ippdbtable{StackObjectThin.XinfoFlag2} (where
\ippdbtable{X} is one of {$grizy$}), and
\ippdbtable{ForcedWarpMeasurement.FinfoFlag2}.

\begin{table*}
\begin{center}
\caption{\label{tab:det_flag_values}
Detection Flag
Values \#1 reported by \texttt{psphot}. These are saved in output catalogs as the field
\texttt{FLAGS}, in the DVO database as
\textit{Measure.photFlags}, and in the public database as
\textit{Detection.infoFlag},
\textit{StackObjectThin.XinfoFlag} (where \textit{X} is one
of {$grizy$}), and \textit{ForcedWarpMeasurement.FinfoFlag}.}
\footnotesize
\begin{tabular}{lrl}
\hline
\hline
{\bf Flag Name} & {\bf Flag Value} & {\bf Description} \\
\hline
 PM\_SOURCE\_MODE\_PSFMODEL            & 0x00000001 & Source fitted with a PSF model (linear or non-linear) \\
 PM\_SOURCE\_MODE\_EXTMODEL            & 0x00000002 & Source fitted with an extended-source model \\
 PM\_SOURCE\_MODE\_FITTED              & 0x00000004 & Source fitted with non-linear model (PSF or EXT; good or bad) \\
 PM\_SOURCE\_MODE\_FAIL                & 0x00000008 & Fit (non-linear) failed (non-converge, off-edge, run to zero) \\
 PM\_SOURCE\_MODE\_POOR                & 0x00000010 & Fit succeeds, but low-S/N or high chi-square \\
 PM\_SOURCE\_MODE\_PAIR                & 0x00000020 & Source fitted with a double PSF \\
 PM\_SOURCE\_MODE\_PSFSTAR             & 0x00000040 & Source used to define PSF model \\
 PM\_SOURCE\_MODE\_SATSTAR             & 0x00000080 & Source model peak is above saturation \\
 PM\_SOURCE\_MODE\_BLEND               & 0x00000100 & Source is a blend with other sources$^1$ \\
 PM\_SOURCE\_MODE\_EXTERNAL            & 0x00000200 & Source based on supplied input position \\
 PM\_SOURCE\_MODE\_BADPSF              & 0x00000400 & Failed to get good estimate of object's PSF \\
 PM\_SOURCE\_MODE\_DEFECT              & 0x00000800 & Source is thought to be a defect \\
 PM\_SOURCE\_MODE\_SATURATED           & 0x00001000 & Source is thought to be saturated pixels (bleed trail) \\
 PM\_SOURCE\_MODE\_CR\_LIMIT           & 0x00002000 & Source has crNsigma above limit \\
 PM\_SOURCE\_MODE\_EXT\_LIMIT          & 0x00004000 & Source has extNsigma above limit \\
 PM\_SOURCE\_MODE\_MOMENTS\_FAILURE    & 0x00008000 & could not measure the moments \\
 PM\_SOURCE\_MODE\_SKY\_FAILURE        & 0x00010000 & could not measure the local sky \\
 PM\_SOURCE\_MODE\_SKYVAR\_FAILURE     & 0x00020000 & could not measure the local sky variance \\
 PM\_SOURCE\_MODE\_BELOW\_MOMENTS\_SN  & 0x00040000 & moments not measured due to low S/N.$^1$ \\
 PM\_SOURCE\_MODE\_BIG\_RADIUS         & 0x00100000 & poor moments for small radius, try large radius \\
 PM\_SOURCE\_MODE\_AP\_MAGS            & 0x00200000 & source has an aperture magnitude \\
 PM\_SOURCE\_MODE\_BLEND\_FIT          & 0x00400000 & source was fitted as a blend \\
 PM\_SOURCE\_MODE\_EXTENDED\_FIT       & 0x00800000 & full extended fit was used \\
 PM\_SOURCE\_MODE\_EXTENDED\_STATS     & 0x01000000 & extended aperture stats calculated \\
 PM\_SOURCE\_MODE\_LINEAR\_FIT         & 0x02000000 & source fitted with the linear fit \\
 PM\_SOURCE\_MODE\_NONLINEAR\_FIT      & 0x04000000 & source fitted with the non-linear fit \\
 PM\_SOURCE\_MODE\_RADIAL\_FLUX        & 0x08000000 & radial flux measurements calculated \\
 PM\_SOURCE\_MODE\_SIZE\_SKIPPED       & 0x10000000 & size could not be determined$^1$ \\
 PM\_SOURCE\_MODE\_ON\_SPIKE           & 0x20000000 & peak lands on diffraction spike \\
 PM\_SOURCE\_MODE\_ON\_GHOST           & 0x40000000 & peak lands on ghost or glint \\
 PM\_SOURCE\_MODE\_OFF\_CHIP           & 0x80000000 & peak lands off edge of chip \\
\hline
\multicolumn{3}{l}{$^1$ Not used for DR1 or DR2.} \\
\hline
\end{tabular}
\end{center}
\end{table*}

\begin{table*}
\caption{\label{tab:det_flag2_values}
Detection Flag Values \#2 reported by \nocode{psphot}.
These are saved in output catalogs as the field
\texttt{FLAGS2}, in the DVO database as
\textit{Measure.photFlags2}, and in the public database as
\textit{Detection.infoFlag2},
\textit{StackObjectThin.XinfoFlag2} (where \textit{X} is one
of $grizy$), and \textit{ForcedWarpMeasurement.FinfoFlag2}.
}
\begin{center}
\footnotesize
\begin{tabular}{lrl}
\hline
\hline
{\bf Flag Name} & {\bf Flag Value} & {\bf Description} \\
\hline
 PM\_SOURCE\_MODE2\_DIFF\_WITH\_SINGLE      & 0x00000001 & diff source matched to a single positive detection \\
 PM\_SOURCE\_MODE2\_DIFF\_WITH\_DOUBLE      & 0x00000002 & diff source matched to positive detections in both images \\
 PM\_SOURCE\_MODE2\_MATCHED          	    & 0x00000004 & source generated based on another image \\
 PM\_SOURCE\_MODE2\_ON\_SPIKE               & 0x00000008 & $> 25\%$ of (PSF-weighted) pixels land on diffraction spike \\
 PM\_SOURCE\_MODE2\_ON\_STARCORE            & 0x00000010 & $> 25\%$ of (PSF-weighted) pixels land on starcore \\
 PM\_SOURCE\_MODE2\_ON\_BURNTOOL            & 0x00000020 & $> 25\%$ of (PSF-weighted) pixels land on burntool \\
 PM\_SOURCE\_MODE2\_ON\_CONVPOOR            & 0x00000040 & $> 25\%$ of (PSF-weighted) pixels land on convpoor \\
 PM\_SOURCE\_MODE2\_PASS1\_SRC              & 0x00000080 & source detected in first pass analysis \\
 PM\_SOURCE\_MODE2\_HAS\_BRIGHTER\_NEIGHBOR & 0x00000100 & peak is not the brightest in its footprint \\
 PM\_SOURCE\_MODE2\_BRIGHT\_NEIGHBOR\_1     & 0x00000200 & $flux_{\rm n} / (r^2 flux_{\rm p}) > 1$ \\
 PM\_SOURCE\_MODE2\_BRIGHT\_NEIGHBOR\_10    & 0x00000400 & $flux_{\rm n} / (r^2 flux_{\rm p}) > 10$ \\
 PM\_SOURCE\_MODE2\_DIFF\_SELF\_MATCH       & 0x00000800 & positive detection match is probably this source \\
 PM\_SOURCE\_MODE2\_SATSTAR\_PROFILE        & 0x00001000 & saturated source is modeled with a radial profile \\
 PM\_SOURCE\_MODE2\_ECONTOUR\_FEW\_PTS      & 0x00002000 & too few points to measure the elliptical contour \\
 PM\_SOURCE\_MODE2\_RADBIN\_NAN\_CENTER     & 0x00004000 & radial bins failed with too many NaN center bin \\
 PM\_SOURCE\_MODE2\_PETRO\_NAN\_CENTER      & 0x00008000 & petrosian radial bins failed with too many NaN center bin$^1$ \\
 PM\_SOURCE\_MODE2\_PETRO\_NO\_PROFILE      & 0x00010000 & petrosian not build because radial bins missing \\
 PM\_SOURCE\_MODE2\_PETRO\_INSIG\_RATIO     & 0x00020000 & insignificant measurement of petrosian ratio \\
 PM\_SOURCE\_MODE2\_PETRO\_RATIO\_ZEROBIN   & 0x00040000 & petrosian ratio in the 0th bin (likely bad) \\
 PM\_SOURCE\_MODE2\_EXT\_FITS\_RUN          & 0x00080000 & we attempted to run extended fits on this source \\
 PM\_SOURCE\_MODE2\_EXT\_FITS\_FAIL         & 0x00100000 & at least one of the model fits failed \\
 PM\_SOURCE\_MODE2\_EXT\_FITS\_RETRY        & 0x00200000 & trailed asteroid model fit was re-tried with new window \\
 PM\_SOURCE\_MODE2\_EXT\_FITS\_NONE         & 0x00400000 & ALL of the model fits failed \\
\hline
\multicolumn{3}{l}{$^1$ Not used for DR1 or DR2.} \\
\hline
\end{tabular}
\end{center}
\end{table*}

\subsection{Image Preparation}
\label{sec:image.preparation}

The first step is to prepare the image for detection of the
astronomical sources.  We need three separate images: the measured
flux (signal image), the corresponding variance image, and a mask
defining which pixels are valid and which should be ignored.  The
signal and variance images are represented internally as 32-bit
floating point values.  The variance and mask images may either
be provided by the user, or they may be automatically generated from
the input image, based on configuration-defined values for the image
gain, read-noise, saturation, and so forth.  \textadd{Within the IPP analysis,
we normally use images which are equivalent to the digital numbers
(scaled by the detrend images), but as long as the variance image is
constructed in a consistent fashion, \ippprog{psphot} can use images
in electron, calibrated flux units or other conventions (though this would
require some tuning of configuration parameters).}  For the function-call
form of the program, the flux image is provided in the API, and
references to the mask and variance are provided in the configuration
information.  As in the stand-alone C-program, the variance and mask may
be constructed automatically by \ippprog{psphot}.

The mask is represented as a 16-bit integer image in which a value of
0 represents a valid pixel.  Each of the 16 bits define different
reasons a pixel should be ignored, \textadd{listed in
  Table~\ref{tab:mask_values}}.
This allows us to optionally
respect or ignore the mask depending on the circumstance.  For
example, in some cases, we ignore saturated pixels completely while in
other circumstances, it may be useful to know the flux value of the
saturated pixel.  In addition, the mask pixels are used to define the
pixels available during a model fit; those which should be ignored for
that specific fit are `marked' by setting a special bit (\code{MARK = 0x8000}).
The initial mask, if not supplied by the user or library calls, is
constructed by default from the image by applying three rules: 1)
Pixels which are above a specified saturation level are marked as
saturated.  The level is specified by the camera format keyword
\code{CELL.SATURATION}, which may specify a value or define a header
keyword which in turn specifies the value in the image header.  In the
case of PS1 PV3, the header keyword \code{MAXLIN} specifies the
saturation level for each chip \citep[see][]{waters2017}. 2) Pixels
which are below a user-defined value (\code{CELL.BAD} = 0 for PV3) are considered unresponsive and
masked as dead.  3) Pixels which lie outside of a user-defined coordinate window
are considered non-data pixels (\eg, overscan) and are marked as
invalid.  (\ippprog{psphot} recipe keywords \code{XMIN}, \code{XMAX},
\code{YMIN}, \code{YMAX}, all set to 0 for PS1 PV3 -- invalid pixels
were specified for PS1 PV3 with a supplied mask image
\citep[see][]{waters2017}.

The library functions used by \ippprog{psphot} understand two types of
masked pixels: ``bad'' and ``suspect''.  Bad pixels are those which
should not be used in any operations, while suspect pixels are those
for which the reported signal may be contaminated or biased, but may
be usable in some contexts.  For example, a pixel with poor charge
transfer efficiency is likely to be too untrustworthy to use in any
circumstance, while a pixel in which persistence ghosts have been
subtracted might be useful for detection or even analysis of brighter
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}.  

{\bf An important point to note is that \ippprog{psphot} does not
  attempt to interpolate or replace bad pixel values in the images
  before processing.  The GPC1 images have quite extensive masking due
  to both defects and natural gaps between detectors and amplifier
  regions.  On average, roughly 71\% of the full useable field-of-view
  is covered with valid pixels (See Paper III for more discussion).
  Any attempt to interpolate bad pixels would be quickly overwhelmed
  by these extensive regions.  Rather than attempt to fill in the bad
  pixels, we rely in the PS1 PV3 processing on the fact that regions
  on the sky were observed many times.  Thus, it should be noted that
  model-fitting measurements (which can naturally ignore masked
  pixels) should generally be more reliable than aperture-like
  measurements for single exposures.  Aperture-like measurements from
  the stacks do not suffer from this masking issue. See also the
  discussion of the \ippmisc{PSF_QF} and \ippmisc{PSF_QF_PERFECT}
  parameters for judging the impact of masking on a particular source
  (Section~\ref{sec:psf.model.choice}).}

\begin{table*}

\caption{\label{tab:mask_values} Pixel values for input GPC1 mask
  images used by \nocode{psphot}.  The table gives the bit value used
  to mark the listed effects.  Bits marked as `dynamic' are set for
  each image based on the contents, such as the locations of bright
  stars.  Bits marked as `suspect' represent effects which do not
  definitely affect the photometry, but users should be careful.  The
  mask image headers also list these values.} % \vspace{-0.5cm}
\begin{center}
\footnotesize
\begin{tabular}{lcccl}
\hline
\hline
{\bf Mask Name} & {\bf Mask Value} & {\bf Dynamic?} & {\bf Suspect?} & {\bf Description} \\
\hline
  DETECTOR & 0x0001 & N & N & A detector defect is present. \\
  FLAT     & 0x0002 & N & N & The flat field model does not calibrate the pixel reliably. \\
  DARK     & 0x0004 & N & N & The dark model does not calibrate the pixel reliably. \\
  BLANK    & 0x0008 & N & N & The pixel does not contain valid data. \\
  CTE      & 0x0010 & N & N & The pixel has poor charge transfer efficiency. \\
  SAT      & 0x0020 & Y & N & The pixel is saturated. \\
  LOW      & 0x0040 & Y & N & The pixel has a lower value than expected. \\
  SUSPECT  & 0x0080 & Y & Y & The pixel is suspected of being bad$^1$. \\
  BURNTOOL & 0x0080 & Y & Y & The pixel contain an burntool repaired streak. \\
  CR       & 0x0100 & Y & N & A cosmic ray is present. \\
  SPIKE    & 0x0200 & Y & Y & A diffraction spike is present. \\
  GHOST    & 0x0400 & Y & Y & An optical ghost is present. \\
  STREAK   & 0x0800 & Y & Y & A streak is present. \\
  STARCORE & 0x1000 & Y & Y & A bright star core is present. \\
  CONV.BAD & 0x2000 & Y & N & The pixel is bad after convolution with a bad pixel. \\
  CONV.POOR& 0x4000 & Y & Y & The pixel is poor after convolution with a bad pixel. \\
  MARK     & 0x8000 & X & X & An internal flag for temporarily marking a pixel. \\
\hline
\multicolumn{5}{l}{$^1$ The SUSPECT bit is generic and only
  used if a specific reason cannot be identified.}\\
\multicolumn{5}{l}{~~~It is overloaded on the same bit as BURNTOOL.}\\
\end{tabular}
\end{center}
\end{table*}

The variance image, if not supplied, is constructed by default from
the flux image using the configuration supplied gain and read noise
values to calculate the appropriate Poisson statistics for each pixel.
The parameters are determined based on the camera format keywords
\code{CELL.GAIN} and \code{CELL.READNOISE}, which in the case of PS1
PV3 refer to the header keywords \code{GAIN} and \code{RDNOISE}.  In
this case, the image is assumed to represent the readout from a single
detector, with well-defined gain and read noise characteristics.  This
assumption is not always valid.  For example, if the input flux image
is the result of an image stack with a variable number of input
measurements per pixel (due to masking and dithering), the variance
cannot be calculated from the signal image alone.  It is necessary in
such a case to supply a variance image which accurately represents the
variance as a function of position in the image.

Some image processing steps introduce cross-correlation between pixel
fluxes.  An obvious case is smoothing, but geometric transformations
which redistribute fractional flux between neighboring pixels also
introduces cross-correlations.  In the noise model, it is necessary to
track the impact of the cross correlations on the per-pixel variance.
In the general case, this would require a complete covariance image,
consisting of the set of cross-correlated pixels for each image pixel.
Since a typical smoothing or warping operation may introduce
correlation between 25 - 100 neighboring pixels, the size of such a
covariance image is prohibitive.

%% \note{describe the way we handle covariance}

%% Within the IPP analysis generally, we carry a simplified
%% representation of the impact of covariance on the variance values
%% used in pixel analysis operations.  Whenever image operations
%% introduce covariance by combining information from multiple pixels,
%% we update a matrix tracking the covariance at the image center for
%% a small range of pixels.  

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 defined by the \ippprog{psphot} recipe values
\code{BACKGROUND.XBIN, BACKGROUND.YBIN}, set to 400 pixels
\textadd{($\sim 100$ arcseconds)} for PV3.  Superpixels of size
\code{BACKGROUND.XSAMPLE, BACKGROUND.YSAMPLE} ($2 \times 2$ for PV3)
times larger than this spacing are used to measure the local
background for each background grid point, thus over-sampling the
background spatial variations.  In the interest of speed, a subset of
\code{IMSTATS_NPIX} (10,000 for PV3) randomly selected {\em unmasked}
pixels in these regions are used to determine the background.  The
background value for each superpixel is determined by fitting a
Gaussian distribution to the histogram of pixels values.  

If the image were empty of stars and only contained flux from a
uniform background sky, we would expect the distribution to be Poisson
distributed, and in general in a high-enough signal range to be
essentially Gaussian.  We fit a symmetric Gaussian to all histogram
bins within 15\% of the peak bin value to determine the mean and
standard deviation values for the background.  

If, however, the sky is not empty of stars or other sources, and we
have correctly masked the large majority of non-responsive pixels,
then we expect the flux distribution of the pixels to be asymmetric
with a Gaussian core representing the sky and a tail to the high end
representing the pixels with astronomical source flux contributions.
We would like to determine the mean of the underlying Gaussian without
suffering bias from the stellar flux.  We thus perform a second
Gaussian fit using an asymmetric subset of the histogram pixels,
fitting those histogram bins which are left of the peak but for which
the bin value is greater than 25\% of the peak bin, or right of the
peak but only using those bins for whch the bin value is greater than
50\% of the peak bin value.

\begin{table}
\caption{\label{tab:sky.offset} Comparison of background
  measurement methods.  Backgrounds were measured for simulated images with the given stellar
  density (at the low-density detection threshold) and known
  background level.  The {\tt psphot} technique is less biased at high
stellar densities.} \vspace{-0.5cm}
\begin{center}
% \footnotesize
\begin{tabular}{cccccc}
\hline
\hline
{\bf Density} & {\bf True} & {\bf Image} & {\bf  Image} & {\bf Gauss} & {\bf \tt psphot} \\
{\bf \footnotesize $\log_{10}(\mbox{deg}^{-2}$)} & {\bf Sky} & {\bf Mean} & {\bf  Median} & {\bf Fit} & {\bf \tt Value} \\
\hline
4.2 & 202.8 & 203.3 & 202.8 & 202.8 & 202.9 \\
4.7 & 202.8 & 204.9 & 203.1 & 203.0 & 203.0 \\
5.2 & 202.8 & 210.6 & 204.0 & 203.5 & 203.5 \\
5.7 & 202.8 & 233.9 & 207.4 & 205.4 & 205.3 \\
6.2 & 202.8 & 300.9 & 219.7 & 211.2 & 210.6 \\
6.7 & 202.8 & 534.6 & 286.2 & 242.8 & 233.9 \\
\hline
%\multicolumn{5}{l}{$^1$ a footnote} \\
\end{tabular}
\end{center}
\end{table}

If the fit to the asymmetric lower fraction of the curve is less than
the symmetric fit, but greater than the above lower-bound of the full
symmetric fit, then the lower fraction value is kept as the true mean
sky value for this superpixel.  Table~\ref{tab:sky.offset} shows a
comparison of this technique to several other methods to measure the
sky background using simulated data with a range of stellar
densities. The stellar density listed in the table is the number of stars per
square degree at the $5\sigma$ detection limit {\em in the
  lowest-density image}. In our simulations, we find that as the
stellar density rises to values typical in the Galactic plane regions,
this technique results in a more accurate estimate of the background,
though it still over-estimates the background compared to the truth.

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 source in the
output catalog.  For more details of the background subtraction, see
the discussion in Section~3.11 of \cite{waters2017}.

\textadd{Since the subtraction of the sky model supresses larger-scale
  structures, features such as large galaxies which are comparable to
  the superpixel size are adversely affected by the subtraction.
  Photometry for galaxies larger than $\sim 30$ arcseconds is
  unreliable as a result.  The superpixel size used for the sky model
  in the PV3 analysis was chosen as compromise between the need to
  follow bright features with small spatial scales and the desire to
  measure photometry of galaxies of sizes up to at least 30
  arcseconds.  Features which we wished to suppress include both
  astronomical sources, such as bright nebulosity and the wings of
  bright stars, and non-astronomical sources, such as moonlight and
  other scattered light sources.  In some contexts, we have used a
  finer spacing for the background model, such as in the dedicated
  analysis of the photometry of the Andromeda Galaxy, where we are
  only interested in stellar sources and the analysis is otherwise
  badly affected by the background from this galaxy.}

\subsection{Initial Source Detection}
\label{sec:source.detection}

\subsubsection{Peak Detection}
\label{sec:peaks}

%% is there a ref I can use for the optimal detection? see SDSS docs?

The initial source detection step is focused on finding and
identifying the brighter point sources.  The goals are two-fold: 1) to
select sources which can be used to model the PSF and 2) to subtract
the brighter sources so that fainter sources may be found throughout
the image .

The sources 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.  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
(\code{PEAKS_NSIGMA_LIMIT} = 20.0 for PV3).  A maximum of
\code{PEAKS_NMAX} (5000 for PV3) are found at this stage.

\textadd{For an image with a Gaussian PSF of the same size, this method
  would represent the optimal detection algorithm, equivalent to a
  matched filter.  At this stage, our goal is simply to
  detect the brighter sources, so the exact size and shape of the PSF
  is not critical. }
The
detection efficiency for the brighter sources is not strongly
dependent on the form of this smoothing function.  \textadd{Instead,
  our goal with the smoothing kernel is to reduce our sensitivity to
  pixel-to-pixel fluctuations in the location of the peak of the
  sources in the image.}.  

The local peaks in the smoothed image are found by first detecting
local peaks in each row.  For each peak, the neighboring pixels are
then examined and the peak is accepted or rejected depending on a set
of simple rules.  \textadd{The rules are defined so that we choose a unique set
of peaks which are not immediately adjacent to other peaks.}  First, any peak which is greater than all 8
neighboring pixels is kept.  Any peak which is lower than any of the 8
neighboring pixels is rejected.  Any peak which has the same value as
any of the other 8 pixels is kept {\em if} the pixel $X$ and $Y$ coordinates
are greater than or equal to the other equal-value pixels.  \textmod{This
last rule means that a flat-topped region will result in peaks at
the maximum $X$ and $Y$ corners of the region.}

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
2D quadratic polynomial, and use a very simple bi-quadratic fit to
these pixels.  We use the following function to describe the peak

\[ f(x,y) = C_{00} + C_{10}x + C_{01} y + C_{11} x y + C_{20} x^2 + C_{02} y^2 \]

and write the chi-square equation:

\[ \chi^2 = \sum_{i,j} (F_{i,j} - f(x,y))^2 / \sigma_{i,j}^2 \]

By approximating the error per pixel as the Poisson error on just the peak,
and pulling that term out of the above equation, and recognizing that
the values $X,Y$ in the $3 \times 3$ grid centered on the peak pixel have values
of only 0 or 1, we can greatly simplify the chi-square equation to a
square matrix equation with the following values:
\[
\left( \begin{array}{cccccc}
9 & 0 & 0 & 0 & 6 & 6 \\ 
0 & 6 & 0 & 0 & 0 & 0 \\ 
0 & 0 & 6 & 0 & 0 & 0 \\ 
0 & 0 & 0 & 6 & 0 & 0 \\ 
6 & 0 & 0 & 0 & 6 & 4 \\ 
6 & 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 $3 \times 3$ 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}
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 resulting peak position, ($x_{min}, y_{min}$), is used as the
default starting coordinate for the source.  Later in the
\ippprog{psphot} analysis, improved measurements of the source positions
are calculated as discussed below.

% uses plots.sh in this directory
\begin{figure}[htbp]
  \begin{center}
 \includegraphics[width=\hsize,clip]{\picdir/peaks.pdf}
  \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. Note that this 1D
    illustration is representative of the full 2D path which may be
    followed from one peak to the next.}
  \end{center}
\end{figure}

\subsubsection{Footprints}
\label{sec: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.  \ippprog{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
(\code{PEAKS_NSIGMA_LIMIT}).  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 \code{FOOTPRINT_GROW_RADIUS} (3
pixels for PV3).  Peaks are assigned to the footprints in which
they are contained (note by construction all peaks must be located in
a footprint since the peaks must be above the detection threshold).

For any peak which is not the brightest peak in that footprint it is
possible to reach the brightest peak by following a sequence of the
highest valued pixels between the two peaks.  The lowest pixel along
this \textadd{(potentially meandering)} 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 for PV3) 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}).  \textadd{If more than one such path is
  possible, the path with the highest key col is used for this test.}
In the vicinity of a saturated star, the rule is somewhat more
aggressive 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.
Sources for which the peak is greater than 50\% of the saturation
value require the col to also be a fixed fraction (5\%) of the
saturation below the peak to avoid being marked as locally
insignificant.

Sometimes it is useful to know if a source has a near neighbor which
may be affecting the photometry.  Three flag bits are used to identify
such possible situations.  Peaks which are {\em not} the brightest peak
within a single footprint have the flag bit
\code{PM_SOURCE_MODE2_HAS_BRIGHTER_NEIGHBOR} set.  This is a fairly
common situation.  We also define the following ratio to compare the
flux of the bright source to the flux of a neighbor scaled by
intervening area: $R = \frac{f_{\rm n}}{r^2 f_{\rm p}}$ where $f_{\rm
  n}$ is the flux of the brightest neighbor in the footprint, $f_{\rm
  p}$ is the flux of the source of interest, and $r$ is the separation
between the two sources.  If $R > 1$, the flag bit
\code{PM_SOURCE_MODE2_HAS_BRIGHT_NEIGHBOR_1} is set. If $R > 10$, the flag bit
\code{PM_SOURCE_MODE2_HAS_BRIGHT_NEIGHBOR_10} is set. 

\subsubsection{Centroid and Higher-Order Moments}
\label{sec:moments}

% /data/kukui.3/eugene/psphot.20161214/mana.sh
\begin{figure}[htbp]
  \begin{center}
  \includegraphics[width=0.95\hsize]{{\picdir/FWHM.smooth.trend.v1.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.  Panels (a) - (e)
    corresponds to a different value of $\sigma_w$, corresponding to
    the window FWHM values as marked.  The solid red line is the true
    FWHM of the PSF used to generate the stars (1.4 arcsec in all
    cases).  The blue solid line is the FWHM of the window function.
    The gray dots are the FWHM derived from the measured second
    moments for stars in the image.  The median of this distribution
    (mag $< -10$) is listed as ``obs''.  The ratio of the median FWHM
    to the FWHM of the window function is listed as ``ratio'', while
    the ratio of the median FWHM to the true stellar FWHM is listed as
    ``bias''.  The dotted blue line is the target (65\% of the window
    function).  In this example, we would choose $\sigma_w$ between
    0.5 and 0.8 arcseconds (FWHM between 2.64 and 3.52 arcseconds), so
    the dotted blue line would match the bright end of the gray dots.
    See discussion in the text for the choice of target window.  }
  \end{center}
\end{figure}

Once a collection of peaks has been identified, a number of basic
properties of the sources related to the first, second, and higher
moments are measured.  \textmod{These moments can be used for a crude
  classification of the sources.  As discussed below, the second
  moments are used to select candidate stellar sources to be used in
  modeling the PSF and the exclude `cosmic rays' and extended sources.
  The radial moment is used in the measurement of the Kron magnitudes \citep{1980ApJS...43..305K}.
  The higher-order moments are provided primarily for image quality
  diagnostics.}

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,
centered on the object, with size $\sigma_w$ chosen to be large
compared to the PSF size, $\sigma_{\rm PSF}$.  This window function
reduces the noise of the measurement of the moments by suppressing the
noisy pixels at high radial distance as well as by reducing the
contaminating effects of neighboring stars.  The choice of $\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 the PSF size,
$\sigma^{\prime}_{\rm PSF}$ \textmod{which is different} from the true value due to
the effect of the window function.  The measured value of the PSF size
will be biased high or low depending on both the signal-to-noise of
the source and the size of the window function compared to the true
PSF size.

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 $\sigma_{\rm PSF}$
corresponding to a FWHM of 1.4 arcseconds.  For bright stars, as the
window function $\sigma_w$ is increased, the measured FWHM rises from
an initially under-estimated value to meet the truth value.  For faint
stars, the measured value of the FWHM is initially under-estimated as
well.  However, as the value of $\sigma_w$ increases, the measured
FWHM for faint stars rises, and then over-shoots the truth value,
while the scatter increases.  Thus, for large values of $\sigma_w$,
the result is both a poorly estimated FWHM for the image and a \textmod{trend
with the} signal-to-noise of the star.  We attempt to minimize the
scatter and trends with instrumental magnitude at the cost of overall
bias.

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
\citep[c.f., ``the Brighter-Fatter''
  effect,][]{2014JInst...9C3048A,2015JInst..10C5032G}.  To make a
robust choice for $\sigma_w$, we choose a value such that the measured
value of $\sigma^{\prime}_{\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 a sequence of values
spanning a range guaranteed to contain any reasonable seeing values.
The values are specified in the \ippprog{psphot} recipe as
\code{PSF.SIGMA.VALUES} and have the following values for PS1 PV3: (1,
2, 3, 4.5, 6, 9, 12, 18) pixels $\approx$ (0.26, 0.51, 0.77, 1.15,
1.54, 2.3, 3.1, 4.6) arcseconds.  For each of these $\sigma_w$ values,
we then select candidate PSF stars based on the distribution of the
measured $\sigma^{\prime}_{\rm PSF}$ in the two principal directions:
$\sigma_{x,x}$ and $\sigma_{y,y}$ (see
Section~\ref{sec:psf.source.selection}, below).  For each test value
of $\sigma_w$, we determine the ratio $\rho_\sigma =
\frac{\sigma_{x} + \sigma{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 $\rho_\sigma$ is expected to be 0.65.  We use
an aperture with a radius of 4$\sigma_w$ to select the pixels for the
measurement of the moments.

%% confirmed: PSF_MOMENTS_RADIUS = 4 * MOMENTS_GAUSS_SIGMA (\sigma_w)
%% factor of 4 is hard-wired in psphotSourceStats.c where MOMENTS_GAUSS_SIGMA is set.
%% PSF_MOMENTS_RADIUS used for: moments analysis, Kron analysis (starting radius),
%% radial profile wings (starting radius)

Once $\sigma_w$ has been determined, moments are measured as defined
below.

\begin{eqnarray}
x_0      & = & \frac{1}{S} \sum_i w_i (f_i - s_i)x_i \\
y_0      & = & \frac{1}{S} \sum_i w_i (f_i - s_i)y_i \\
M_{xx}   & = & \frac{1}{S} \sum_i w_i (f_i - s_i)(x_i - x_0)^2 \\
M_{xy}   & = & \frac{1}{S} \sum_i w_i (f_i - s_i)(x_i - x_0)(y_i - y_0) \\
M_{yy}   & = & \frac{1}{S} \sum_i w_i (f_i - s_i)(y_i - y_0)^2 \\
M_{xxx}  & = & \frac{1}{S} \sum_i \frac{w_i}{r_i} (f_i - s_i)(x_i - x_0)^3 \\
M_{xxy}  & = & \frac{1}{S} \sum_i \frac{w_i}{r_i} (f_i - s_i)(x_i - x_0)^2(y_i - y_0) \\
M_{xyy}  & = & \frac{1}{S} \sum_i \frac{w_i}{r_i} (f_i - s_i)(x_i - x_0)(y_i - y_0)^2 \\
M_{yyy}  & = & \frac{1}{S} \sum_i \frac{w_i}{r_i} (f_i - s_i)(y_i - y_0)^3 \\
M_{xxxx} & = & \frac{1}{S} \sum_i \frac{w_i}{r^2_i} (f_i - s_i)(x_i - x_0)^4 \\
M_{xxxy} & = & \frac{1}{S} \sum_i \frac{w_i}{r^2_i} (f_i - s_i)(x_i - x_0)^3(y_i - y_0) \\
M_{xxyy} & = & \frac{1}{S} \sum_i \frac{w_i}{r^2_i} (f_i - s_i)(x_i - x_0)^2(y_i - y_0)^2 \\
M_{xyyy} & = & \frac{1}{S} \sum_i \frac{w_i}{r^2_i} (f_i - s_i)(y_i - y_0)(y_i - y_0)^3 \\
M_{yyyy} & = & \frac{1}{S} \sum_i \frac{w_i}{r^2_i} (f_i - s_i)(y_i - y_0)^4
\end{eqnarray}
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 sums are performed
over all (unmasked) 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.

\textadd{The motivation of measuring these higher order moments was to
  select exposures with image quality problems.  For example, trefoil
  caused by errors in the collimation and alignment can in principle
  be detected with the third-order moments.  In our experience, these
  statistics can be used to select some images with such problems, but
  we have not been able to use these values to exclude poor images
  from the data processing.  If we were to reject images based on
  these moments, we would reject too many images with image quality
  issues that are not so poor as to preclude a useful analysis.  A
  future machine-learning based analysis starting with these moments
  might potentially provide a better rejection statistic, but such
  work is beyond the scope of this article.}

For sources with peak flux above the saturation limit, the moments are
generally poorly measured if the aperture defined by $\sigma_w$ is
used.  For these sources, the quality of the measurment is compromised
by the saturation.  However, it is still useful to estimate the first
and second moments of the source in order to allow a crude measurement
of the brightness from the wings of the source.  In this case, a
larger aperture, 3 times the standard aperture, is used to make a
crude estimate.  For such sources, the flag bit
\code{PM_SOURCE_MODE_BIG_RADIUS} is set and the source is ignored in
all analyses below except for the analysis applied to very bright
stars (Section~\ref{sec:very.bright.star}).

If the measured centroid coordinates ($x_0, y_0$) differ from the peak
coordinates \textmod{by} a large amount (1.5$\sigma_w$), then the peak is
identified as being of poor quality and is skipped in further
analyses; the flag bit
\code{PM_SOURCE_MOMENTS_FAILURE} is set for such sources.  In such
a case, it is likely that the `peak' was identified in a region of
flat flux distribution or many saturated or edge pixels.  During the
analysis of the moments, the background (``sky'') model is also examined for the
location of each source.  The value of the background and the variance
of the background are recorded for each source.  In some cases, the
sky model or the variance is not well defined at the location of a
specific sources (e.g., due to an extrapolation failure).  In these
cases, the flag bits \code{PM_SOURCE_SKY_FAILURE} or
\code{PM_SOURCE_SKYVAR_FAILURE} are set as appropriate and the
measurement of the moments is skipped.  

In addition to the moments above, the 1st and half-radial moments,
$M_r$ and $M_h$ as defined below, are calculated:
\begin{eqnarray}
\label{eqn:first.radial.moment}
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. 

With the first radial moment, we can calculate a preliminary Kron
radius and magnitude.  \textadd{Kron magnitudes are provided as an option for
galaxy photometry.  In addition, the comparison of Kron and PSF
magnitudes is useful as a star-galaxy separator.} The Kron radius
\citep{1980ApJS...43..305K} is defined the be 2.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\sigma_w$ if that cannot be
  determined.  $R_{\rm max}$ is set to the size of the moments
  aperture, $4\sigma_w$.  These Kron measurements are performed for
  all sources with a valid set of moments.  At this stage, the
  measurement of the Kron parameters are preliminary since the
  aperture has been chosen as a fixed size relative to the size of the
  PSF.  At a later stage, higher-quality Kron parameters appropriate
  to galaxies are measured with more care paid to the exact aperture
  used (Section~\ref{sec:kron.mags}).

% $\sigma_w$ is saved as MOMENTS_GAUSS_SIGMA
% the aperture radius is saved as PSF_MOMENTS_RADIUS

\subsection{PSF Determination}
\label{sec:PSF.Model}

\subsubsection{PSF Model vs Source Model}
\label{sec:Source.Model}

The point spread function (PSF) of an image describes the shape of all
unresolved sources in the image.  In a typical wide-field image, the
shape of unresolved sources varies as a function of position in the
image.  The full PSF thus needs to include a model with parameters
which vary across the image.

The PSF used by \ippprog{psphot} consists of an analytical function
combined with a pixelized representation of the residual differences
between the analytical model and the true PSF.  Both the shape
parameters of the analytical model and the pixelized residual
differences are allowed to vary in two dimensions across the images.

Within \ippprog{psphot}, several analytical models may be used to
describe the smooth portion of the PSF, but all share a few common
characteristics.  As an example, a simple model consists of a 2-D
elliptical Gaussian:
\begin{eqnarray}
\label{eqn:2d.gaussian}
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}
Here the model parameters consist of the centroid coordinates ($x_o,
y_o$), the elliptical shape parameters ($\sigma_x, \sigma_y,
\sigma_{\rm xy}$), the model normalization ($I_o$) and the local value
of the background ($S$).  

A specific source will have a particular set of values for the model
parameters, some of which depend on the PSF model and the position of
the source in the image, while the rest are unique to the individual
source.  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$), though as
noted below (Section~\ref{sec:nonlinear.psf.model}), in practice we do
not allow the sky to be fitted independently since we subtract the
background model.  Thus the shape parameters are each a function of
the source centroid coordinates:
\begin{eqnarray}
\sigma_x    & = & f_1(x_{\rm ccd},y_{\rm ccd}) \\
\sigma_y    & = & f_2(x_{\rm ccd},y_{\rm ccd}) \\
\sigma_{xy} & = & f_3(x_{\rm ccd},y_{\rm ccd}).
\end{eqnarray}
\ippprog{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
sources within a subregion of the image.  In the latter case, the
value at a specific coordinate in the image is determined \textmod{via
  bi-linear} 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 \textadd{(grid cells of
  size $\sim 3.4$ arcminutes)}.  For the earlier PV2 analysis, the
maximum grid sampling was $3\times 3$ per GPC1 chip image
\textadd{(grid cells of size $\sim 6.9$ arcminutes)}.  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.

Figure~\ref{fig:iq.exposure} illustrates the 2D variations in the PSF
shapes seen in PS1 data.  This figure shows the FWHM, $e_1$, and $e_2$
polarizations of the stars as a function of position in 4 exposures.
For images with good image quality, variations of the PSF shape due to
the optical aberrations can be see.  The optical aberrations vary as
the active collimation and alignment are adjusted and as the focus
changes.  These aberrations are coupled to the piston of the chips,
which have been adjusted to crudely follow the focal surface
\citep{chambers2017}.  During regular operations, image with large
PSFs are usually caused by the atmosphere (seeing) or by telescope
tracking errors, both of which result in common shapes across the
field of the camera.  In the figure, the top panel shows the
circularization of the PSF due to the atmosphere washes out the
lower-level variations caused by the optics.

Examples of 2D PSF variations.
    Each row represents an exposure.  The left-most column shows the
    distribution of FWHM across the camera; the median value in
    arcseconds is given in the inset.  The middle column gives the
    $e_1$ polarization measured from second moments (see
    Section~\ref{sec:lensing.params} while the right column gives the
    $e_2$ polarization.

% \note{write up the fitting process to define the grid?}

Several analytical functions which are likely candidates to describe
the smooth portion of the PSF are available in \ippprog{psphot}:
\begin{itemize}
\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]}, $\alpha > 1.25$
\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}
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.
\textadd{For these PSF models, the functions are evaluated at the pixel center.
Unlike some galaxy model representations (see
Section~\ref{sec:galaxy.conv.fit} ), the first derivatives of these
functions approach zero as the radius approaches zero, so sub-pixel
integration is not necessary.}
A user may choose to try more than one analytical function for a given
image.  As discussed below (Section~\ref{sec:psf.model.choice}),
\ippprog{psphot} can automatically choose the best model based on the
quality of the PSF fits.

% Figure 3: ** repaired PDF text **
% pueo:psphot.iq.20200413/mana.sh : show.e12 (iq.exposures.pdf)
\begin{figure}[htbp]
  \begin{center}
  \includegraphics[width=\hsize]{{\picdir/iq.exposures}.\plotext}
  \caption{\label{fig:iq.exposure} Examples of 2D PSF variations.
    Each row represents an exposure.  The left-most column shows the
    distribution of FWHM across the camera; the median value in
    arcseconds is given in the inset.  The middle column gives the
    $e_1$ polarization measured from second moments (see
    Section~\ref{sec:lensing.params}) while the right column gives the
    $e_2$ polarization. }
  \end{center}
\end{figure}

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.

The analytical models in \ippprog{psphot} are written with a high degree
of code abstraction making it relatively easy to add different
analytical models to the software.  The same portion of code used to
describe the analytical portion of the PSF sources is also used to for
galaxy models. 

% moffat : 1969A&A.....3..455M
% buonanno : 1983A&AS...51...83B

% Figure 4:
% /data/kukui.3/eugene/psphot.20161214/mana.sh
\begin{figure}[htbp]
  \begin{center}
  \includegraphics[width=\hsize]{{\picdir/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 red and blue points are individual pixel
    values.  The black line shows the PSF model with radial trend of
    the form $(1 + \kappa r^2 + r^{3.33})^{-1}$.  The models use a 1D
    average of the 2D analytical portion of the PSF models fitted
    to these specific stars in their standard analysis.  }
  \end{center}
\end{figure}

Once the smooth component of the PSF has been fitted with an
analytical model, a pixel representation of the residuals is
generated.  This representation is constructed as an image of the
expected residuals for any position in the image.  The value of each
pixel in the image model is determined from 2D fits to the measured
residuals of the PSF stars.  

The residual model is calculated using the residuals for all PSF
stars.  The residuals (and their errors) for each star are
renormalized by the flux of the star to put them on a consistent flux
scale.  For each PSF star, all pixels within a user-specified radius
(\code{PSF.RESIDUALS.RADIUS = 9}) are selected for the measurement.  \textmod{For a
given pixel in the model, the value is calculated from the 4 closest
pixels in the PSF stars via bi-linear interpolation.} Pixels may be used in
this analysis if their signal-to-noise exceeds a user-defined limit.
For the PV3 $3\pi$ analysis, we allowed all pixels within the
user-specified radius, not limiting on the basis of the
signal-to-noise.

Pixels for a given star which are more than a number of sigmas
(\code{PSF.RESIDUALS.NSIGMA = 3.0}) deviant from the median value of
the pixels from all stars are rejected.

If no spatial variation is allowed, the mean or median value is
calculated for the model pixel based on the user-specified mean
statistic (\code{PSF.RESIDUALS.STATISTIC = ROBUST_MEDIAN}).

If spatial variation is requested, then the pixel values are fitted to
a linear model:
\[
\begin{array}{lll}
R[(x_{\rm mod},y_{\rm mod})][(x_{\rm ccd},y_{\rm ccd})] & = & R_o[(x_{\rm mod},y_{\rm mod})] \\
& + & R_x[(x_{\rm mod},y_{\rm mod})] x_{\rm ccd} \\
& + & R_y[(x_{\rm mod},y_{\rm mod})] y_{\rm ccd} \\
\end{array}
\]
where $R[(x_{\rm mod},y_{\rm mod})][(x_{\rm ccd},y_{\rm ccd})]$ is the
\textmod{value of the residual for model} pixel $(x_{\rm mod},y_{\rm mod})$ for a star with
centroid at image pixel $(x_{\rm ccd},y_{\rm ccd})$.  \textmod{The parameters
$R_o, R_x, R_y$ are the elements of the 2-D linear fit for each pixel $(x_{\rm mod},y_{\rm mod})$
in the model. }

\subsubsection{Candidate PSF Source Selection}
\label{sec:psf.source.selection}

% Figure 5:
% /data/kukui.3/eugene/psphot.20161214/mana.sh
\begin{figure}[htbp]
  \begin{center}
  \includegraphics[width=\hsize]{{\picdir/moment.class}.\plotext}
  \caption{\label{fig:moment.class} Illustration of PSF star selection
    using the second moments. \textadd{Each point represents the
      second moments in the $X_{\rm ccd}$ and $Y_{\rm ccd}$ directions
      for sources measured in one chip (XY32) from a particular PS\,1
      exposure (o6065g0428o)}.  The dominant clump is located in this
    diagram \textadd{to identify the stars.}  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}

The first stage of determining the PSF model for an image is to
identify a collection of sources in the image which are {\em likely}
to be unresolved (i.e., stars).  \ippprog{psphot} uses the source sizes as
estimated from the second moments to make the initial guess at a
collection of unresolved sources.  At this point, the program has
measured the second order moments for all sources identified by their
peaks, as well as an approximate signal-to-noise ratio, above the
bright threshold.  All sources with a S/N ratio greater than a
user-defined parameter (\code{PSF_SN_LIM} = 20.0 for PS1 PV3) are
selected by \ippprog{psphot}, though sources which have more than a
certain number of saturated pixels are excluded at this stage.  The
program then examines the 2-D plane of $M_{x,x}, M_{y,y}$ in search
of a concentrated clump of sources (see
Figure~\ref{fig:moment.class}).  To do this, it constructs an
artificial image with pixels representing the value of $M_{x,x},
M_{y,y}$, using $0.1 \sigma^2_w$ as the size of a pixel in this
artificial image.  The binned $M_{x,x}, M_{y,y}$ plane is then
examined to find a significant peak.  Unless the image is extremely
sparse, such a peak will be well-defined and should represent the
sources which are all very similar in shape.  Other sources 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, sources which have second-moments very close to 0 are
ignored.  For these sources, the flag bit \code{PM_SOURCE_MODE_DEFECT}
is set.

Once a peak has been detected in this plane, the centroid and second
moments of this peak are measured.  All sources which land within 2
pixels of this centroid are selected as candidate PSF sources in the
image.

When the second moments are measured, \ippprog{psphot} also counts the
number of saturated pixels within the analysis aperture.  If more than
a single saturated pixel is found, and if the second moments of that
object are more than one standard deviation larger than the clump
identified above, this source is identified as a highly saturated star
and marked with the flag bit \code{PM_SOURCE_MODE_SATSTAR}.  Sources
which have more than a single saturated pixel, but for which the
second moments do not exceed the above limits are marked as likely
saturated regions (e.g., bleed trails).  These sources are skipping in
most additional analyses and are marked with the flag bit
\code{PM_SOURCE_MODE_SATURATED}.

\subsubsection{Candidate PSF Source Model Fits}
\label{sec:psf.model.choice}
% \note{link to psLibADD}

% Madsen:
%% http://www2.imm.dtu.dk/pubdb/views/edoc_download.php/3215/pdf/imm3215.pdf

All candidate PSF sources are then fitted with the selected source
model, allowing all of the parameters (PSF and independent) to vary in
the fit.  The software uses the Levenberg-Marquardt minimization
technique \citep[e.g.,][]{1992nrca.book.....P,Madsen} for the non-linear fitting.  Non-linear
fitting can be very computationally intensive, particularly if the
starting parameters are far from the minimization values.  The first
and second moments are used to make a good guess for the centroid and
shape parameters for the PSF models.  Any sources which fail to
converge in the fit are flagged as invalid.

{\bf
To generate the initial guess, the second moments are converted to the equivalent sigma values for a
2D elliptical Gaussian contour using the following transformations
inspired by \cite{sextractor,1980JBIS...33..323S}.  
First, we calculate the sigma values in the major ($\sigma_a$) and
minor ($\sigma_b$) axis directions, along with the position angle
$\theta$ from the moments using:
\begin{eqnarray}
\theta   & = & \frac{1}{2} \arctantwo (2 M_{xy}, g_2) \\
\sigma_a & = & \sqrt{\frac{g_1 + g_3}{2}} \\
\sigma_b & = & \sqrt{\frac{g_1 - g_3}{2}}
\end{eqnarray}
where the function $\arctantwo (y,x)$ returns the arctangent in the
proper quadrant (e.g,. as implemented by the \code{atan2(y,x)}
function in C) and the intermediate values $g_1$, $g_2$, $g_3$ are given by:
\begin{eqnarray}
g_1 & = & M_{xx} + M_{yy} \\
g_2 & = & M_{xx} - M_{yy} \\
g_3 & = & \sqrt{g_2^2 + 4 M_{xy}^2}
\end{eqnarray}
Since the moments may be noisy, the calculated value of $\sigma_b$ can
be numerically invalid if $g_3 > g_1$, a situation which is especially
likely for highly ellongated sources.  We avoid this situation by
limiting the axial ratio to a maximum of 20 (setting $\sigma_b$ to
$\sigma_a / 20$ if the expected axial ratio would be greater than this
limit).  The selected value of 20 is somewhat ad-hoc, chosen based on
failures in real images.  A more careful examination of the trade-off
space would be worthwhile in the future.

With $\sigma_a$, $\sigma_b$, $\theta$ in hand, we can now transform
these values to the parameters of our fits, $\sigma_x$, $\sigma_y$,
$\sigma_{\rm xy}$ (Eqn~\ref{eqn:2d.gaussian} above).  This transformation
can be determined by rotating the 2D Gaussian equation, yielding:
\begin{eqnarray}
\sigma_x^{-2}  & = & \sigma_a^{-2} \cos^2 \theta + \sigma_b^{-2}\sin^2 \theta \\
\sigma_y^{-2}  & = & \sigma_b^{-2} \cos^2 \theta + \sigma_a^{-2}\sin^2 \theta \\
\sigma_{\rm xy} & = & \frac{1}{2} \sin (2 \theta) (\sigma_b^{-2} - \sigma_a^{-2})
\end{eqnarray}
In fact, since the calculated second moments have been measured with a
window function applied (see discussion in Section~\ref{sec:moments}), we instead
use the measured value of $M_r$ (Eqn~\ref{eqn:first.radial.moment}), the first
radial moment as the major axis size for the Gaussian ($\sigma_a$), retaining
the position angle and axial ratio from the calculation above.  We use
these guess parameters for all version of the PSF analytical models,
despite the fact that for the versions which are not approximations of
Gaussians these guess values will be systematically incorrect.  
It would be worthwhile in the future to tweak the guesses for
the different model version to speed up the convergence.
}

% https://www.astromatic.net/pubsvn/software/sextractor/trunk/doc/sextractor.pdf

For the resulting collection of source model parameters, the
PSF-dependent parameters of the models are all fitted as a function of
position using either the 2-D polynomial or the gridded 
representation described above.  The maximum order of these fits depends on the number
of PSF sources (see Table~\ref{tab:psf.order.nstars}).  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 sources
which are not well described by the PSF model (e.g., galaxies which
snuck into the sample).  Sources whose model parameters are rejected
by this iterative fitting technique are also marked as invalid PSF
sources and ignored in the later PSF model fitting stages.  Sources
which are actually used to define the PSF model for a given image have
the flag bit \code{PM_SOURCE_MODE_PSFSTAR} set.

The order of the fit or number of grid samples is modified if the
number of stars available for the fit is insufficient to justify the
highest value.  Regardless of the requested order, if the number of
stars is below the following limits, the order is limited as shown in
Table~\ref{tab:psf.order.nstars}.  Note that the number of grid cells
in one dimension is one greater than the equivalent polynomial order.

\begin{table}
\caption{\label{tab:psf.order.nstars} Minimum number of stars required
  for a given order of the PSF 2D variations, or for the given number of grid cells.} % \vspace{-0.5cm}
\begin{center}
\begin{tabular}{llll}
\hline
\hline
{\bf Minimum }    & {\bf Order} & {\bf Number of}  & {\bf Cell Size} \\
{\bf \# of Stars} &             & {\bf Grid Cells} & {\bf (arcmin) } \\
\hline
 16 &  1 &  4 & 10.3 \\
 54 &  2 &  9 &  6.9 \\
128 &  3 & 16 &  5.1 \\
300 &  4 & 25 &  4.1 \\
576 &  5 & 36 &  3.4 \\
\hline
\end{tabular}
\end{center}
\end{table}

All of the PSF-candidate sources are then re-fitted using the PSF
model to specify the PSF-dependent model parameter values for each
source.  For example, in the case of the elliptical Gaussian model,
the shape parameters ($\sigma_x, \sigma_y, \sigma_{xy}$) for each
source are set by the coordinates of the source centroid and fixed
(not allowed to vary) in the fitting procedure.  The resulting fitted
models are then used to determine a metric which tests the quality of
the PSF model for this particular image.

% For each model test, the above
% corrected ApResid scatter is measured.  The PSF model function with
% the smallest value for the ApResid scatter is then used by
% \ippprog{psphot} as the best PSF model for this image.  

{\bf \ippprog{psphot} allows a collection of PSF model functions to be
tried on all PSF candidate sources.  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 tested.  The metric used by \ippprog{psphot} to assess the
PSF model is the scatter in the differences between the aperture and
fit magnitudes for the PSF sources.  This difference is a critical
parameter for any PSF modeling software as it is a measurement of how
well the PSF model captures the flux of the star.  Aperture photometry
is measured for a circular aperture with a radius of
\code{PSF_APERTURE_SCALE} (4.5 for PV3) times $\sigma_w$
(Section~\ref{sec:moments}).  The average aperture correction ($m_{\rm
  AP} - m_{\rm PSF}$) is measured and, if multiple PSF model types are
selected, the PSF model with the minimum clipped scatter in this
statistic is chosen for the image.  For the PV3 analysis, however, only the
\code{PS1_V1} model function was used.}

An approximate aperture correction is measured at this stage, with a
more detailed correction measured after all source analysis is
performed (see Section~\ref{sec:aperture.correction}).  Sources for
which the aperture magnitude is measured have the flag bit
\code{PM_SOURCE_MODE_AP_MAGS} set.  These aperture magnitudes are
stored in the DVO field \ippdbtable{Measure.Map} and exported to the
PSPS as a flux in Janskies in the field \ippdbtable{Detection.apFlux}.
The radius (in arcseconds) of the aperture used for each exposure is
reported in PSPS as \ippdbtable{Detection.apRadius}, while the
unmasked fraction of the aperture is reported in PSPS as
\ippdbtable{Detection.apFillF}.

\textmod{As noted above (Section~\ref{sec:image.preparation}), we do not
attempt to replace or interpolate masked pixel values.  Aperture
photometry measurements of objects which include masked pixels are
thus inaccurate.  For a stellar object, the amount of error is a
function of how close the masked pixels are to the core of the PSF.
To provide guidance, when the PSF and aperture photometry for a source
is measured, two additional quantities are measured which are useful
to assess the impact of masking.}  First, the mask image is examined
and the number of unmasked pixels is summed, weighted by the
normalized PSF model.  The resulting quantity, \code{PSF_QF} has a
value between 0.0 (totally masked) and 1.0 (totally unmasked).
Elsewhere in the IPP system, we use this value to filter out
detections which are unreliable due to the masking.  For a generous
cut, leaning toward completeness at the cost of some lower quality
measurements, \code{PSF_QF} $> 0.85$ is used in some contexts; in
other cases, we require \code{PSF_QF} $> 0.95$ to ensure a
high-quality measurement \citep[see for example the calculation of
  average photometry in][]{magnier2017.calibration}.  The second
quantity is related to the first: \code{PSF_QF_PERFECT} uses all mask
values to assess the quality factor, while \code{PSF_QF} uses only the
``bad'' mask bit values (see Section~\ref{sec:image.preparation}).

Several flag bits are raised based on statistics which are similar to
the \code{PSF_QF} measurement.  First, \ippprog{psphot} calculates the
normalized, PSF-weighted fraction of pixels which are masked due to
one of the following four causes: a diffraction spike (\code{SPIKE}),
the core of a saturated star (\code{CORE}), burntool-subtracted region
(\code{BURNTOOL}), or a pixel for which, due to interpolation or
convolution, a significant fraction of the pixel flux comes from a
masked pixel.  These masking conditions are all treated as ``suspect''
by \ippprog{psphot}, which means they are {\em included} in the analysis of
the source pixels.  However, since they may potentially affect the
photometry (or astrometry), it is useful to note of a source has a
non-trivial fraction of these poor mask pixels.  If the normalized
PSF-weighted fraction of pixels masked due to any of these four
conditions exceeds 25\%, then one of the following bits is raised for
the corresponding condition:  \code{PM_SOURCE_MODE2_ON_SPIKE},
\code{PM_SOURCE_MODE2_ON_STARCORE},
\code{PM_SOURCE_MODE2_ON_BURNTOOL},
\code{PM_SOURCE_MODE2_ON_CONVPOOR}.  In addition, the following flag
bits may also be raised if the central pixel
of a source lands on a pixel masked for a diffraction spike
(\code{PM_SOURCE_MODE_ON_SPIKE}), an optical ghost
(\code{PM_SOURCE_MODE_ON_GHOST}), or off the active pixels of the CCD (\code{PM_SOURCE_MODE_OFF_CHIP}).

\subsection{Bright Source Analysis}

Once a PSF model has been determined, the brighter sources in the
image may be analyzed in detail.  The goals in this stage are (1) to
determine the fluxes and positions of the bright stellar sources with
high precision appropriate to their high signal-to-noise and (2) to
characterize the bright source flux profiles sufficiently well that
they may be subtracted from the image to allow for the clean detection
of the fainter sources.  Note that as the analysis proceeds, there are
several stages in which the 2D flux models for all sources are
subtracted from the image, and individual sources are replaced in the
image for a particular analysis step and then removed again.  The flux
limit for this analysis stage is user-defined as a signal-to-noise
value.  In the PV3 analysis of the $3\pi$ survey data, this limit was
set to a signal-to-noise ratio of 20.0.

% maybe drop this discussion? too much detail?
In order to allow for multiple threads to process a single image, the
pixels in an image are divided into a grid of superpixels \textadd{(note that
these superpixels are not the same as those used for either the
background model or the PSF parameter variations)}.  The
superpixels are assigned to one of four groups so that each superpixel
in a group is well separated from the other superpixels of that group.
The analysis of the image proceeds in 4 steps, one for each of these
groups.  Each of the superpixels in the first group is assigned to a
single thread until all threads are assigned.  A single thread is
responsible for the analysis of sources which land within their
current superpixel, as determined by the centroid coordinates.  Since
the superpixels in a given thread group are not contiguous by
construction, sources near the edge of a superpixel can be analysed by
considering the nearby pixels from neighboring superpixel (guaranteed
not to be in the current thread group).

% \note{explain number of superpixels (psphotThreadTools.c)}

As the threads complete their analysis, they are assigned the next
unfinished superpixel in the active group.  When all superpixels in
one group have been processed, then the superpixels in the next group
can start.  This strategy allows the threading to process sources
which may be extended without the danger that two threads are actively
touching the same pixels.  For the PV3 analysis, 4 threads were used
for most processing tasks.

\subsubsection{Very Bright Stars}
\label{sec:very.bright.star}

The standard \ippprog{psphot} PSF modeling code fails to fit the wings of
highly saturated stars, especially if the core of the star is too
contaminated by saturated pixels.  For stars with more than a single
saturated pixel, we model the radial profile of the logarithmic
instrumental flux in logarithmically spaced radial bins.  For each
radial bin, we determine the median of the log-flux.  This median
profile is then interpolated to generate the full radial flux
distribution.  Note that in the case of very saturated stars, pixels
in the central regions are largely masked, because they are
saturated.  Thus in these cases, the PSF-weighted masked fraction (see
Section~\ref{sec:psf.model.choice}) is generally quite low or 0.0.
Sources for which this radial profile is subtracted have the flag bit
\code{PM_SOURCE_MODE2_SATSTAR_PROFILE} set.

% logRdel = 0.1
% logRmax = log(320)

\subsubsection{Fast Ensemble PSF Fitting}
\label{sec:ensemble.fitting}

Before the detailed analysis of the sources is performed, it is
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.  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 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
results in a matrix equation:
\[ 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 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.  Direct inversion of the matrix would
be computationally very slow.  However, an iterative 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 process is repeated several times
to achieve convergence.  Convergence is quick (a few iterations)
because of the highly diagonal matrix with small off-diagonal terms:
the dot product of source $i$ and source $j$ is 1 where $i = j$ and
much less than 1 where $i \neq j$.

Once a solution set for $A_i$ is found, all of the sources are
subtracted from the image by applying these values to the unit-flux
PSF.  Sources for which a PSF model has been fitted (whether or not
this is retained as the best model in the end) has the flag field
\code{PM_SOURCE_MODE_PSFMODEL} set.  All sources which are included in
the ensemble linear fit have the flag bit
\code{PM_SOURCE_MODE_LINEAR_FIT} set, including those for which the
model is not the PSF.

\subsubsection{Radial Profile Wings}
\label{sec:radial.profile}

We attempt to measure the radial profile of sources in order to find
the radius at which the flux of the source 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 source 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 \textmod{output FITS catalog files}, but
not sent to the PSPS.  The sky radius value is used below in the
calculation of the Kron magnitude.

\subsubsection{Kron Magnitudes}
\label{sec:kron.mags}

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, following the algorithm described by \cite{sextractor},
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 to 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
1/2$\times$ the effective surface brightness of a point 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 source has 180\degree\ symmetry, this operation has no impact.
However, if one of the two pixels is unusually high, the value will be
suppressed 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.  \textadd{We found it necessary to apply this filter because,
although the source models have been subtracted, at this point in the
analysis, only PSF models have been used.  Thus extended objects
(galaxies) can leave behind significant amounts of flux to contaminate
the neighbors.}

\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, \ippprog{psphot} uses these two
measurements, along with some additional pixel-level analysis, \textmod{for
classification based on source size.}  Sources identified as
extended will be fitted with a galaxy model (or possibly another type
of extended source model in special cases).  \textadd{If the source 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 $\delta M_{\rm KP} = m_{\rm Kron} - m_{\rm PSF}$,
the difference between the PSF and Kron magnitudes, is calculated for
each source.  The median of $\delta M_{\rm KP}$ is calculated for the
PSF stars.  This median is subtracted from $\delta M_{\rm KP}$ 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 the configuration value \code{PSPHOT.EXT.NSIGMA.LIMIT} (3.0 for PV3), the source is
considered to be extended and the flag bit
\code{PM_SOURCE_MODE_EXT_LIMIT} is set for the source.

\textmod{We decided to use the $\delta M_{\rm KP}$ metric for this
  assessment after we tested several possible star-galaxy separation
  statistics.  We found that the Kron-PSF comparison was more reliable
  than second-moment and first-radial-moment based measurements.  In
  addition, since we needed a statistic which could be calculated
  relatively quickly on every detected source, we rejected using a
  galaxy model fit for the star-galaxy separator.}

Cosmic rays are identified by a combination of the Kron magnitude and
the second-moment width of the source 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} < 0.8$ pixels$^2$ and the signal-to-noise of the
flux measured in the moments analysis $> 7$, then the source is
identified as a cosmic ray and the associated pixels are masked.
These values are tuned empirically for the PV3 analysis based on
cosmic rays identified in the GPC1 images.  Sources which are
determined to the a cosmic ray in this manner have the flag bit
\code{PM_SOURCE_MODE_DEFECT} set.

The pixels of any suspected cosmic ray identified above are examined
in additional detail to make a final judgement.  The Laplacian edge
detection algorithm based on \cite{2001PASP..113.1420V} is used to
check for sharp edges in the flux distribution.  If the sharpness
exceeds a defined limit, then the pixels are masked and the flag bit
\code{PM_SOURCE_MODE_CR_LIMIT} is set for the source.

% Mminor < 0.8 && SN > 7

% for dynamic CR parameters, use object with Mminor < 1.2 and Mkron <
% -5.5 to assess the distribution

\subsubsection{Full PSF Model Fitting}
\label{sec:nonlinear.psf.model}

% gaussSigma = MOMENTS_GAUSS_SIGMA from recipe (initially)
% gaussSigma = Sigma used for window function, $\sigma_w$

% fitRadius, apRadius = (fitScale, apScale) * gaussSigma
% fitScale = 7.0
% apScale = 4.5

\textadd{At this point, we have a PSF model for the image, we have an
assessment of the size (PSF-like, extended, or cosmic-ray) for each
object, and we have fitted the PSF model for the normalization to each
source (Section~\ref{sec:ensemble.fitting}).  However, the positions
for the sources have been fixed to the position determined from the
peak detection stage (Section~\ref{sec:peaks}) or the centroid from
the analysis of the moments (Section~\ref{sec:moments}).  A better
position, and thus a better normalization, can be determined by
simultaneously fitting for all three parameters.  We therefore go
through the image and re-fit the PSF model to each source
one-at-a-time with all other sources subtracted based on the earlier
fit.}

\textmod{This re-fitting analysis is performed for all of the sources
  with signal-to-noise ratio greater than a user-defined limit.  In
  the PV3 analysis of the $3\pi$ survey data, this limit was set to a
  signal-to-noise ratio of 20.0 for the \ippstage{chip} and
  \ippstage{stack} analysis stages.  In these fits, the dependent
  parameters are fixed by the PSF model and only the 3 independent
  source model parameters (position in $X$ and $Y$ and flux
  normalization) are allowed to vary in the fit.  Note that we do {\em
    not} allow the local sky to be fitted as a free parameters.  Since
  we have subtracted a model for the background, allowing the sky to
  be again at this stage is redundant.  In fact, in our testing, we
  found that allowing the sky to float resulted in higher scatter for
  the flux normalizations.  For the non-linear fitting,
\ippprog{psphot} again uses the Levenberg-Marquardt technique.}  The
sources are fitted in their S/N order, starting with the brightest and
working down to the user-specified limit, with the other sources
subtracted as discussed above.  All sources for which a non-linear PSF
model has been attempted have the flag bit
\code{PM_SOURCE_MODE_FITTED} set, regardless of the quality of that
fit.

Since the PSF model describes the variation of the PSF across the
image, the parameters used to fit a specific object are drawn from the
model at the position corresponding to the object centroid.
Occasionally, a PSF model for an image may not be well determined in
all regions of the image.  For example, not enough bright stars were
available across the full range of the image to model the PSF and the
resulting fitted parameters yield non-sensical solutions in areas
where detected (fainter) sources are found.  In such cases, the PSF
fitting is skipped and the flag bit \code{PM_SOURCE_MODE_BADPSF} is set.

For the PSF model fitting, only pixels within a circular aperture
scaled based on the seeing are used.  The radius of the circular
aperture is set to be a fixed multiple (\code{PSF_FIT_RADIUS_SCALE})
of $\sigma_w$, the width of the Gaussian window function determined
based on the analysis of the second moments (see
Section~\ref{sec:moments}).  For the PV3 $3\pi$ analysis, the PSF fit
window radius is $7 \times \sigma_w$.

%% Once a solution has been achieved for a source, \ippprog{psphot} attempts to
%% judge the quality of the PSF model as a representation of the source
%% 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 source is well represented by the PSF, but
%% large if the PSF is not a good representation of the source flux.  The
%% model quality is judged by the change in the two shape parameters
%% which represent the 2D size of the source.  For the case of the
%% elliptical Gaussian, these two parameters are $\sigma_x$ and
%% $\sigma_y$.  For a generic model, the shape parameters may be defined
%% differently, but there should always be two parameters which scale the
%% source size in two dimensions.  Currently, \ippprog{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 source in question
%% and the value of the shape parameter itself.  The expected standard
%% deviation on the shape parameter is, eg, $\sigma_x / {\rm S/N}$.  If
%% the source is well-represented by the PSF, then the shape parameter
%% values should be close to their minimization value.  We can thus ask,
%% for each source, given the measured amplitude of the Gauss-Newton
%% step, how many standard deviations from the expected value (of 0.0) is
%% this particular value?  Sources for which the variation in the shape
%% parameters is a large positive number of standard deviations are
%% likely to be better represented by a larger flux distribution than the
%% PSF (eg, a Galaxy or Comet, etc).  Sources for which the variation in
%% the shape parameters is a large negative number of standard deviations
%% are likely to be better represented by a smaller flux distribution
%% than the PSF (ie, a cosmic ray or other defect).  A user-defined
%% number of standard deviations is used to select these two cases, and
%% to flag the source as a likely galaxy (really meaning `extended') or
%% as a likely defect.  

After the PSF model is fitted to each object, \ippprog{psphot} makes an
assessment of the quality of the PSF fits.  First, it checks that the
non-linear fitting process has converged with a valid fit.  The fit
for an object can fail if there are too few valid pixels, due to
masking or proximity to an edge, or if the parameters are driven to
extreme values which are not permitted.  In addition, it is possible
for the peak finding algorithm to identify peaks in locations which
are not actually a normal peak.  Some of these cases are in the edges
of saturated, bleeding columns from bright stars, in the nearly flat
halos of very bright stars, and so on.  In these cases, a local peak
exists, with a lower nearby sky region.  However, the fitted PSF model
cannot converge on the peak because it is very poorly defined (perhaps
only existing in the smoothed image).  In these cases, \ippprog{psphot}
flags the object with the bad bit \code{PM_SOURCE_MODE_FAIL}.  It is
also possible in this type of case for the fit to result in a very low
or negative value for the flux normalization parameter.  Sources for
which the peak is less than 0.02 are also marked as failing the
non-linear PSF fit (\code{PM_SOURCE_MODE_FAIL}).

Poor fits are also identified by the signal-to-noise and the $\chi^2$
value of the resulting fit.  If a source has a PSF S/N ratio lower
than a user-defined cutoff (set to 2.0 for the PV3 analysis of the
$3\pi$ survey), the non-linear PSF fit will be rejected.  If the
$\chi^2$ per degree of freedom is greater than a user-defined limit
(set to 50.0 for the PV3 analysis of the $3\pi$ survey), the
non-linear PSF fit will be rejected.  These sources are marked with
the flag bit (\code{PM_SOURCE_MODE_POOR}).

Sources which pass the above tests are marked as having a valid
non-linear PSF model fit with the flag bit
\code{PM_SOURCE_MODE_NONLINEAR_FIT}.  Among these sources, those for
which the peak flux is greater than the saturation limit (see
Section~\ref{sec:image.preparation}) are marked as saturated stars
(\code{PM_SOURCE_MODE_SATSTAR}).  These model fits should be
considered with caution, but the fluxes and positions may have some
validity.

% \citep[see the discussion in][regarding the masking of saturated
% pixels]{waters2017}

As the sources are fitted to the PSF model, those which survive the
exclusion stage are subtracted from the image.  The subtraction
process modifies the image pixels (removing the fitted flux) but does
not modify the mask or the variance images.  The signal-to-noise ratio
in the image after subtraction represents the significance of the
remaining flux.  If the subtractions are sufficiently accurate models
of the PSF flux distribution, \textmod{the remaining flux should be
  normally distributed about zero with a standard deviation of 1
  $\sigma$}.  In practice the cores of bright stars are poorly
represented and may have larger residual significance.

\subsubsection{Double and Blended Sources}

\textmod{In fields with high stellar density, the non-linear source fitting can
be adversely affected by close neighbors.  We implemented two
modifications of the non-linear fitting code to address this issue for
different scales to the nearby neighbors.  One version addresses the
case of nearby sources which are separately detected in the
peak-detection stage; the other version of the analysis attempts to
fit a pair of PSFs for sources which are apparently extended.  {\em
  Note that for DR1 \& DR2, neither of these options were used because
  they tended to prevent galaxies from being fitted as extended
  objects; these rules for distinguishing blended stars and galaxies will
  be revisited in future re-analyses.}  We outline the strategy below
although it was not used for these data releases.}

{\TEXTADD Sources which are blended with other sources may be fitted together as
a set of PSFs.  Blended objects are identified by first searching for
objects for which the PSF fit windows overlap.  For a group of such
neighboring objects, a contour is determined in the flux image at
$25\%$ of the peak of the brightest source in the group.  All objects
lying within this contour are treated as blends of this brightest
source.  If other objects in this group exist, the brightest object
not already assigned to a blend is used to define the contour for
blends of this next object.  All objects in the image are tested as
possible blends.  A single multi-source fit is performed on each group
of blended peaks.  Sources which are identified as members of a
blended group have the flag bit \code{PM_SOURCE_MODE_BLEND} set, while
those for which a blended PSF fit succeeds have the flag bit
\code{PM_SOURCE_MODE_BLEND_FIT} set.  For sources in groups of blended
stars, the resulting fits are evaluated independently.  Any which are
determined to be valid PSF fits are subtracted from the image and kept
for future analysis.}

{\TEXTADD Sources which are judged to be non-PSF-like are confronted with two
possible alternative choices: double-star or extended source model
(see next section).  For the double-star model, the assumption is made
that there are two neighboring PSF-like sources, but the peaks are not
resolved.  The initial guess for the two peaks is made by splitting
the flux of the single source in half and locating the two starting
peaks at +/- 2 pixels from the original peak along the direction of
the semi-major axis of the sources, as measured from the second
moments.  In order for the two-source model to be accepted, both
sources must be judged as a valid PSF source.  Otherwise, the
double-PSF model is rejected and the source is fitted with the
available non-PSF model or models.  Sources for which a double-PSF
model is fitted have the flag bit \code{PM_SOURCE_MODE_PAIR} set. }

\subsubsection{Non-PSF Sources}
\label{sec:nonlinear.galaxy.model}

\textmod{Once every source (above the S/N cutoff) has been confronted with the
PSF model, the sources which are thought to be extended (resolved) can
now be fit with an appropriate model (e.g., galaxy profile or other
likely extended shapes).}  Again, the fitting stage starts with the
brightest sources (as judged by the rough S/N measured from the
moments aperture) and working to a user defined S/N limit.

{\TEXTADD The choice of extended source model or models is set by the user for a given
analysis.  In the configuration system, the keyword \code{EXT_MODEL}
is set to the model of interest.}  All suspected extended sources are
fitted with the model, allowing all of the parameters to float.  The
initial parameter guesses are critical here to achieving convergence
on the model fits in a reasonable time.  The moments and the pixel
flux distribution are used to make the initial parameter guess.  Many
of the source parameters can be accurately guessed from the first and
second moments.  The power-law slope can be guessed by measuring the
isophotal level at two elliptical radii and comparing the ratio to
that expected.

For each type of extended source model (in fact for all source
models), a function is defined which examines the fit results and
determines if the fit can be considered as a success or a failure.  The
exact criteria for this decision depends on the details of the model,
and so this level of abstraction is needed.  For example, in some
case, the range of valid values for each of the parameters must be
considered in the fit assessment.  In other cases, we may choose to
use only the parameter errors and the fit chi-square value.

All extended source model fits which are successful are then
subtracted from the image as is done for the successful PSF model
fits.  The background flux is retained, with the result that only the
source is subtracted from the image.  At this stage, the variance
image is not modified.  

For the single exposure (\ippstage{chip}) and \ippstage{stack} image
analysis, these galaxy model fits are only used internally to generate
a clean object-subtracted residual image.  For the PV3 analysis of the
$3\pi$ survey, these model fits were saved in the output catalog
files, but not loaded to the public database.  The \code{QGAUSS}
extended source model was used for the PV3 analysis (see
Section~\ref{sec:Source.Model}).  The convolved galaxy model fits (see
Section~\ref{sec:galaxy.conv.fit}) and the forced galaxy model fits
(see Section~\ref{sec:galaxy.forced.fit}) provide more reliable and
physically-motivated galaxy models.

For the difference image analysis, a trailed object model is used for
the extended sources; these model fit parameters are passed to the
Moving Object Processing System \citep[MOPS][]{2013PASP..125..357D}.

Any source which is fitted with the extended source model has the flag field
\code{PM_SOURCE_MODE_EXTMODEL} set.

\subsection{Faint Source Analysis}
\label{sec:faint.psf.model}

After a first pass through the image, in which the brighter sources
above a high threshold level have been detected, measured, and
subtracted, \ippprog{psphot} optionally begins a second pass at the
image.  \textadd{Some of the steps described in the previous sections
  are repeated in this analysis, though with some modifications as
  discussed below.}

\textadd{The source detection steps described in
  Section~\ref{sec:source.detection} are repeated.  To start, the
  sources detected in the previous steps are subtracted from the
  image, and the variance enhanced by adding the variance predicted by
  the model to the variance image, doubling the effective variance at
  the location of previously detected sources.  As in
  Section~\ref{sec:peaks}, the image is smoothed, but in this pass it
  is convolved with the PSF determined above, not a place-holder
  Gaussian.  The new peaks are detected on the smoothed image. The
  peak detection process again uses the variance image to test the
  significance of the individual peaks.  All peaks with a significance
  greater than a user-defined minimum threshold are accepted as
  sources of potential interest.  Footprints are again generated as in
  Section~\ref{sec:footprints}.}

\textadd{Next, moments are measured as in Section~\ref{sec:moments}.
  In this pass, however, the size of the window function applied for
  the measurement of the moments is fixed at the value determined from
  the bright source analysis.  All sources, including those measured
  in the bright-source analysis (which are re-added to the image and
  their variance reset), are then simultaneously fit for their flux
  normalizations as in Section~\ref{sec:ensemble.fitting}.  In this
  step, the `best' model is used for each source, either a PSF model
  or the unconvolved extended source model determined in
  Section~\ref{sec:nonlinear.galaxy.model}.  For the newly detected
  sources, the PSF model is used, with the position set by the
  centroids.}

\textadd{After the flux-normalization is calculated, the radial
  profile is measured (Section~\ref{sec:radial.profile}) and the
  moments are used to calculate the preliminary Kron radius and flux
  (see Section~\ref{sec:kron.mags}).  These are in turn used to assess
  the source sizes as in Section~\ref{sec:source.size}.  However, the
  non-linear fitting steps for the PSF model fits
  (Section~\ref{sec:nonlinear.psf.model}) and the extended source
  model fits (Section~\ref{sec:nonlinear.galaxy.model}) are not
  performed for these faint sources.  These steps are skipped for two
  reasons: First, the non-linear fitting steps are costly in terms of
  computation time and the faint sources usually far out-number the
  brighter sources.  Second, since these are faint sources they do not
  have the signal-to-noise to constrain models with many additional
  parameters.  In addition, the positions (for PSF sources) are not
  much improved using the non-linear fitting compared with the
  non-parametric centroid measurement for these faint sources. }

The PV3 threshold for the bright source analysis is a signal-to-noise
of 20.  The flag bit \code{PM_SOURCE_MODE2_PASS1_SRC} is raised for
sources detected in this initial analysis stage.  The lower limit
cutoff for the faint source analysis in PV3 is a signal-to-noise of
5.0.

In the \ippprog{psphotStack} version of the code, the 5 filter images
are processed together.  In this case, any source which is detected in
at least two of the five filters are then also measured on the other
filter images in which it was not detected above the signal-to-noise
limit.  The position in the other stack images is fixed based on the
pixel coordinates in the images in which the source was detected.
Detection in two filters is required in order to avoid excessive
forced photometry of spurious detections.  There is an interesting
class of astronomical objects which are extremely red (\eg, brown
dwarfs and high-redshift quasars).  Such sources are expected to be
detected only in the reddest filter (\yps).  For the $3\pi$ PV3
processing, we therefore also force the photometry in all filters for
sources which are only detected in \yps.  All sources which are forced
on one image based on detections in other images have the flag bit
\code{PM_SOURCE_MODE2_MATCHED} set.

\subsection{Aperture Correction and Total Aperture Fluxes}
\label{sec:aperture.correction}

A PSF model will always fail to describe the flux of the stellar
sources at some level.  For high-precision photometry, we need to be
able to correct for the difference between the PSF model fluxes and
the total flux of the sources.  In the end, all astronomical
photometry is in some sense a relative measurement between two images.
Whether the goal is calibration of a science image taken at one
location to a standard star image at another location, or the goal is
simply the repetitive photometry of the same star at the same location
in the image, it is always necessary to compare the photometry between
two images.  If this measurement is to be consistent, then the
measurement must represent the flux of the stars in the same way
regardless of the conditions under which the images were taken, at
least within some range of normal image conditions.  So, for example,
two images with different image quality, or with different tracking
and focus errors, will have different PSF models.  To the extent the
PSF model is inaccurate, the measured flux of the same source in the
two images will be different (even assuming all other atmospheric and
instrumental effects have been corrected).  The amplitude of the error
will by determined by how inconsistently the models represent the
actual source flux.

Aperture photometry attempts to avoid these problems, but introduces
other difficulties.  In aperture photometry, if a large enough
aperture is chosen, the amount of flux which is lost will be a small
fraction of the total source flux.  Even more importantly, as the
image conditions change, the amount lost will change by an even
smaller fraction, at least for a large aperture.  
%
% This can be seen by
% the fact that the dominant variations in the image quality are in the
% focus, tracking and seeing.  All of these errors initially affect the
% cores of the stellar images, rather than the wide wings.  The wide
% wings are largely dominated by scattering in the optics and scattering
% in the atmosphere.  The amplitude and distribution of these two
% scattering functions do not change significantly or quickly for a
% single telescope and site.  
%
Aperture photometry can then be used to
correct the PSF photometry.

The difficulty for aperture photometry is the need to make an accurate
measurement of the local background for each source.  As the aperture
grows, errors in the measurement of the sky flux start to become
dominant.  If the aperture is too small, then variations in the image
quality are dominant.  The brighter is the source, the smaller is the
error introduced by the large size of the aperture.  However, the
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.  

In order to thread the needle between these effects, \ippprog{psphot}
measures the aperture photometry on a modest-sized aperture, and then
uses the PSF model to extrapolate to a large aperture.  When the PSF
fluxes are calculated, the aperture flux for the modest-sized aperture
is also determined.  The aperture is a circular aperture with radius
set to a fixed multiple (\code{PSF_APERTURE_SCALE}) of $\sigma_w$, the
width of the Gaussian window function determined based on the analysis
of the second moments (see Section~\ref{sec:moments}).  For the PV3
$3\pi$ analysis, the aperture window radius is $4.5 \times \sigma_w$,
while the large reference aperture radius is set to 25 pixels
(\code{PSF_REF_RADIUS} = 6\farcs4).  These corrected aperture
magnitudes are saved in the output catalogs as \code{AP_MAG}, the
uncorrected aperture magnitudes are saved as \code{AP_MAG_RAW}, and
the radius used to measure the raw aperture flux is saved as
\code{AP_MAG_RADIUS}.  The corresponding flux and the flux error are
saved as \code{AP_FLUX} and \code{AP_FLUX_SIG}.

% Figure 6:  ** repaired PDF text **
% pueo:/home/real/eugene/ppsim.20200407/tap_psphot_deteff.pro : all.complete
\begin{figure}[htbp]
  \begin{center}
 \includegraphics[width=\hsize,clip]{\picdir/{completion.ppsim}.pdf}
  \caption{\label{fig:complete.ppsim} Completeness as a function of
    magnitude (blue curves) for different stellar densities in
    simulated data.  The curves are labeled with the logarithm of the
    stellar density at the detection threshold of the low-density
    image.  The dotted red line shows the detection limit expected for
    the sky level and seeing.  The solid red curve shows the
    completeness estimated for the low-density image based on
    injection and recovery.}
  \end{center}
\end{figure}

With these aperture magnitudes in hand, it is now possible to make an
average correction to the PSF magnitudes to bring the PSF and aperture
magnitudes to the same system.  This correction is measured using the
same stars from which the PSF model is measured, as long as the PSF
magnitude error for the star is less than 0.03 mag.  The correction is
calculated using the weighted average of the values $m_{\rm AP} -
m_{\rm PSF}$.  Since the PSF may vary across the image, the correction
is determined as a function of position in the image.  Like the PSF
model, the 2D variations of the aperture correction may be modeled as
a polynomial or via interpolation in a grid.  For the $3\pi$ PV3
analysis, a grid with a maximum of $6\times 6$ samples per GPC1 chip
image was used.  The reported PSF magnitudes for all objects have this
aperture correction applied.  \textadd{Note that an initial aperture correction was
measured during the initial steps of the analysis before the PSF model
was chosen.  However, since the sources in the image were not yet
measured and subtracted, that aperture could be contaminated by
neighbors.  The analysis here is performed one fairly bright star at a
time with all other sources subtracted in order to minimize such contamination.}

% growth curve analysis in psphot:
% in psphotChoosePSF : call psphotMakeGrowthCurve
% in psphotMakeGrowthCurve : boolean GROWTH_FROM_SOURCES, use
%% pmGrowthCurveGenerateFromSources or
%% pmGrowthCurveGenerate (uses PSF model only)
%% GROWTH_FROM_SOURCES is set to TRUE for default recipe

%% ApTrend:
%% in psphotApResid, called by psphotReadout near the end of the
%% analysis
%% ApTrend = f(x,y) for (apMag - psfMag) for psfMagErr <= 0.03
%% apMag is growth curve corrected
%% psfMag is raw

%% raw psfMag and raw apMag are measured
%% apMag = apMagRaw + growth curve correction (from apRadius to 25 pix
%% = PSF_REF_RADIUS)
%% psfMag = psfMagRaw + aperture trend (<ap - psf> + growth curve)

% How important is this effect?  Consider a typical bright source 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 source
% 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}.

%%% \ippprog{psphot} measures the aperture correction ({\em ApResid}) for every PSF
%%% candidate source, 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 source and the average bias inherent
%%% in the sky measurement for the image.  The scatter of the
%%% PSF-candidate source 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 source.  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.

\subsection{Completeness \& Contamination}

{\TEXTADD At the end of the \ippprog{psphot} analysis of the sources in the
image, an analysis is performed to test the detection efficiency.  A
number of fake PSF sources are injected into the image and the peak
detection analysis (Section~\ref{sec:peaks}) is use to determine if
these sources would have been recovered.  The PSF model fluxes are
measured for the source which are detected.  For a given image, the
detection threshold is predicted based on the median image variance
and the seeing.  A series of brightness bins straddling the threshold
are defined and a number of sources are injected with magnitudes
corresponding to each of these bin values.  The \ippprog{psphot}
recipe value \code{EFF.NUM} specifies the number of sources in each
brightness bin (500 the PV3), and the value \code{@EFF.MAG} specifies
the bins as magnitudes above and below the threshold.  For PV3, the 13
magnitude offsets were (-2.0, -1.0, -0.5, -0.25, -0.1, -0.05, 0.0,
0.05, 0.1, 0.25, 0.5, 1.0, 2.0), providing fine sampling near the
limit, but more coarse coverage further away.  Poisson noise
appropriate to the photon counts of the injected sources are included
in the image.  Injected sources are uniformly distributed across the
image in $X$ and $Y$ pixel coordinates {\em without any consideration
  of the masked regions}.  This last point means the recovered
fraction in the bright bins can be used to test the masking fraction.}

{\TEXTADD As the stellar density increases, the completeness suffers due to
crowding and confusion.  Since the injection and recovery analysis of
the fake sources operates on the source-subtracted image and does not
attempt to fully discovery the sources, this analysis over-estimates
the completeness in crowded fields.  To explore the completeness in
crowded field images, we generate a series of simulated images using a
Gaussian PSF with FWHM = 1\arcsec\ for a range of stellar densities.
We generate fake stars with fluxes as faint as $\frac{1}{5}$ of the
flux as the low-density detection limit, with densities ranging from
\approx 14,000 stars per square degree at low-density detection limit
to \approx 4.8 million stars per square degree at the low-density
detection limit.  The latter is comparable to observed densities in
the Galactic plane.  We run the \ippprog{psphot} analysis on these
simulated images and compare the detected stars to those injected to
calculate the completeness for each image as a function of the true
magnitude of the stars.  Figure~\ref{fig:complete.ppsim} shows the measured
completeness for each of the six simulated images, labeled by the
logarithm of their faint-end stellar density. The red dashed line
shows the expected detection limit based on the background and seeing,
while the red curve shows the completeness curve calculated
automatically by \ippprog{psphot} using the injection and recovery
analysis.}

{\TEXTADD For low-density fields, the completeness function determined by
injection and recovery is similar to that measured by the simulation,
with the 50\% completeness threshold roughly 0.3 magnitudes too faint.
As the stellar density increases, the true 50\% completeness magnitude
rises relative to the value estimated by injection and recovery.}

{\TEXTADD Ideally, all sources detected by \ippprog{psphot} would correspond to
real astrophysical objects.  In reality, many sources are detected in
the images which do not correspond to real sources in the sky.  In the
very simplified simulations discussed above, which do not include
realistic detector artifacts, we find that the fraction of bogus
detections is extremely low, even at the very faint end.  In real
data, bogus detections are due to a variety of typical instrumental
features including cosmic rays, diffraction spikes, satelite tracks,
glows, non-Gaussian noise, variance mis-estimation, etc.  See paper III
for extensive discussion of instrumental artifacts in the Pan-STARRS images.}

% Figure 7: ** repaired PDF text **
% pueo:/home/real/eugene/psphot.complete.20200407/complete.sh : full.figure.all
\begin{figure}[htbp]
  \begin{center}
 \includegraphics[width=\hsize,clip]{\picdir/{psphot.complete.pv3}.pdf}
  \caption{\label{fig:complete.pv3} Completeness and bogus fraction
    as a function of magnitude for different stellar densities in real
    PS1 exposures.  Each panel represents an exposure at different
    Galactic latitudes towards anti-center, labeled by the density of
    stars at the detection limit of the low-density exposure.  In each
    panel, the completeness (compared to deep stack data) and fraction
    of false detections (bogus fraction) is shown for a series of
    cuts.  The gold curves show all detections in the exposures.  The
    dotted black curve shows the impact of cutting detections
    identified by {\tt psphot} as cosmic rays.  The blue curve
    excludes cosmic rays and detections with {\tt PSF\_QF} $< 0.95$
    while the red curve excludes cosmic rays and detections with {\tt
      PSF\_QF\_PERFECT} $< 0.95$.}
  \end{center}
\end{figure}

% Figure 8:  ** repaired PDF text **
% /data/kukui.1/eugene/psphot.examples.20190423/compare.sh : mkfigure
\begin{figure*}[htbp]
  \begin{center}
 \includegraphics[width=\hsize,clip]{\picdir/{mag.resid.psf.v1}.\plotext}
  \caption{\label{fig:mag.resid.psf} PSF Photometry demonstration.
    Panel (d) shows the difference of the measured PSF photometry for
    stars in the first image of an image sequence with constant
    pointing compared to the next 17 images, after correction for a
    relative zero point, as a function of the instrumental magnitudes
    above the detection threshold.  Gray dots are from stars for
    which both measurements have {\tt PSF\_QF} $> 0.95$, while light red
    dots have lower {\tt PSF\_QF} values.  The top three panels (a) -
    (c) show histograms in three magnitude ranges for the magnitude
    difference divided by the reported measurement error: $N\sigma =
    (m_0 - m_1) / \sqrt{\sigma_0^2 + \sigma_1^2}$.  The red curves are
    Gaussian fits to these histograms, with the measured standard
    deviations in the upper-right corners of the plots.  The magnitude
    ranges are listed in the upper-left corners of the three plots and
    the boundaries are marked as vertical red lines in the lower plot.
  }
  \end{center}
\end{figure*}

{\TEXTADD Figure~\ref{fig:complete.pv3} illustrates the completeness and bogus
detection fraction for a set of 4 real PS1 exposures from the $3\pi$
Survey.  This figure uses \ips-band exposures with Galactic longitude
roughly 200\degrees and latitudes of 0, 10, 30, 90 degrees.  We
identify the real astrophysical sources in these fields by comparing
with the deeper stack exposures and counting as real any source
detected in both \rps\ and \ips.  We correct for the masking fraction
in the exposures (which is roughly 80\%) in the case of GPC1 and plot
the completeness fraction for all detections in 0.5 magnitude wide
bins from the saturation limit to below the detection limit.  We also
show the bogus fraction, calculated as $1 - f_{\rm pure}$, where
$f_{\rm pure}$ is the ratio of real detections to all detections for
the given sample.  We then apply three cuts to remove certain kinds of
bogus sources.  First, we exclude cosmic rays identified by
\ippprog{psphot} by rejecting sources with the flag bit
\code{PM_SOURCE_MODE_CR_LIMIT} (see Section~\ref{sec:source.size}).
Next, we also remove detections with \ippmisc{PSF_QF} less than 0.95.
Because this cut removes detections with heavy masking, it exclude a
number of bogus detections due to glows and edge defects.  Finally, we
also exclude  detections with \ippmisc{PSF_QF_PERFECT} less than
0.95.  This cut removes detections on residual persistent glows and
diffraction spikes.}

% Figure 9: ** repaired PDF text **
% /data/kukui.1/eugene/psphot.examples.20190423/compare.sh
\begin{figure*}[htbp]
  \begin{center}
 \includegraphics[width=\hsize,clip]{\picdir/{mag.resid.aper.v1}.\plotext}
  \caption{\label{fig:mag.resid.aper} Aperture Photometry
    demonstration.  The plots show identical measurements to those in
    Figure~\ref{fig:mag.resid.psf}, but for aperture photometry, as discussed in
    Section~\ref{sec:aperture.correction}, rather than PSF photometry.}
  \end{center}
\end{figure*}

{\TEXTADD For the exposures at high-Galactic latitude, with a relatively low
density of sources, the cosmic rays represent a significant
contamination, as seen in the excess of bogus sources with \ips-band
magnitudes in the range 17 - 19.  These are efficiently removed with
the cosmic ray cut listed above without noticable impact on the
completeness.  The other two cuts remove significant numbers of bogus
detections, especially at the faint end, but at a significant cost in
completeness at even brighter magnitudes.  The completeness impact of
these cuts is more significant at low-Galactic latitude, likely
because the chance of having a source lie on the diffraction spikes or
persistence glows is greatly increased at higher stellar densities.
The impact of the crowding on the completeness is also clear in this dataset.}

\subsection{Stellar Photometry Example}
\label{sec:phot.example}

To illustrate the quality of the stellar photometry as measured with
PSF and aperture magnitudes, we show the results of an analysis of a
set of 18 images obtained by PS1 19 February 2010.  These images were
obtained for the stellar transit survey ``Pan-Planets''
\citep{2016AA...587A..49O} and thus target a relatively dense
Galactic plane field.  The observations were obtained with
approximately consistent pointing, reducing our sensitivity to
small-scale variations in the flat-field structures.

Figures~\ref{fig:mag.resid.psf} and ~\ref{fig:mag.resid.aper} show
comparisons of the PSF and aperture photometry measured for these 18
images.  In these figures, the photometry has been measured using the
configuration for \ippprog{psphot} as used for the full PV3
\ippstage{chip} analysis.  The first image of the sequence is compared
to the remaining 17 images.  A relative zero-point correction is
applied, measured as the median of the photometry difference for stars
with signal-to-noise greater than 50.  The combined error is reported
and used to generate the histograms shown in the figures.  From these
two figures, one can observe the trade-off between PSF and aperture
photometry.  For the brightest instrumental magnitudes, corresponding
to signal-to-noise ratios greater than roughly 300, aperture
magnitudes provide a more consistent measurement of the stellar
fluxes, while the PSF magnitudes are more reliable for fainter
sources.  Catastophic failures or extreme outliers are also reduced
for the PSF photometry.

We largely attribute the behavior on the bright end to systematic
errors in the photometry due to our inability to perfectly represent
the shape of the PSF.  The PSF of stars at the bright end will depend
on the brightness because of the ``Brighter-Fatter'' effect
\citep[][]{2014JInst...9C3048A,2015JInst..10C5032G}, in which the
charge already present in the pixels will force the newly arriving
photoeletrons to be systematically pushed away from the accumulating
stellar image, but we do not include a brightness term in our PSF
model.  Detector or electronic non-linearity may also affect the PSF
shape and thus the PSF photometry, though non-linearity will affect
the reported photometry for both PSF and aperture magnitudes.

% Figure 10: ** repaired PDF text **
% on pueo ~eugene
% /data/kukui.1/eugene/psphot.examples.20190423/compare.sh : figure.resids
\begin{figure}[t]
  \begin{center}
 \includegraphics[width=\hsize,clip]{\picdir/{bright.mag.resid}.\plotext}
  \caption{\label{fig:mag.resid.stdevs} Demonstration of photometric
    accuracy using the image sequence from
    Figure~\ref{fig:mag.resid.psf}. Using only bright stars (7 - 8
    magnitudes above the detection threshold), we calculate the
    difference between the magnitudes in the first image and the other
    17 images.  The plotted dots show for each pair the mean
    difference vs the standard deviation of the difference.  Red dots
    show the PSF magnitudes and blue dots show aperture
    magnitudes. Despite real transparency variations of 0.4 over the
    50 minutes of this sequence, magnitudes are consistent at the few
    millimagnitude level.  Aperture magnitudes have scatter in
    the 2 - 7 millimagnitude range, while the PSF magnitudes have
    scatter of 7 - 14 millimagntiudes.  
}
\end{center}
\end{figure}

We believe the observed behavior at the faint end is primarily a
result of the increased crowding.  Aperture photometry is more
adversely affected by close neighbors than PSF photometry.  Compared
to the formal errors, the faint PSF photometry is the most reliable,
with the aperture photometry degrading rapidly as the flux of the star
decreases.  

% Figure 11: ** repaired PDF text **
% on pueo ~eugene/zpts.20200406/mana.sh : go.zpt.stats i
\begin{figure*}[tb]
  \begin{center}
 \includegraphics[width=\hsize,clip]{\picdir/{zpt.mjd.v0.i}.\plotext}
  \caption{\label{fig:zpt.iband} Historical \ips-band zero points.
    Blue dots are the individual exposure zero points, corrected to
    airmass at the zenith.  Red dots are the median of zero points
    from images groups in bins of 10 nights.  The grey line is a
    spline fit to these median values.  }
\end{center}
\end{figure*}

{\TEXTADD The figures above show the relative photometric accuracy for
  observations at a consistent pointing compared to the photon
  counting statistics. A related question is to ask how consistent is
  the photometry of the very brightest stars in terms of magnitudes.
  Figure~\ref{fig:mag.resid.stdevs} shows the accuracy of the
  brightest stars in these images for both PSF and aperture
  magnitudes.  The relative zero point between the 1st image in the
  sequence and each of the remaining images was calculated and the
  standard deviations were measured using stars 7 to 8 magnitudes
  brighter than the detection threshold, for which the photon noise is
  less than 1 millimagnitude.  Significant zero-point differences
  between the images are observed, largely due to the atmospheric
  transparency variations.  Even so, the relative zero points
  calculated from the aperture magnitudes have standard deviations of
  2.4 - 7.4 millimags with a median of 3.5 millimags, while for PSF
  magnitudes, the standard deviations are in the range 6.7 - 14.2
  millimags, with a median of 9.2. }

{\TEXTADD Our ultimate ability to accurately measure the brightness of
  individual sources depends on a few factors: the accuracy of the
  flat-field response, the consistency of the flux measurement across
  the image (either due to the accuracy of the PSF model or the
  accuracy of the aperture correction), and the accuracy of our
  correction for any zero point changes.  Our ability to accurately
  measure the zero point of each exposure depends in part on the
  characteristics of the observing site.  In hazy conditions, the
  transparency of the atmosphere may vary substantially in time but be
  relatively stable across the field-of-view of the camera, as is
  shown in Figure~\ref{fig:mag.resid.stdevs}.  Conversely, thin patchy
  clouds can result in small average transparency changes but
  substantial localized variations.  If the site experiences more
  patchy clouds than smooth haze, photometric calibration will be
  difficult.  A large fraction of time with cloudless conditions will
  benefit the calibration.}

{\TEXTADD To examine the Pan-STARRS site characteristics, we extracted
  \ips\ zero points for the lifetime of the observatory (2009 June -
  2020 April), shown in Figure~\ref{fig:zpt.iband}.  These zero points
  were measured as part of the PV3 analysis of the $3\pi$ Survey, and
  from the nightly data analysis after the end of the $3\pi$ Survey,
  in both cases using the Pan-STARRS-based reference catalog.  The
  zero points vary from night-to-night and over long periods.  Over
  the 11 years of PS1 operation, the observed \ips-band zero point
  (for data in good weather, extrapolated to the zenith), has varied
  over 0.175 magnitudes (see Figure~\ref{fig:zpt.iband}).  The
  long-term variations are believed to be due mostly to dust
  accumulation on the primary mirror and occasional cleaning, though
  the effect of the atmosphere cannot be ruled out.}

{\TEXTADD Figure~\ref{fig:zpt.resid.hist} shows a log-scale histogram
  of the \ips-band zero points after subtracting a smoothly varying
  spline fit to the median of groups of 10 nights.  A Gaussian fit to
  this distribution has $\sigma = 26.6$ millimags.  If we
  alternatively subtract a median zero point for each night, the
  standard deviation is reduced to 17.6 millimags.  These values can be
  compared to the results of \cite{2012ApJ...756..158S} in which only
  photometric nights were included, yielding a standard deviation of
  9.0 millimags.  On short time scales, weather (e.g., clouds \& haze)
  causes the deviations to lower zero point values.  A small fraction
  of positive deviations also seen in Figure~\ref{fig:zpt.resid.hist}
  which are not expected from the normal effects of weather.  We
  believe these are largely due to aperture correction errors.}

% Figure 12: ** repaired PDF text **
% on pueo ~eugene/zpts.20200406/mana.sh
\begin{figure}[b]
  \begin{center}
 \includegraphics[width=\hsize,clip]{\picdir/{zptres.hist.v0.i}.\plotext}
  \caption{\label{fig:zpt.resid.hist} Historical \ips-band zero-point
    residual variations.  Log-histogram (black line) of the
    per-exposure zero points, corrected to the zenith, after
    subtracting a spline fit to the median of image groups in bins of
    10 nights.  The inset shows the core of the distribution.  In
    both, the red line is a Gaussian fit to the distribution.  The
    large negative tails are due to clouds and haze.  }
\end{center}
\end{figure}

\subsection{Basic Analysis Summary}

\textadd{This section is focused on the basic analysis of the image
  for point-source detection and measurement.  This analysis is
  applied as described to the invidual exposures in the
  \ippstage{chip}-stage analysis and the measurements are exposed in
  the public release PSPS database in the \ippdbtable{Detection}
  table.  The same analysis is applied to the individual skycells in
  the \ippstage{stack}-stage analysis and the resulting values are
  presented in the PSPS \ippdbtable{StackObjectThin} and
  \ippdbtable{StackObjectAttribute} tables, with the later presenting
  values in instrumental units and the former giving calibrated
  values.  The detection efficiency information determined from the
  injection and recovery analysis is stored in the
  \ippdbtable{ImageDetEffMeta} and \ippdbtable{StackDetEffMeta} tables
  for the \ippstage{chip} and \ippstage{stack} stage analysis.  }

\section{Extended Source Analysis}
\label{sec:extended.source}

After the initial, fast analysis of the image relying primarily on the
PSF model, a complete analysis of the extended source properties may
be performed.  For PS1 processing, this step is skipped in the nightly
(PV0) analysis of individual exposures and only performed for the
stacks in the major reprocessings.

The extended source analysis consists of the following types of
measurements: 1) an analysis of the radial profile of the surface
brightness of the source; 2) measurement of the Petrosian radius and
magnitude; 3) convolved galaxy model fits; and 4) photometry in
several fixed-sized apertures, both raw and convolved to a defined
PSF size.  \textadd{The motivation for these measurements is to
  provide options to the end users for galaxy photometry and reliable
  galaxy colors.  The photometric redshift analysis of
  \cite{2012ApJ...746..128S}, for example, uses the convolved,
  fixed-size aperture photometry.}  

%% NOTE: This is NOT true: extended source analysis applied to both
%% PASS 1 and PASS 2 sources.  
%%
%% In order for a source to be included in the extended source
%% analysis, it much have been detected in the `bright source' analysis
%% step ($S/N > 20$, Section~\ref{sec:xxxx}).  

The extended source analysis is not applied to all \textmod{objects} which may be
galaxies.  Several restrictions are possible within the software.  For
example, it is possible to limit which objects are processed by their
apparent magnitudes, by their signal-to-noise, by an indication if they
are in fact extended, by the local stellar density, or by the galactic
latitude.  Some of these selections may be defined differently for the
galaxy model fits and the Petrosian parameters.

For the $3\pi$ PV3 processing, both an apparent magnitude cut and a
Galactic latitude cut were applied.  The apparent magnitude limits for
the galaxy model fits are applied to the measured Kron magnitude and
depend on the filter as follows: (\grizy) = (21.5, 21.5, 21.5, 20.5,
19.5).  These values were chosen to have roughly similar
signal-to-noise in a typical stack image for objects with colors of
typical galaxies. The magnitude limits for the Petrosian parameters
were set to 25.0 for all filters, far below the detection limits and
effectively not limiting the analysis based on apparent magnitude. For
both galaxy model fits and Petrosian parameters, the Galactic latitude
cut 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.  See
Figure~\ref{fig:galplanecut} for an illustration of the cut used for PV3.  \textadd{The
  Galactic plane cut is made on an object-by-object basis.}  This
contour avoids the denser portions of the Galactic plane and bulge,
limiting 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 sources).

% galaxy model fits performed based on limits set in psphotChooseAnalysisOptions.c

% petrosian analysis performed on same objects as galaxy model fits
% if EXTENDED_SOURCE_PETROSIAN == TRUE (TRUE for PV3 stack - STACKPHOT)

% galaxy model fits are performed on :
% all if (PSPHOT.EXT.FIT.ALL.SOURCES == TRUE) (FALSE for PV3 stack)
%   (even so, if density is higher than PSPHOT.EXT.FIT.ALL.THRESH, skip)

% only extended sources (based on EXT.NSIGMA) if EXT.NSIGMA.LIMIT.USE
% == TRUE (FALSE for PV3 stacks)

% fit sources / measure petrosian to fixed flux limit if limits are
% defined (they are for PV3)

% mag limits by filter, e.g., : petro 25, extfit 21.5
% are translated to flux in counts and compared to Kron flux

% SN limit is used only if fixed flux limits are not set
% SN limit set to EXTENDED_SOURCE_SN_LIM (10.0 for PV3)
% S/N limit for Kron flux, 

% S/N lim values set to 0.0 for all models in PV3

% galaxy coordinate limits:
% if |b| > 20.0 + 15.0 exp(-long^2 / (2 * 50^2))

\subsection{Radial Profiles}
\label{sec:radial.profile.v2}

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 described with a single elliptical
contour and radial profile.  

% Figure 13
% uses plots.sh in this directory
\begin{figure}[b]
 \begin{center}
 \includegraphics[width=\hsize,clip]{\picdir/galplanecut.pdf}
  \caption{\label{fig:galplanecut} Illustration of the Galactic Plane
    cut used for PV3, in Galactic coordinates.  Objects within the red
    contours are skipped for galaxy model fits and Petrosian parameters.}
  \end{center}
\end{figure}

% Figure 14  ** repaired PDF text **
% on pueo ~eugene/sdss.psphot.2020414/mana.sh : go.figure
\begin{figure*}[htbp]
 \begin{center}
 \includegraphics[width=\hsize,clip]{\picdir/{petrosians.mags}.pdf}
  \caption{\label{fig:petrosians} Comparison of PS1 ({\tt psphot}) and
    SDSS Petrosian parameters for objects identified as galaxies by
    SDSS.  Panel (a) shows the difference in the measured Petrosian
    magnitudes as a function of the Petrosian magnitude.  Panel (b)
    shows the magnitude difference as a function of the measured
    difference in the Petrosian radius.  }
  \end{center}
\end{figure*}

In order to facilitate the Petrosian photometry analysis below, \ippprog{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 then 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.
If fewer than 4 radial surface-brightness values are available for the
analysis, the source is skipped and the flag bit
\code{PM_SOURCE_MODE2_ECONTOUR_FEW_PTS} is set.  Some apparently
extended sources are in fact bright stars with central saturation.
These sources show up in this analysis as having many NAN-valued
pixels in the central regions.  During the radial profile analysis,
such sources are flagged with the bit
\code{PM_SOURCE_MODE2_RADBIN_NAN_CENTER} and are skipped from the rest
of the analysis.

The surface brightness profiles are then used to define the azimuthal
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 FITS table as a
vector (\code{PROF_SB}).  The flux at each radial position and the
fill-factor (fraction of pixels used to the total possible) are 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
output file FITS header (\code{RMIN_NN}, \code{RMAX_NN}).  \textadd{These
measurements are saved in the catalog FITS files generated by
\ippprog{psphot}, but they are not currently exported to the PSPS
database for easy access.}

\subsection{Petrosian Radii and Magnitudes}
\label{sec:petrosian}

\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 from the observer.  \textmod{Since the surface
brightness profile} in a resolved source is conserved as a function of
distance, 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 brightness
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
20\% of the average surface brightness at that point.  If the profile
falls below the Petrosian ratio for the first radial bin, the flag bit
\code{PM_SOURCE_MODE2_PETRO_RATIO_ZEROBIN} is set to note that the
Petrosian radius may be poorly determined.

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.  If the Petrosian ratio does not fall below
0.2 for any of the measured radii, the annulus for which the ratio
falls to the lowest (yet still significant) value.  In such a case,
the flag bit \code{PM_SOURCE_MODE2_PETRO_INSIG_RATIO} is set.  

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.  Sources for which the Petrosian
parameters are successfully measured have the flag bit
\code{PM_SOURCE_MODE_EXTENDED_STATS} set.  Sources for which the Petrosian
parameters were attempted, but for which the radial profile analysis
failed have the flag bit
\code{PM_SOURCE_MODE2_PETRO_NO_PROFILE} set.  \textadd{These measurements are
available from the PSPS \ippdbtable{StackPetrosian} table.}

Our implementation of the Petrosian apertures and fluxes is designed
to match the SDSS implementation \citep{2002AJ....123..485S} and
therefore the measured parameters should be quite comparable between
the two surveys.  Figure~\ref{fig:petrosians} compare the Petrosian
magnitudes and radii as measured by \ippprog{psphot} on the $3\pi$
Survey observations and the values measured by SDSS for the same
objects.  Objects identified by SDSS as galaxies ({\tt probPSF\_r} $<
0.5$) near the Galactic north pole ($\alpha$ = 180\degrees\ to
190\degrees, $\delta$ = 25\degrees\ to 35\degrees) are selected from
the PS1 $3\pi$ Survey dataset base on positional coincidence.  The
figure shows the difference in the $r$-band Petrosian magnitudes as a
function of the Petrosian magnitude and as a function of the
difference in the measured Petrosian radii.  Differences in the
measured magnitudes are driven by differences in the size estimates
from the two datasets and analysis methods.  The PS1 analysis tends to
find larger radii for the same objects than the SDSS analysis, with
a mean difference of 0.3 arcseconds.  The larger aperture results in
more flux captured in the aperture and thus brighter magnitudes for
the same object: the mean difference is -0.23 magnitude in the sense
of larger fluxes for the PS1 measurements.

\subsection{Convolved Galaxy Model Fits}
\label{sec:galaxy.conv.fit}

In the galaxy model fitting stage, sources which meet certain
criteria are fitted with analytical models for galaxies.  Three
traditional analytical galaxy models are implemented in \ippprog{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 S\'ersic \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 S\'ersic profile has the additional free parameter of
the S\'ersic 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.

Sources which passed the extended source restrictions described above
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}.  All
sources for which the extended source model fits were attempted the
flag bit \code{PM_SOURCE_MODE2_EXT_FITS_RUN} set.  If any of the
attempted model fits failed, then the flag bit
\code{PM_SOURCE_MODE2_EXT_FITS_FAIL} is set.  If all model fits
failed, then the flag bit \code{PM_SOURCE_MODE2_EXT_FITS_NONE} is
set.  

% \note{need a discussion of the detector saturation behavior?}

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 \cite[][Table1]{2005PASA...22..118G}.  They
quantify the relationships between the first radial moment used to
calculated a Kron Magnitude and the effective radius for different
S\'ersic index values, $n$.  Since the Exponential and DeVaucouleur
models are equivalent to S\'ersic 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, \textadd{as described for PSF sources
in Section~\ref{sec:psf.model.choice}}.  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 S\'ersic models and determining the
relationship between the central intensity and the integrated flux as
a function of the S\'ersic index.

The PSF-convolved galaxy model fitting analysis 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.

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  \\
2 \nabla   \chi^2  & = & \sum_p f_p \frac{\partial f_p}{\partial a_m} \\
\nabla^2 \chi^2  & \approx & H_{m,n} \\
2 H_{m,n}  & = & \sum_p \frac{\partial f_p}{\partial a_m} \frac{\partial f_p}{\partial a_n}
\end{eqnarray}
where we define
\begin{eqnarray}
f_p (a_m)         & = & \frac{1}{\sigma_p} (I_p - M_p \otimes \mbox{PSF}) 
\end{eqnarray}
and 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} \\
2 \nabla \chi^2    & = & -\sum_p \frac{I_p - \mathcal{M}_p}{\sigma_p} \mathcal{M}^\prime_{p,m} \\
2 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
technique 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 or the change in the $\chi^2$ value is small.

In the analysis, convolved galaxy fits, 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 \approx 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?  show timing comparisons?}

For the Exponential and DeVaucouleur fits, all parameters are fitted
in the non-linear minimization stage.  For the S\'ersic 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 S\'ersic index, the $R_{\rm eff}$ is
related to the 1st radial moment by the scale factor specified 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,
1.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$), S\'ersic 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 quadratic
interpolation to find the index $n$ which produces the locally minimum
$\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.  Sources for which a convolved galaxy
model fit was successful have the flag bit
\code{PM_SOURCE_MODE_EXTENDED_FIT} set.

% Graham & Driver : Graham A. W., Driver S. P.  2005, PASA 22, 118
% DOI: https://doi.org/10.1071/AS05001

The central pixel of the S\'ersic, DeVaucouleur, and Exponential
models requires special handling.  When comparing an analytical model
to the pixelized image recorded by a CCD, one normally treats the
value in a pixel as equivalent to the value of the model at the center
of the pixel.  However, in reality, the number of counts observed in a
pixel represents the integral of the surface brightness across the
area of the pixel.  This average will be equal to the central surface
brightness times the area of a pixel as long as the second and higher
derivatives of the analytical model are zero.  However, if the first
and second derivatives are non-zero, the curvature of the function
within the pixel will make the integral differ from the central
surface brightness times a fixed pixel area.  If the curvature of the
model function is sufficiently large, this difference will have a
significant impact on the evaluation of the model.   This situation is
particularly true for the central portion of the S\'ersic-like model
functions. 

%% this can be seen by writing the taylor expansion of the function
%% about the center of the pixel.  do this?

\textmod{In order to accurately compare the observed galaxy flux distribution
to a model, it is necessary to integrate the pixel flux for a given
set of model parameter values.  In the \ippprog{psphot}
implementation, we currently use a brute-force numerical evaluation of
the integral, dividing the central pixel into a grid of subpixels,
with the sampling set by the S\'ersic index of the model being
evaluated as $N_{\rm sub} = 2 Integer(6n / R_{\rm min})$ where $N_{\rm sub}$
is subpixel scale $n$ is the S\'ersic index and $R_{\rm min}$ is the
size of the minor axis in pixel units.  The value of $N_{\rm sub}$ is
constrained to be in the range 11 to 121, so the number of subpixels
evaluations ranges from 121 to $121^2 = 14,641$.  Faster
approximations to this analysis were explored but they resulted in
unsatisfactory results.  This is definitely an area where
\ippprog{psphot} could benefit from some of the lessons in the
literature \citep[e.g.][]{2013PASP..125..719H}.}

%% This could be done in a numerical
%% fashion, but in practice brute-force evaluation of the numerical
%% integral is computationally very expensive, requiring many evaluations
%% of the model function.  Within \ippprog{psphot}, we bypass this
%% problem by defining a set of pre-calculated images for the central 9
%% pixels (the $3 \times 3$ grid of pixels centered on the peak).  These
%% pixel images are defined at higher resolution, with 11 subpixels per
%% real CCD pixel.  The pre-calculated images are generated for a series
%% of values for the following parameters: S\'ersic index, effective
%% radius, axial ratio.  We then select the closest image to our specific
%% case, and integrate over the true sub-pixels relevant for our position
%% and model.  We have thus turned the problem from thousands of
%% evaluations of the full galaxy model to \approx 100 straight
%% additions, or up to $6 \times$ that number if we interpolate between
%% any of the parameters.

The convolved galaxy model fit results are available in one of three
PSPS database tables: \ippdbtable{StackModelFitExp},
\ippdbtable{StackModelFitDeV}, \ippdbtable{StackModelFitSer} for the
Exponential, DeVaucouleur, and S\'ersic models, respectively.

\subsection{Fixed Aperture Photometry}
\label{sec:fixed.aperture.photom}

%% some notes:
%% PV3 3pi used the STACK_THREEPI recipe for ppStack
%% this set (PSF.INPUT.MAX = 10.0)

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.  \ippprog{psphot} determines
a set of PSF-matched radial aperture flux measurements in order to
minimize the impact of the stack image quality variations.

In \ippprog{psphotStack}, the stack analysis version of \ippprog{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 fixed circular apertures are measured.  In the first
set, the fluxes in the apertures are measured using the raw stack
images.  The centers of the apertures for each source 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 source is subtracted from the image pixels for all
sources excluding the source in consideration.  The `best model' is
determined based on the minimum $\chi^2$ value for the model fits.

In addition to the raw fixed circular apertures, the stack images are
each convolved with a circular Gaussian with $\sigma$ chosen to yield
an image with a typical FWHM of 6 pixels (1.5\arcsec).  The full set
of circular apertures are again measured on these convolved images.
Again, the best source models are subtracted from the image for
sources 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 pixels
(2\arcsec).  \textadd{Note that we do not attempt to restrict the
stack image quality to match these convolution targets.  If the
stack has an effective FWHM larger than either the 6 or 8 pixel
targets, the convolution does not take place and the resulting
analysis is performed on the raw stack.  In such cases, the smallest
apertures are sensitive to seeing variations and should be avoided.
Of the individual images, the fraction with FWHM large than 8
pixels is (\grizy) = (9.8\%, 5.1\%, 4.9\%, 3.4\%, 3.7\%), and the
fraction of stacks with an effective FWHM larger than this limit will
be much smaller.}

%% filter : fraction with FWHM < the given value:
%%   2 arcsec (7.78 pix) : 10 pixels
%% g : 0.902134021127 0.976721262999
%% r : 0.949197642856 0.988031447337
%% i : 0.950975372093 0.987227497659
%% z : 0.965690586534 0.992015231573
%% y : 0.963653738055 0.990022484542

For the PV3 analysis of the $3\pi$ survey data, the fluxes are
measured for a set of up to 9 circular apertures with sizes chosen to
match the similar circular apertures measured by the SDSS analysis.
These apertures have radii of (4.16, 7.04, 12.0, 18.56, 29.76, 45.68,
72.80, 112.80, 176.88) pixels = (1.04, 1.76, 3.00, 4.64, 7.44, 11.42,
18.20, 28.20, 44.22) arcseconds.  If the object is too faint, the
larger apertures will be largely noise and the computation is
wasteful.  We only calculate the circular apertures out to the second
aperture larger than the ``sky radius'' (defined in
Section~\ref{sec:radial.profile}), but we calculate photometry for at
least the smallest 4 apertures.  Sources for which photometry in these
fixed aperture are calculated have the flag bit
\code{PM_SOURCE_MODE_RADIAL_FLUX} set.  \textadd{Although these aperture are
chosen to match the SDSS apertures, the SDSS images are measured on
unconvolved images.  Since the median seeing for the SDSS images is
$\sim 1.4$ arcseconds in $r$-band \citep{2007ApJS..172..634A}, our 1.5
arcsecond aperture photometry should generally compare well to the
SDSS aperture magnitudes.}

\textadd{The measurements described in this subsection are presented
  in the PSPS database (Paper VI) in the 
  \ippdbtable{StackApFlxExGalUnc}, \ippdbtable{StackApFlxExGalCon6},
 \ippdbtable{StackApFlxExGalCon8}, and \ippdbtable{StackApFlx} tables.
 The first three tables present measurements for all apertures from
 the unconvolved, 6, and 8-pixel FWHM convolved images (respectively)
 while the last table presents a subset of the radii from all three
 sets of measurements joined together for ease of access.}

% \note{test SDSS radial apertures?}

% at least out to aperture # RADIAL_AP_MIN (= 4), but no further than
% the aperture with radius > skyRadius
% last bin is first with inner radius >= skyRadius

\subsection{Galaxy Model Simulations}

To test the galaxy model analysis, we have generated a series of
simulated images containing both stars and galaxies on which we have
run the \ippprog{psphot} PSF-convolved galaxy model fitting analysis.
The images generated for this analysis have dimensions of $4000 \times
4000$ pixels with a spatial scale of 0.25 arcseconds per pixel.  The
images are generated using an effective exposure time of 30 seconds,
with zero points matching the PS1 \rps-filter, and a realistic sky
brightness of 20.86 magnitudes per square arcsecond.  The stars are
injected into these images with fluxes drawn from a realistic stellar
luminosity function and random spatial locations.  For each image, the
same underlying simulated stellar population was used.  Galaxies are
injected into the image with positions on a regularly spaced grid with
separation of 120 pixels.  A total of 1089 galaxies are injected in
each image, with all galaxies in an image having the same total
magnitude.  Galaxies were injected for 51 magnitude bins, ranging from
17.0 to 22.0.  Given the parameters of the image, these values span
the range from high signal-to-noise ($> 100$) to undetectable (S/N $<
1$).  Note that we are not attempting to represent a realistic
astronomical distribution of galaxies, but rather attempting to
understand our ability to recover galaxies of a certain type, size,
and shape in a given image.

The galaxies are injected using Exponential and DeVaucouleur profiles
in separate simulation runs.  The major axis values are randomly
distributed between 1 and 10 pixels (0.25 - 2.5 arcseconds) while the
aspect ratios are randomly chosen in a range from 0.25 to 1.0.  The
position angles are set by the sequence in the image and allowed to
vary from 0 to 180 degrees.  The images are then convolved with a
circular PSF model using the \code{PS1_V1} profile (\textadd{FWHM = 1.0
arcseconds, }$\kappa = 0.2$) and noise is added using Poisson
statistics for the detected photons.

For the figures below, we present results as a function of the (input)
instrumental magnitude of the galaxy minus the instrumental magnitude
corresponding to the stellar $5 \sigma$ detection limit.  We make the
simplifying assumption that the stellar detection threshold
encapsulates enough information about the sensitivity of the images
that this magnitude difference may be used to compare the results
shown here to images with other depths.  Thus this and subsequent
figures may be compared with the reported detection limits from the
PS1 $3\pi$ survey.  Note for reference that the typical stellar
detection limits in the PS1 $3\pi$ stack images (Paper I) are (\grizy) = (23.3,
23.2, 23.1, 22.3, 21.4).  The minimum Kron magnitudes for which galaxy
model fits were performed for the PV3 analysis
(Section~\ref{sec:extended.source}) thus correspond to -1.6 to -1.8 in
these plots.

Figure~\ref{fig:galaxy.complete} shows completeness for the detection
of the Exponential and DeVaucouleur model galaxies.  This analysis
does not indicate if the galaxy was detected {\em as a galaxy} (\ie,
was the extended nature of the source sufficiently clear), only if
the source was detected by the peak-finding algorithm.  As expected,
the more compact galaxies are more likely to be detected; Exponential
profile galaxies, with a broader light distribution for the same
effective radius, are less likely to be detected for the same
magnitude than DeVaucouleur profile galaxies.  This completeness
should be compared to our earlier work \citep{2013MNRAS.435.1825M} in
which we injected a realistic population of simulated galaxies into
real PS1 images.  That work found that the 50\% completeness for the
typical galaxy was roughly 0.5 magnitude brighter than the 50\%
stellar completeness limit, somewhat fainter than the completeness
shown in Figure~\ref{fig:galaxy.complete}.  However, that previous
work did not explore the depedency of the completeness on the galaxy
size or profile.  The difference suggests that the galaxies in the
earlier work were generally compact.

% Figure 15: ** repaired PDF text **
% was /data/kukui.1/eugene/galaxies.20190425/tap_psphot_galaxies.pro : go.bigtest.ckgalaxy
% is pueo:galaxies.20190425/tap_psphot_galaxies.pro : go.bigtest.ckgalaxy
\begin{figure}[htbp]
  \begin{center}
 \includegraphics[width=\hsize,clip]{\picdir/{galaxy.exp.complete}.\plotext}
 \includegraphics[width=\hsize,clip]{\picdir/{galaxy.dev.complete}.\plotext}
  \caption{\label{fig:galaxy.complete} {\bf Panel (a):} Completeness
    curves for simulated galaxies with Exponential profiles.  {\bf
      Panel (b):} Completeness curves for simulated galaxies with
    DeVaucouleur profiles.  The curves are shown as a function of the
    difference between the injected instrumental magnitude of the
    galaxy and the magnitude corresponding to the $5\sigma$ detection
    threshold for a PSF-like source.  The black curves shows the
    compleness for all galaxies.  The three colored curves show the
    completeness for three major axis ranges. Compact galaxies are
    more likely to be detected since peaks are detected after
    convolution with the PSF. \textadd{The simulated images have seeing of
    1\arcsec, equal to pixels.}}
  \end{center}
\end{figure}

Figures~\ref{fig:exp.params} and \ref{fig:dev.params} demonstrate the
recovery of galaxy parameters in these simulations for galaxy using
the Exponential and DeVaucouleur models, respectively.  
Both figures show the reliability of the measured magnitudes, major
and minor axis sizes, and \textadd{ellipticities, which we represent as
$\frac{R_{\rm Major} - R_{\rm Minor}}{R_{\rm Major} + R_{\rm
    Minor}}$. Galaxies are grouped in 0.1 magnitude bins based
on their injected magnitudes.  
For all recovered
parameters, the standard deviation of the difference between the
measured parameter and the truth value is shown for all galaxies in
each magnitude bin, as well as for subsets in major-axis ranges.}
The mean of the difference, illustrating any biases, is also given
for all galaxies.  The comparison for the major and minor axis sizes
are shown as absolute measurements (in pixels) and as a fraction of
axis in question.

Some overall patterns can be observed.  First, the parameters measured
for the exponential profile galaxies are consistently more reliable
than those measured for the DeVaucouleur profiles.  Second, the errors
in the estimated magnitudes are primarily driven by errors in the size
measurements: If the sizes are over-estimated, the fluxes are also
over-estimated.  Finally, the magnitudes and major axes are more
accurate for the smaller galaxies, but the ellipticities are more
accurate for the larger galaxies.

% Figure 16 ** repaired PDF text **
% was /data/kukui.1/eugene/galaxies.20190425/tap_psphot_galaxies.pro : go.bigtest.ckgalaxy
% is pueo:galaxies.20190425/tap_psphot_galaxies.pro : go.bigtest.ckgalaxy
\begin{figure*}[htbp]
  \begin{center}
 \includegraphics[width=\hsize,clip]{\picdir/{galaxy.exp.params}.\plotext}
  \caption{\label{fig:exp.params} Parameter recovery for simulated
    galaxies with Exponential profiles.  In each panel, we show
    several statistics for the difference between the truth parameter
    and the measured value for galaxies in a series of brightness
    bins, ranging from 5 magnitudes to 0.5 magnitudes brighter than
    the PSF $5-\sigma$ detection threshold.  Statistics, as shown in
    the legend, include the average difference, the standard deviation
    for all galaxies in the bin, and the standard deviations for 3
    major axis ranges.  Panel (a) through (d) show in order the
    recovered magnitudes, ellipticity, major axis size; minor axis
    size.  Panels (e) and (f) show statistics for the major and minor
    axis difference as a fraction of the truth value.  The average
    magnitude difference is plotted so that positive numbers mean the
    fitted flux is brighter than the injected flux.  }
  \end{center}
\end{figure*}

% Figure 17 ** repaired PDF text **
% was /data/kukui.1/eugene/galaxies.20190425/tap_psphot_galaxies.pro : go.bigtest.ckgalaxy
% is pueo:galaxies.20190425/tap_psphot_galaxies.pro : go.bigtest.ckgalaxy
\begin{figure*}[htbp]
  \begin{center}
 \includegraphics[width=\hsize,clip]{\picdir/{galaxy.dev.params}.\plotext}
  \caption{\label{fig:dev.params} Parameter recovery for simulated
    galaxies with DeVaucouleur profiles.  See
    Figure~\ref{fig:exp.params} for complete explanation.} 
  \end{center}
\end{figure*}

\section{Forced Warp Analysis}
\label{sec:psf.forced.fit}

Traditionally, projects which use multiple exposures to increase the
depth and sensitivity of the observations have generated something
equivalent to the stack images produced by the IPP analysis,
\textadd{as done for example by the Canada-France-Hawaii Telescope (CFHT) Legacy Survey
  \citep{2006ApJ...647..116H} or the Cosmic Evolution Survey
  \citep[COSMOS][]{2007ApJS..172...99C}}.  In theory, the photometry
of the stack images produces the ``best'' photometry catalog, with
best sensitivity and the best data quality at all magnitudes
\citep[see e.g., the discussion of]{2017ApJ...836..187Z}.  In
practice, these images have some significant limitations due to the
difficulty of modeling the PSF variations.  This difficulty is
particularly severe for the Pan-STARRS $3\pi$ survey stacks due to the
combination of the substantial mask fraction of the individual input
exposures, the large intrinsic image quality variations within a
single exposure, and the wide range of image quality conditions under
which data were obtained and used to generate the $3\pi$ PV3 stacks.

% CFHTLS release doc:
% http://www.cfht.hawaii.edu/Science/CFHLS/T0007/CFHTLS_T0007-TechnicalDocumentation.pdf

For any specific stack, the point spread function at a particular
location is the result of the combination of the point spread
functions for those individual exposures which went into the stack at
that point.  Because of the high mask fraction, the exposures which
contributed to pixels at one location may be somewhat different just a
few tens of pixels away.  In the end, the stack images have
a effective point spread function which is not just variable, but
changing significantly on small scales in a highly textured fashion.

Any measurement which relies on a good knowledge of the PSF at the
location of an object needs to determine the PSF variations present in
the stack image or the measurement will be somewhat degraded.  The
highly textured PSF variations make this a very challenging problem:
not only would such a PSF model require an unusually fine-grained PSF
model, there would likely not be enough PSF stars in a given stack
image to determine the model at the resolution required.  The IPP
photometry analysis code uses a PSF model with 2D variations using a
grid of at most $6\times 6$ samples per skycell, a number reasonably
well-matched to the density of stars at most moderate Galactic
latitudes.  This scale is far too large to track the fine-grained
changes apparent in the stack images.

As a result, PSF photometry as well as convolved galaxy models in the
stack are degraded by the PSF variations.  Aperture-like measurements
are in general not as affected by the PSF variations, as long as the
aperture in question is large compared to the FWHM of the PSF.

%% The IPP team initially explored the option of convolving each input
%% warp to a single target PSF chosen to match the worst of the input
%% images for a given stack.  

The IPP analysis solves this problem by starting with the sources
detected in the stack images and performing \textmod{a highly
  constrained analysis on the individual warp images used to generate
  the stack, and then combining the resulting measurements to
  determine a high-quality average value.  We consider all of these
  measurements to be ``forced'' because of the strong prior
  constraints from the stack.}  This analysis is performed using the
\ippprog{psphotFullForce} variant of \ippprog{psphot}.

\subsection{Forced PSF Photometry}

\textmod{PSF photometry is measured for all objects detected in the
  stack using the positions determined by the stack photometry
  analysis.}  Candidate PSF stars are pre-identified as those stars
used to generate the PSF model in the stack photometry analysis.  A
PSF model is generated for each input warp image based on those stars;
PSF stars which are excessively masked on a particular image are not
used to model the PSF.  The \textadd{normalization of the} PSF model
is fitted to all of the known source positions in the warp images to
determine the PSF fluxes. This measurement is performed simultaneously
for all sources in the image at once using the method described above
(Section~\ref{sec:ensemble.fitting}).

\textmod{Aperture fluxes, Kron fluxes}, and moments are also measured
at this stage for each warp.  \textmod{For the Kron fluxes, the radii
  are fixed to the value determined in the analysis of the stack.
  Fluxes are also measured in 3 of the fixed apertures discussed in
  Section~\ref{sec:fixed.aperture.photom}: those with 3.00, 4.64,
  and 7.44 arcsecond radii.}
  Note that the flux measurement for a faint,
but significant, source from the stack image may be at a low
significance (less than the $5\sigma$ criterion used when the
photometry is not run in this forced mode) in any individual warp
image; the measured flux may even be negative due to statistical
fluctuations.  When combined together, these low-significance
measurements result in a significant measurement as the
signal-to-noise increases with the combination of more data.

Individual warp images are processed independently with separate
executions of the \ippprog{psphotFullForce} program.  Sources which
are loaded by \ippprog{psphotFullForce} for analysis are marked with
the flag bit \code{PM_SOURCE_MODE_EXTERNAL}.  This bit is also used to
mark user-supplied sources loaded for analysis by the regular version
of \ippprog{psphot}.  Once all of the forced photometry measurements
(for point sources as well as galaxies, discussed below) are completed
for all of the warps which contributed to a stack image, the
measurements are combined together by other portions of the IPP
system.  The PSF photometry measurements are combined in the context
of the DVO database system \citep{magnier2017.datasystem}, including
recalibration of the zero points for the individual warp.  \textadd{These
measurements for each warp are available from the PSPS database
\ippdbtable{ForcedWarpMeasurement} and \ippdbtable{ForcedWarpExtended}
tables, the latter containing the three fixed-aperture fluxes.  The
average values calculated over the warps are found in the
\ippdbtable{ForcedMeanObject} tables.}

% Figure 18 ** repaired PDF text **
% on pueo ~eugene/sdss.psphot.2020414/photcompare.sh : go.figure
\begin{figure}[htbp]
 \begin{center}
 \includegraphics[width=\hsize,clip]{\picdir/{compare.mags}.pdf}
  \caption{\label{fig:compare.mags} Comparison of {\tt psphot} average
    chip photometry, average forced-warp photometry, and stack
    photometry from $3\pi$ Survey data to average forced-warp
    photometry from the Pan-STARRS\,1 Medium-Deep Survey field MD06
    At bright magnitudes, average chip photometry is the most
    accurate while the stack photometry is degraded by the
    highly-texturd PSF.  At faint magnitudes, average chip magnitudes
    are biased to artifically bright values.}
  \end{center}
\end{figure}

{\TEXTADD With the inclusion of the forced-warp photometry, we have three
distinct methods for measuring the PSF photometry of stars in the
Pan-STARRS survey data: the average of the \ippstage{chip}-stage
photometry from the individual exposures; the measurement from the
stacks, and the average of the forced-warp photometry described here.
It is worth considering which of these should be used in which
circumstance.  Figure~\ref{fig:compare.mags} shows a comparison of
these three different methods to deeper data from the Medium Deep
Survey observations (MD06 field).  Our conclusion from this and other
analysis is that the average \ippstage{chip}-stage photometry is the
best (most accurate) measurement for brighter objects, where the
signal-to-noise is roughly 10 or more.  This is the photometry source
which was used for the global photometry solution discussed by
\cite{2012ApJ...756..158S} and used in the overall calibration (see
Paper V).}

{\TEXTADD As can be clearly seen in the figure, the average from the forced-warp
photometry is slightly worse than the chip photometry, while the stack
PSF photometry is significantly degraded.  We attribute the latter
effect to the highly-textured PSF observed in the stack images due to
the combination of variable PSFs in each exposure and significant
masking fraction in the PS1 camera.  At the faint end, the chip
photometry is significantly worse that both average warp and stack
photometry.  First, in order to have a measurement, a source must be
detected above the detection threshold in at least one of the
exposures, limiting the depth possible of the average chip
photometry. Second, at the faint end, only bright fluctuations will be
detected, resulting in a bright bias. This latter effect is clearly
seen in Figure~\ref{fig:compare.mags} as the average chip magnitudes
diverge from the deeper Medium Deep photometry measurements.  As has
been noted elsewhere \citep{2018ApJS..234....1B}, the warp and stack
photometry is also degraded for objects which have significant proper
motion over the course of the $3\pi$ Survey since the position is held
constant for all epochs, while the average chip photometry is
calculated on detections which are cross-matched in the database.
Thus, warp and stack photometry should be avoided for sources with
proper motion greater than roughly 100 milliarcseconds per year.}

\subsection{Forced Galaxy Models}
\label{sec:galaxy.forced.fit}

The convolved galaxy models are also re-measured on the warp images by
the \ippprog{psphotFullForce} analysis.  In this analysis, the galaxy
models determined from the stack image analysis are used to seed the
analysis in the individual warp images.  The motivation of this
analysis is the same as the forced PSF photometry: the PSF of the
stack image is poorly determined due to the masking and PSF variations
in the inputs.  Without a good PSF model, the PSF-convolved galaxy
models are of limited accuracy.

In the forced galaxy model analysis, we assume that the galaxy
position and position angle, along with the S\'ersic index if
appropriate, have been sufficiently well determined in the analysis of
the stack image.  In this case, the goal is to determine the best
values for the major and minor axis of the elliptical contour and at
the same time the best normalization corresponding to the best
elliptical shape, and thus the best galaxy magnitude value.
\textadd{This technique is similar to the joint fitting of multiple
  exposures performed by the CFHT Lensing Survey team \citep{2013MNRAS.429.2858M}.}

For each warp image, the stack values for the major and minor axis are
used as the center of a grid search of the major and minor axis
parameter values.  The grid spacing is defined as a function of the
signal-to-noise of the galaxy in the stack image so that bright
galaxies are measured with a much finer grid spacing than faint
galaxies.  For the PV3 $3\pi$ analysis, a $5 \times 5$ grid was used;
values in both the major and minor axis directions of ($1 -
\frac{3.0}{S/N}$, $1 - \frac{1.5}{S/N}$, 1.0, $1 + \frac{1.5}{S/N}$,
$1 + \frac{3.0}{S/N}$) times the dimension are tested.  For each grid
point, the major and minor axis values at that point are used to
generate the model.  The model is then convolved with the PSF model
for the warp image at that point.  The resulting convolved model is
then compared to the warp pixel data values and the best fit
normalization value is determined.  The integrated flux, flux error,
and the $\chi^2$ value for each grid point are recorded.
\textadd{This analysis is performed on the warp images one object at a
  time with the other objects in the image subtracted according to
  the best model currently known.}

For a given galaxy, the result is a collection of $\chi^2$ values,
fluxes, and flux errors for each of the grid points spanning all
warp images.  A single $\chi^2$ grid can then be made by
combining each grid point across the inputs.  The combined $\chi^2$
for a single grid point is simply the sum of all $\chi^2$ values at
that point.  If, for a single warp image, the galaxy model
is excessively masked, then that image will be dropped for all grid
points for that galaxy.  The reduced $\chi^2$ values can be determined
by tracking the total number of pixels used across all inputs to
generate the combined $\chi^2$ values.  From the combined grid of
$\chi^2$ values, the point in the grid with the minimum $\chi^2$ is
found.  Quadratic interpolation is used to determine the major, minor
axis values for the interpolated minimum $\chi^2$ value.  The errors
on these two parameters is then found by determining the contour at
which the $\chi^2$ increases by 1.

In this way, the forced galaxy model analysis uses the PSF information
from each warp image to determine a best set of convolved galaxy
models for each galaxy model measured for the stack image.  The
results of these galaxy model fits are available from the PSPS
database \ippdbtable{ForcedGalaxyShape} table.

\subsection{Galaxy Lensing Parameters}
\label{sec:lensing.params}

Weak-lensing studies frequently use non-parametric measurements of the
ellipticities of galaxies to quantify the strength of gravitational
lensing, and thus directly measure mass distributions in the Universe.
The classic approach was originally described by
\cite[KSB]{1995ApJ...449..460K} and applied to a set of deep HST
observations.  The details of the technique were further refined by
\cite[HFK]{1998ApJ...504..636H}; in the discussion below we primarily use
their notation, though we explicitly cast their integrals as sums over
discrete pixels.

The KSB-style analysis of object ellipticities has also been used by
several authors to search for marginally-resolved binary stars
in wide-field imaging data.  The use of the lensing statistics for
this application was described by \cite{2005ApJ...626.1070H} in the
context of vetting planet transit events in data from the Optical
Gravitational Lensing Experiment (OGLE).  \cite{2013ApJS..206...18T} 
applied the techinique to PTF data to search for binary stars and
\cite{2017MNRAS.468.3499D} used the same technique to search for
binary companions to known ultracool dwarfs using Pan-STARRS $3\pi$
data.  The work by \cite{2017MNRAS.468.3499D} used images and their
own analysis of the pixels with the program Sextractor
\citep{sextractor}.

For the Pan-STARRS $3\pi$ PV3 analysis, we have measured the full set
of KSB lensing parameters for all objects with measured second moments
(i.e., excluding saturated stars, suspected cosmic rays, and other
likely defects) of the data to enable both lensing studies and binary
/ multiple star searches.  Here we describe the measurements as
performed within \ippprog{psphot}, reviewing the mathematical
framework as described by \cite{1995ApJ...449..460K} and
\cite{1998ApJ...504..636H}.  \textadd{Just like the forced PSF
  photometry and the constrained galaxy models above, this analysis is
  performed by measuring the KSB lensing parameters on the individual
  warp images and averaging these measurements for each object.}

The goal of the KSB technique is to measure the intrinsic ellipticity
of objects (i.e., galaxies, in the case of weak lensing studies) as
would be observed sky on the without instrumental effects and to
determine the impact weak graviational lensing would have on the
observed shapes, after correction for the instrumental effects.  The
analysis starts with the observed ellipticity of objects as represented
by the two polarization components derived from the second moments
(see Section~\ref{sec:moments}):
\begin{eqnarray}
\label{eqn:polarization}
  e_1 = \frac{M_{xx} - M_{yy}}{M_{xx} + M_{yy}} \\
  e_2 = \frac{2 M_{xy}}{M_{xx} + M_{yy}}. \\
\end{eqnarray}
These two quantities have values which range from -1 to +1, with a
circularly-symmetric object having $(e_1,e_2) = 0.0,0.0$.  The two
polarization components vary sinusoidally with twice the position
angle of the object: an elongated object aligned with the $x$-pixel
axis will have positive values of $e_1$ and $e_2$ values near zero,
while the same object aligned with the $y$-pixel axis will negative
$e_1$ values.  An object with a position angle on the 45\degree\ lines
between the pixel axes will have large positive or negative values of
$e_2$ and low absolute values of $e_1$.

Note that in our analysis of the second moments, we are applying a
Gaussian window function to down-weight the noise contributions from
pixels at high radii and low flux (see Section~\ref{sec:moments}).
This type of window function is also assumed in the KSB formalism, and
is represented in the equations below as $W$.

The measured ellipticity of an object observed in a real instrument
will be affected by the point spread function of the instrument.  To
first order, the effect on the polarization components can be
described as a combination of the circularly symmetric seeing disc,
which smears the observed shapes (driving $e_1,e_2$ to low absolute
values) and the shearing effect of the anisotropic component of the
PSF, in which the observed shape is stretched in one direction
relative to the others (driving $e_1,e_2$ to larger absolute values).

KSB and HFK quantify the change in the observed polarization due to
the smearing effect of the PSF with
\begin{equation}
  \delta e^{\rm sm}_\alpha = P^{\rm sm}_{\alpha, \beta} p_{\beta}
\end{equation}
$p_\beta$ is a measurement of the
anisotropy of the PSF (see below), and $P^{\rm sm}_{\alpha,\beta}$ is
the ``Smear Polarizability'' of the object, defined as  
\begin{eqnarray}
  P^{\rm sm}_{\alpha \beta} = X^{\rm sm}_{\alpha \beta} - e_\alpha e^{\rm sm}_\beta
\end{eqnarray}  
where 
\begin{eqnarray}
X^{\rm sm}_{1,1} &=& \frac{1}{T} \sum f \left[ W + 2W^\prime r^2 + W^{\prime \prime} (x^2 - y^2)^2 \right] \\
X^{\rm sm}_{1,2} &=& \frac{1}{T} \sum f \left[ 2W^{\prime\prime} (x^2 - y^2) x y \right] \\
X^{\rm sm}_{2,2} &=& \frac{1}{T} \sum f \left[ W + 2W^\prime r^2 + 4W^{\prime \prime} x^2 y^2 \right]
\end{eqnarray}
and  
\begin{eqnarray}
e^{\rm sm}_1 &=& \frac{1}{T} \sum f \left[ 2W^\prime + W^{\prime \prime} (x^2 + y^2) \right] (x^2 - y^2) \\
e^{\rm sm}_2 &=& \frac{1}{T} \sum f \left[ 2W^\prime + W^{\prime \prime} (x^2 + y^2) \right] 2 x y.
\end{eqnarray}
In these equations, $T = M_{xx} + M_{yy}$ and $W$ is the window
function applied when measuring the second moments.  The terms
$W^\prime$ and $W^{\prime \prime}$ are the derivatives of the window
function with respect to $r^2 = x^2 + y^2$.  Since the window function
is a circularly-symmetric Gaussian with width $\sigma_w$, the
derivatives are simply $W^\prime = -\frac{1}{2\sigma^2_w} W$ and
$W^{\prime \prime} = \frac{1}{4\sigma^4_w} W$.

The elements of the equations above can be written in terms of the second and higher-order
moments calculated in Section~\ref{sec:moments}:
\begin{eqnarray}
X^{\rm sm}_{1,1} &=& \frac{1}{T} \left[ 1 - \frac{R_2}{\sigma^{2}} + \frac{(M_{xxxx} - 2 M_{xxyy} + M_{yyyy})}{4 \sigma^{4}} \right] \\[0.1in]
X^{\rm sm}_{1,2} &=& \frac{1}{T} \left[ \frac{(M_{xyyy} - M_{xxxy})}{2 \sigma^{4}} \right] \\[0.1in]
X^{\rm sm}_{2,2} &=& \frac{1}{T} \left[ 1 - \frac{R_2}{\sigma^{2}} + \frac{ M_{xxyy}}{\sigma^{4}} \right]
\end{eqnarray}
and  
\begin{eqnarray}
e^{\rm sm}_1 &=& \frac{1}{T} \left[ \frac{M_{xx} - M_{yy}}{\sigma^{2}} + \frac{M_{xxxx} - M_{yyyy}}{4 \sigma^{4}} \right] \\[0.1in]
e^{\rm sm}_2 &=& \frac{1}{T} \left[ \frac{(M_{xxxy} + M_{xyyy})}{2\sigma^{4}} - \frac{2 M_{xy}}{\sigma^{2}} \right]
\end{eqnarray}
where $R_2 = M_{xx} + M_{yy}$.

KSB and HFK use the observed ellipticities of stars and the smear
polarizability of the stars to estimate the anisotropy due to the PSF:
\begin{eqnarray}
p_\alpha = \frac{e^*_{\alpha}}{P^{{\rm sm},*}_{\alpha \alpha}}
\end{eqnarray}
where the terms with the $*$ represent parameters measured on stars.

%% \begin{eqnarray}
%%   p_1 &=& M_{xx} - M_{yy} \\
%%   p_2 &=& 2 M_{xy}
%% \end{eqnarray}

Similarly, the impact of shear can be quantified by the ``Shear
Polarizabilty'' in a similar fashion:
\begin{equation}
  \delta e^{\rm sh}_\alpha = P^{\rm sh}_{\alpha, \beta} p_{\beta}
\end{equation}
where now the shear polarizability $P^{\rm sh}_{\alpha \beta}$ is
defined as
\begin{eqnarray}
  P^{\rm sh}_{\alpha \beta} = X^{\rm sh}_{\alpha \beta} - e_\alpha e^{\rm sh}_\beta
\end{eqnarray}  
where
\begin{eqnarray}
X^{\rm sh}_{1,1} &=& \frac{1}{T} \sum f \left[ 2W(x^2 + y^2) + 2W^\prime (x^2 - y^2)^2 \right] \\
X^{\rm sh}_{1,2} &=& \frac{1}{T} \sum f \left[ 4W^\prime(x^2 - y^2) x y \right] \\
X^{\rm sh}_{2,2} &=& \frac{1}{T} \sum f \left[ 2W(x^2 + y^2) + 8W^\prime x^2 y^2 \right]
\end{eqnarray}
and
\begin{eqnarray}
e^{\rm sh}_1 &=& 2 e_1 + \frac{2}{T} \sum f W^\prime (x^2 + y^2) (x^2 - y^2) \\
e^{\rm sh}_2 &=& 2 e_2 + \frac{2}{T} \sum f W^\prime (x^2 + y^2) 2 x y.
\end{eqnarray}

Re-writing in terms of the second and higher-order moments calculated
in Section~\ref{sec:moments}, we find:
\begin{eqnarray}
X^{\rm sh}_{1,1} &=& \frac{1}{T} \left[ 2 R_2 - \frac{(M_{xxxx} - 2 M_{xxyy} + M_{yyyy})}{\sigma^{2}} \right] \\
X^{\rm sh}_{1,2} &=& \frac{1}{T} \left[ \frac{2(M_{xyyy} - M_{xxxy})}{\sigma^{2}} \right] \\
X^{\rm sh}_{2,2} &=& \frac{1}{T} \left[ 2 R_2 - \frac{4 M_{xxyy}}{\sigma^{2}} \right]
\end{eqnarray}
and  
\begin{eqnarray}
e^{\rm sh}_1 &=& \frac{1}{T} \left[ 2 (M_{xx} - M_{yy}) + \frac{( M_{yyyy} - M_{xxxx})}{\sigma^{2}} \right] \\
e^{\rm sh}_2 &=& \frac{1}{T} \left[ 4 M_{xy} - \frac{2 (M_{xxxy} + M_{xyyy})}{\sigma^{2}} \right] 
\end{eqnarray}

In the Pan-STARRS PV3 analysis, we have measured the elements of the
smear polarizability ($X^{\rm sm}_{\alpha \beta}$, $e^{\rm
  sm}_\alpha$) and the shear polarizability ($X^{\rm sh}_{\alpha
  \beta}$, $e^{\rm sh}_\alpha$) for all objects on each of the warp
images.  \textmod{We have also selected only the PSF stars from the images and
interpolated a smoothed version of these parameters from the PSF stars
to the locations of
all objects, using the grid described above to interpolate the 
parameters.  We then determine the ellipticities of the stars}
($e^*_1, e^*_2$) from the equivalent smooth grid.  Thus, for every
object in the $3\pi$ survey, we are able to report the PSF and object
elements of the KSB analysis.  These lensing parameters are measured
for each of the warps, and then averaged over all warps for each of
the filters.  The average values are calculated by including only
measurements from the same warp detection used in the average
photometry (nominally, the primary skycell; see Paper V, Section
5.4.4) and excluding any measurements for which the \code{PSF_QF} or
\code{PSF_QF_PERFECT} is less than 0.85.

The lensing parameters measured for individual warps are available
from the PSPS database \ippdbtable{ForcedWarpLensing} table while the
average values calculated over the warps is found in the
\ippdbtable{ForcedMeanLensing} tables.  \textadd{Although the software used
here was not involved in any of the GRavitational lEnsing Accuracy
Testing (GREAT) challenges, it is similar to the code of the EPFL\_KSB
team \citep{2015MNRAS.450.2963M} and likely to perform similarly.}

% \note{example of using the lensing elements for binaries?}

\section{Difference Image Photometry}

Among the primary science drivers for Pan-STARRS are the detection of
moving objects (e.g., asteroids) and explosive transient sources
(e.g., supernovae).  For both of these situations, difference images
are commonly used to remove the clutter of the static stars and
galaxies.  In the Pan-STARRS system, difference images are generated
using the PSF-matching technique described by
\citep[e.g.,][]{1998ApJ...503..325A}.  \textmod{The description of the
Pan-STARRS implementation is given by \cite{price2017} and uses an
implementation of cross-convolution based on the description of
\cite{2008ApJ...677..808Y}.  The analysis of the sources detected in
these difference images uses a portion of the \ippprog{psphot} code
embedded in the program, \ippprog{ppSub}, which generates those image.
Difference images are generated from three different possible image
combinations: 1) pairs of individual exposures are differenced using
the warp images; 2) warps for individual exposures
are differenced against deep stacks; 3) stacks made from multiple
exposures of the same field within a night are differenced against
deep stacks.  Note that this article is limited to the analysis of the
difference image detections, and that significant additional work is
needed to distinguish real detections from false positives, and
further to classify the detections as objects of scientific interest.
Within the Pan-STARRS science community, the Moving Object Processing
System \citep[MOPS][]{2013PASP..125..357D} is dedicated to the effort
of identifying asteroids and other solar system objects.  Multiple
teams have focused on the identification of supernovae
\citep{2014ApJ...795...44R,2015MNRAS.449..451W}, including the use of
machine-learning techniques to filter the good detections from the bad
detections.}

The analysis of the difference image follows the same basic steps as
other \ippprog{psphot} versions with some minor modifications (see
Table~\ref{tab:measurements}), as follows.  The background subtraction
is performed before the PSF matching and image subtraction is
performed.  The PSF model construction stage is not possible in the
difference image due to the lack of valid sources.  Instead, the PSF
model is generated from the positive image, after PSF-matching
but before the subtraction is performed.  Because we do not expect to
have a large number of sources, only a single source detection pass is
performed, and at the lowest signal-to-noise threshold.  Only linear
PSF model fitting is performed using the centroid determined from the
analysis of the source moments.  

For the difference images, the galaxy model analysis is not relevant.
In a properly-constructed difference image, galaxies are unlikely to
remain behind as significant sources.  Most real sources in the
difference image will be PSF-like and will consist of photometrically
variable sources (flare stars, supernovae, etc) or astrometrically
variable sources (high-proper motion stars or solar-system bodies).
There are three likely classes of sources which will not be well
represented by the PSF model, as discussed below.

Fast-moving solar-system objects will appear as short streaks.  For
example, a fast solar system object may have an apparent rate of 0.5
degrees per hour, translating to 15 arcseconds in a 30 second
exposure.  Even a main belt asteroid at roughly 1 AU has reflex motion
of approximately 1 degree per day, equivalent to 1.25 arcsec in a 30
second exposure, and may be noticeably smeared and non-PSF-like.  In
\ippprog{psphot}, we use a trailed-star model to characterize these
types of sources.  This model is fitted in the same portion of the
code which performs the unconvolved galaxy model analysis.

% \note{describe the trailed analytical model}.

In some cases, the stars in the two images may be somewhat offset.
For specific stars, this offset may be due to differential chromatic
aberration from the atmosphere or the optics, or from modest proper
motion.  If the astrometric solution for one of the two images is
insufficiently accurate, all stars in large portions of the images may
be noticeably displaced.  In both of these situations, the stars will
appear as PSF dipoles in the difference images.  The positive and the
negative images will have stellar profiles, but they will be offset
and will not subtract well.  The two components may not have the same
amplitude.  In theory, a PSF-dipole model could be used to fit these
types of sources, with free parameters of the two centroids and the
two fluxes.  In practice in \ippprog{psphot}, we use a number of
non-parametric pixel-level statistics in an attempt to detect these
cases.

For the difference images, we measure the following quantities for
each of the detections, using only pixels within the photometry
aperture.  First, we count the number of masked pixels (\code{nMask}),
the number of pixels with positive flux (\code{nGood}), and the number
of pixels with negative flux (\code{nBad}).  We also add the total
flux in positive pixels (\code{fGood}) and total absolute value of the
flux in negative pixels (\code{fBad}).  Using these values, We report
the following quantities:
\begin{itemize}
\item \code{nGood}
\item \code{fRatio} = \code{fGood} / (\code{fGood} + \code{fBad})
\item \code{nRatioBad} = \code{nGood} / (\code{nGood} + \code{nBad})
\item \code{nRatioMask} = \code{nGood} / (\code{nGood} + \code{nMask})
\item \code{nRatioAll} = \code{nGood} / (\code{nGood} + \code{nMask} + \code{nBad})
\end{itemize}

We also attempt to place the difference image detections in the
context of the input images, both the positive (subtrahend) and
negative (minuend) images.  We identify the closest source in both the
positive and negative images to the detection in the difference image,
out to a maximum of \code{INPUT.MATCH.RADIUS} (50 pixels for PV3), but only
if the source in those images has a signal-to-noise greater than
\code{INPUT.MATCH.MIN.SN} (10 for PV3).  If there is a close neighbor in the
positive image, and the difference in the magnitudes of the source in
that image and the source in the difference image is less than 5
$\sigma$, then the bit \code{PM_SOURCE_MODE2_DIFF_SELF_MATCH =
  0x00000800} is raised in \code{mask2} as these two detections are
likely the same flux (\ie, detection of an isolated source).  

If the difference image detection is matched to a nearby source in the
positive image, then the signal-to-noise of the neighbor is saved as
\code{DIFF_SN_P} and the distance in pixels between the difference
detection and positive detection is saved as \code{DIFF_R_P}.
Similarly, for a neighbor in the negative image, these values are
saved as \code{DIFF_SN_M} and \code{DIFF_R_M}.  Additional
\code{mask2} bits are also raised: if the difference detection is only
associated with one of the two input images, then the bit
\code{PM_SOURCE_MODE2_DIFF_WITH_SINGLE = 0x00000001} is raised, while
a difference detection which has a match in both input images has
\code{PM_SOURCE_MODE2_DIFF_WITH_DOUBLE = 0x00000002} raised. 

Comets appear in the difference images as a non-PSF sources.  Their
2-D structure includes both the flux from the coma (with a typical
power-law profile) and flux from the tail (with a more complex flux
distribution).  We use the Kron magnitudes to identify possibly
extended objects which may be cometary in nature.

% \note{need some info from MOPS folks on what is used}

For a difference image, both positive and negative sources will be
present.  The basic peak detection algorithm will only trigger for the
positive sources.  In the \ippprog{ppSub} program, both the $A - B$
and the $B - A$ images are sent to the \ippprog{psphot} routine for
source detection and characterization.

Note that the variance image for a difference image must be generated
from the two positive images used to construct the difference.  It is
possible to run \ippprog{psphot} as an external program on a
difference image generated previously.  In this case, the variance
image and the PSF model must be supplied as well as the difference
image.

\section{Conclusions}

The Pan-STARRS Image Processing Pipeline has used the \ippprog{psphot}
software to detect and characterize astronomical sources in images
from both the PS\,1 and PS\,2 telescopes since 2008.  This software
system has produced highly-reliable stellar photometry and astrometry
measurements as demonstrated by the high-quality data products
released as part of the Pan-STARRS Data Releases 1 and 2.  This
configurable software system has been used for stellar photometry in
direct detection model, difference image mode, and forced photometry
mode, as well as galaxy photometry and morphology analysis.  To date
(2019 May), over 900 billion detections of sources in more than 2
million PS\,1 exposures have been characterized (some representing
repeated measurements of the same exposures).  

There is always room for improvement, however.  A number of
possible improvements to \ippprog{psphot} have been identified which
could result in more reliable measurements for either stars or
galaxies.  Here we discuss improvements beyond simply tuning
parameters for a specific dataset.

In general, the improvements we identify share the characteristic of
making use of external information in the analysis.  As described
above, essentially all operations of \ippprog{psphot}, except in the
context of forced photometry, approach each image with no prior
knowledge.  This was necessary in the early stages of the Pan-STARRS
project when we had not yet observed the sky with our instrument and
comparable observations were only available in the SDSS Galactic cap
regions.  However, the sky is now much better known, not only from
PS1, but also for example due to Gaia.

Several improvements to the \ippprog{psphot} analysis could be made by
including as much information from external catalogs about the
positions and characteristics of sources in the images as possible.
For example, known stars (e.g., based on proper motions from Gaia or
colors and morphology from PS1) could be used for PSF sources.  In
areas of high density, especially in known globular or even open
clusters, existing high-resolution imagery could be used to provide a
constraint on location of stars.  External information could also be
used to control the scale on which the background is modelled: a finer
sampling is helpful in regions of known nebulosity and large galaxies
such as M31.  Finally, the galactic latitude or the externally-defined
stellar density could be used to control the choice of fitting double
stars or galaxy models.  This would be a step beyond the current
capability of choosing to fit galaxy models as a function of galactic
latitude.

% PS2 reference:
% https://ui.adsabs.harvard.edu/abs/2012SPIE.8444E..0HM/abstract

\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}

\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?

% programs mentioned in this text:
% psphot
% psphotStack
% psphotFullForce
% ppImage
% ppSub

