Next: Joint Processing of Radio Data Produced by the SSRT Together With Data of Other Spectral Ranges
Up: Data Analysis and Processing Techniques
Previous: The Pixon Method of Image Reconstruction
Table of Contents - Subject Index - Author Index - Search - PS reprint - PDF reprint

Mighell, K. J. 1999, in ASP Conf. Ser., Vol. 172, Astronomical Data Analysis Software and Systems VIII, eds. D. M. Mehringer, R. L. Plante, & D. A. Roberts (San Francisco: ASP), 317

Algorithms for CCD Stellar Photometry

Kenneth J. Mighell
Kitt Peak National Observatory, National Optical Astronomy Observatories, P.O. Box 26732, Tucson, AZ  85726-67332

Abstract:

Basic algorithms for the production of stellar photometry from CCD observations are described and reviewed. Simple algorithms for detecting stars, finding their positions, and estimating the nearby sky background will be presented along with references to the extensive literature devoted to these topics. The basic techniques of aperture and PSF-fitting CCD stellar photometry are discussed in detail. Recent advances in PSF-fitting CCD stellar photometry are described and used to show how accurate V-band stellar photometry with the Next Generation Space Telescope could be accomplished with a 8-m primary mirror that is diffraction limited at near-infrared wavelengths.

1. Introduction to CCD Stellar Photometry

The problem of doing accurate stellar photometry with calibrated charge-coupled device (CCD) data amounts to that of condensing the intensity values of hundreds, thousands, or even millions of picture elements (pixels) into a list containing the magnitude and position of each star on the CCD image.

If the stellar field is not crowded, the astronomer can measure the magnitude of each star by doing the digital equivalent of aperture photometry. Stellar photometry becomes much more difficult if one wishes to study crowded stellar fields. The basic assumptions underlying the simple techniques of CCD aperture photometry are frequently not valid in crowded stellar fields and more sophisticated methods like Point-Spread-Function (PSF) model fitting must be used in order to achieve accurate photometry.

2. CCD Aperture Stellar Photometry

The process of determining the apparent magnitude of a star can be surprisingly complex even with a technique as simple as CCD aperture photometry. We start by assuming that our CCD observations have already been flat-fielded and calibrated. The basic process of aperture stellar photometry is, in principle, quite simple:

  1. Find the star.
  2. Center an aperture of NA pixels on the star.
  3. Add up the electrons within the aperture: SA (units: electrons e-).
  4. Determine the nearby background flux: B  (units: e-pixel-1).
  5. Determine the instrumental magnitude: $m^\prime = -2.5\log\left[({S_A - N_AB})/1\,\mbox{e$^-$}\right]$.
  6. Determine the aperture correction: $\Delta m^\prime$  (units: mag).
  7. Compute the aperture-corrected instrumental magnitude: $m = m^\prime + \Delta m^\prime$.

The important process of transforming instrumental magnitudes to a standard system is extensively described in the literature and will not covered in this contribution due to space limitations.

2.1 Finding Stars

There are many techniques available for the detection of astronomical objects in CCD observations (see, e.g., Fischer & Kochanski 1994, Secker 1995, and references therein). The following is a short introduction to some of the techniques that can be used to detect stars in CCD observations.

The signal of a noisy digital image can frequently be enhanced by suppressing high spatial frequency noise in the image. For modern CCD stellar observations, this generally means suppressing photon and CCD readout noise. This can frequently be accomplished by using small digital low-pass filters like

\begin{displaymath}
\mbox{\bf {LP}}_{3\times3}
=
\left(
\begin{array}{ccc}
1/16 ...
.../16\\
1/8 & 1/4 & 1/8\\
1/16 & 1/8 & 1/16
\end{array}\right)
\end{displaymath}

which is commonly used in digital image processing for this purpose. Small digital low-pass filters are useful in finding stars in digital stellar observations that have had the background sky removed (Irwin 1985). The background sky of an image can frequently be approximated by using small median filters. Combining low-pass filters with median filters can be very useful; the LPD (low-pass difference) filter,

\begin{displaymath}
\mbox{\bf {LPD}}_{3\times3}({\cal F})
=
\mbox{\bf {LP}}_{3\t...
...{5\times5}
\left[
\mbox{\bf {LP}}_{3\times3}({\cal F})
\right]
\end{displaymath}

works well with Hubble Space Telescope WF/PC and WFPC2 images, ${\cal F}$, for the purpose of detecting stars and other point sources (Appendix A of Mighell & Rich 1995). The LPD filter is a high-pass frequency filter which is effectively a digital analog of Malin's unsharp photographic masking technique (Malin 1977, 1981, and references therein).

The following listing, for example, is a FORTRAN implementation of a simple peak detector algorithm which identifies any pixel that is greater than any of its eight neighbors:


C
C A Simple Peak Detector Algorithm 
C
C Copyleft (L) 1998 Kenneth J. Mighell (Kitt Peak National Observatory)
C
      SUBROUTINE PEAKER(IMAGE,NX,NY,PEAKMIN,PEAKMAX)
      INTEGER NX, NY, X, Y, XX, YY
      REAL    IMAGE(NX,NY), PEAKMIN, PEAKMAX, PIXEL, NEIGHBOR
      LOGICAL BINGO
      DO Y = 2,(NY-1)
      DO X = 2,(NX-1)
      PIXEL = IMAGE(X,Y)
      IF ((PIXEL.GE.PEAKMIN).AND.(PIXEL.LT.PEAKMAX)) THEN
        BINGO = .TRUE.
        DO YY = (Y-1),(Y+1)
        DO XX = (X-1),(X+1)
          IF (BINGO) THEN
            NEIGHBOR = IMAGE(XX,YY)
            IF (NEIGHBOR.GT.PIXEL) THEN
              BINGO = .FALSE.
            ELSE IF (NEIGHBOR.EQ.PIXEL) THEN
              IF ((XX.NE.X).OR.(YY.NE.Y)) THEN
                IF (((XX.LE.X).AND.(YY.LE.Y))
     &             .OR.((XX.GT.X).AND.(YY.LT.Y))) BINGO = .FALSE.
              ENDIF
            ENDIF
          ENDIF
        ENDDO
        ENDDO
        IF (BINGO) WRITE (*,*) 
     & 'Found a peak at position (',X,',',Y,') with a value of',PIXEL
      ENDIF
      ENDDO
      ENDDO
      RETURN 
      END
C23456
This algorithm works particularly well with CCD stellar observations that have had the background sky removed (e.g., LPD-filtered images).

2.2 Determining the Centers of Stars

Many CCD aperture photometry programs require the user to give the position of the star on the CCD frame. Most of these programs offer a choice of centering algorithms to determine the center of a star when provided with only a rough estimate. The review article of Stone (1989) compares the performance of five different digital centering algorithms under a wide range of atmospheric seeing and background-level conditions. It may be useful to create your own centroid algorithm based upon the special requirements of your particular analysis problem. The following listing, for example, is a FORTRAN implementation of a centroid algorithm which produces robust estimates without requiring an estimate of the the nearby ``sky'' background:


C
C A Simple Centroid Algorithm
C
C A modification of the Modified Moment Method (Stone 1989,AJ,97,1227)
C that works well with small apertures.
C
C Copyleft (L) 1998 Kenneth J. Mighell (Kitt Peak National Observatory)
C
      SUBROUTINE CENTROID(X,Z,NN,N,X_IN,X_OUT)
      INTEGER          N, NN ! <-- assumes that 1 <= N <= NN
      REAL             X(NN), Z(NN), X_IN, X_OUT, BIG
      DOUBLE PRECISION DELTA, SUM1, SUM2, DIFF, XX
      DOUBLE PRECISION INTENSITY, POSITION, MINIMUM
      INTEGER          I, J , NITERATIONS
      PARAMETER        (BIG=1E30,NITERATIONS=10)
      MINIMUM = BIG
      DO I = 1,N
        INTENSITY = Z(I) ! <-- Z(I) is the intensity at the position X(I)
        MINIMUM = MIN(INTENSITY,MINIMUM)
      ENDDO
      XX = X_IN
      DELTA = 0D0
      DO J = 1,NITERATIONS
        XX = XX + DELTA
        SUM1 = 0D0
        SUM2 = 0D0
        DO I = 1,N
          POSITION = X(I)
          INTENSITY = Z(I)
          DIFF = MAX((INTENSITY-MINIMUM),0D0)
          SUM1 = SUM1 + (POSITION-XX)*DIFF
          SUM2 = SUM2 + DIFF
        ENDDO
        DELTA = SUM1/SUM2
      ENDDO
      X_OUT = XX
      RETURN
      END
C23456

This algorithm is a simplified version of the one used by my CCDCAP aperture stellar photometry package. If precise absolute or relatives positions from CCD observations are required, then one should investigate the extensive literature dedicated to CCD astrometry.

2.3 Determining the Nearby Background

The background (``sky'') flux associated with a star is generally determined by analyzing the distribution of intensities of nearby background pixels. In the case of circular aperture photometry, the background flux is typically determined by analyzing the pixels in an annulus beyond the stellar aperture. The inner radius of the sky annulus is typically several FWHM1distant from the center of the aperture in order to avoid the inclusion of contaminating light from the star itself. The width of the annulus is typically large enough so that the annulus contains between $\sim$50 and a few hundred pixels.

Many aperture photometry programs allow the user to set the background flux to be the modal value of the background intensity distribution. The mode of the background distribution is frequently estimated by using the following useful approximation,

\begin{displaymath}
\mbox{mode}
\approx
3\!\times\!\mbox{median}
-
2\!\times\!\mbox{mean}
~~~\mbox{(for median$\,\leq\,$mean)}
,
\end{displaymath}

(Wells 1979; Kendall & Stuart 1958; Haldane 1942; Pearson 1895). This approximation, unfortunately, is known to produce background estimates that are biased towards higher values (e.g., Newberry 1992). Better methods for the estimation of the background are available. Modal estimates of the background distribution after the iterative rejection of outlier pixel intensities beyond 2.5-3.0 standard deviations of the mean generally produce reasonable results (Da Costa 1992).

Many algorithms are available and most aperture photometry programs offer the user a choice of several methods to determine the background flux. For example, the popular APPHOT aperture photometry IRAF package offers the user the choice of 11 different ways the background can be estimated (Davis 1987). As always when using analysis software, the astronomer is strongly advised to carefully read the user documentation in order to understand which methodology will produce the best results for a given CCD observation or application.

2.4 Determining the Stellar Signal within an Aperture

The fundamental task of a CCD aperture photometry program is to accurately measure all the electrons, SA, that fall within an aperture placed on a CCD image. Unless the background flux, B, is exactly zero electrons per pixel, one can never directly measure the number of electrons, $S_\star$, from a star within the aperture. One measures, instead, the quantity $S_A = S_\star + N_{\!A} B$, where $N_{\!A}$ is the aperture of the aperture (in pixels). The instrumental magnitude of an aperture measurement of a CCD observation of a star can thus be defined as

\begin{displaymath}
m
\equiv
-2.5 \log\left[ \frac{S_\star}{1\,\mbox{e$^-$}} \ri...
....5 \log\left[ \frac{S_A - N_{\!A} B}{1\,\mbox{e$^-$}} \right]
\end{displaymath}

where the final form on the right is in terms of observable quantities. By this definition, a stellar intensity of one electron would have an instrumental magnitude of zero. The estimated error of the instrumental magnitude is approximately

\begin{displaymath}
\Delta m
\approx
1.0857\left[\frac{\Delta S_\star}{S_\star}\right]
\nonumber
\end{displaymath}

where $\Delta S_\star$ is the measurement error (in electrons) of the estimated number of electrons from the star within the aperture.

It is clearly important to determine the observable quantities SA, B, and $N_{\!A}$ as accurately as possible. Let us consider the aperture area, $N_{\!A}$, first. If aj is the area of the jth aperture pixel, then the total area of the aperture, $N_{\!A}$, is simply

\begin{displaymath}
N_{\!A} = \sum_{j = 1}^{A} a_j\ ,
\end{displaymath}

where A indicates that the summation includes all pixels or partial pixels within the aperture. The area of the jth aperture pixel will be exactly one pixel only when the pixel is completely within the aperture, otherwise only a fraction of the pixel lies within the aperture and the value of aj is between zero and one ( $0\!<\!a_j\!<\!1$ pixel). The area of a circular aperture is $A = \pi r^2$ pixels where the radius of the aperture is r pixels. Although it is not difficult to exactly determine the partial pixel areas with square pixels and circular apertures (e.g., Fig. 1), several standard CCD aperture photometry programs only approximate the partial pixel areas. For example, the popular PHOT task of the APPHOT package approximates the circular aperture by an irregular polygon (Davis 1987). While this approximation is generally fine with large apertures, it can sometimes produce large systematic measurement errors for small aperture radii.

Figure 1: A circular aperture on an array of square CCD pixels. The black pixels are entirely within the aperture and the gray pixels are only partially within the aperture.
\begin{figure}
\epsscale{0.4}
\plotone{mighellkj1.eps}
\end{figure}

Let us now consider the sum of all electrons within the aperture (SA). If zj is the intensity (in electrons) of the jth aperture pixel, then one simple approximation of SA is

\begin{displaymath}
S_A = \sum_{j = 1}^{A} a_jz_j\ .
\end{displaymath} (1)

This approximation linearly weights the pixel intensity with the area of the pixel within the aperture. The PHOT task uses a linear pixel weighting algorithm that is very similar to this approximation. The use of Equation (1) implicitly assumes that the Point Spread Function is nearly flat at the edge of the aperture. The PSF is generally nearly flat only at large distances from the center of a star where the encircled-energy function is nearly equal to one (i.e., 100%). The use of Equation (1) is thus generally appropriate for large aperture radii. Wherever the PSF is a rapidly changing function of radius (e.g., at small radii of seeing-optimized CCD observations) the use of Equation (1) is likely to produce large systematic measurement errors.

One way of reducing this systematic measurement error is to split each pixel into subpixels by using a bilinear pixel interpolation algorithm. One such algorithm is to use the two-dimensional analog of the sinc function: $J_1(\pi r)/2r$. This function, like many others, has the unfortunate effect of degrading the original image by spreading photons (electrons) beyond the original pixel. I created created the QUADPX bilinear pixel interpolation algorithm which splits a pixel into 4 subpixels whose sum is always equal to that of the original pixel (Appendix B of Mighell & Rich 1995). Numerical experiments have shown that using the QUADPX algorithm can reduce the systematic measurement error of critically-sampled Gaussians by a factor of $\sim$6. For example, the photometric error of an aperture radius of 2.0 pixels went from 0.068 mag using Equation (1) to 0.011 mag with the CCDCAP package which implements the QUADPX algorithm.

2.5 Optimal Aperture Size and Aperture Corrections

The best (smallest) stellar photometric errors (i.e., the largest signal-to-noise ratios) are generally obtained with relatively small apertures (see, e.g., Fig. 6 of Howell 1989). Analysis of theoretical CCD signal-to-noise-ratio equations (see, e.g., Newberry 1991, Howell 1992, Merline & Howell 1995, Howell et al. 1996, and references therein) shows that large apertures can have large photometric errors when the the total number of stellar photons in the aperture becomes comparable with the total number of background photons in the aperture. Furthermore, a measurement error for the background flux as small as just 1 electron per pixel can by itself produce large photometric uncertainties at large aperture radii. Small apertures, however, can be too small when they allow such a small fraction of the star light to be found within the aperture that the photometric error becomes dominated by small-number (a.k.a. counting or Poisson) statistics because little or no signal has been measured.

A Gaussian is a good model for the Point Spread Function of a ground-based CCD observation since the central core of a ground-based stellar profile is approximately Gaussian (King 1971). One can easily show that the optimum signal-to-noise ratio for a Gaussian PSF is obtained for a circular aperture radius of $\sim$1.6$\sigma$ (i.e.,  $r \approx 0.68\,$FWHM) which contains about 72% of the encircled-energy. Pritchet & Kline (1981) note that the signal-to-noise ratio is fairly insensitive to radius near the ``optimal'' radius value $\sim$1.6$\,\sigma$ for a Gaussian PSF; deviations from the optimal radius by as much as $\pm$50% generally make little difference. Since centering errors will be more critical for smaller apertures than for larger apertures, it is only prudent to err on the larger side by using apertures with radii which are larger than $\sim$0.68$\,$FWHM.

An aperture radius of $r\approx\,$FWHM makes an excellent practical compromise between concerns about systematic centering errors and diminishing signal-to-noise ratios typically obtained with larger aperture radii. By analyzing the CCD signal-to-noise ratio equations, one can show that brighter stars will have larger optimal aperture sizes than do fainter stars. If one must use only one aperture size, then it is clearly advantageous to chose a global aperture size which produces the smallest photometric errors for the faintest stars (i.e., use $r\approx\,$FWHM).

Small apertures frequently do not contain all the flux from a star. The amount of the missing star light can found by determining the appropriate aperture correction by measuring nearby bright isolated stars. Howell (1989) and Stetson (1990), among others, describe the process how aperture corrections can be accurately determined using the aperture growth-curve method.

3. CCD PSF-Fitting Stellar Photometry

Consider a ground-based CCD observation of two stars whose stellar images overlap. Assuming we already know the Point Spread Function of the observation, a simple model of the observation will have seven parameters: peak intensities (I1,I2), positions (X1,Y1,X2,Y2), and the background sky level B which is assumed to be the same for both component images. One finds that the parameters are not independent for overlapping stars with the presence of photon and readout noise. The conservation of photon flux will require that if I1 increases then I2 must decrease and vice versa for a given value of B. The most accurate photometry possible is obtained when these dependent parameters are fitted simultaneously. Any reasonable model of two overlapping stellar images will be a non-linear function when the positions and peak intensities are to be determined simultaneously. The technique of non-linear least-squares fitting was developed to provide for the simultaneous determination of dependent or independent parameters of non-linear model functions.

Assume that we have a calibrated CCD observation with N pixels and that zi is the intensity in electrons (e-) of the ith pixel at (xi,yi) with an error of $\sigma_i$. Let $\mbox{\em {model}}( x, y,;$a1,$\ldots,$aM) be a model of the intensity values that has two coordinates (x,y) and M parameters. For notational convenience, let the vector ri represent the coordinates (xi,yi) and the vector a represent all the parameters [i.e.,  $\mbox{\boldmath$a$}\equiv (a_1,\ldots,a_M)$]. Thus the model of intensities will normally be written as $\mbox{\em {model}}(\mbox{\boldmath$r$}_i;\mbox{\boldmath$a$})$.

The measure of the goodness of fit between the data and the model, called chi-square, is defined as

\begin{displaymath}
\mbox{$\chi^2(\mbox{\boldmath$a$})$}\equiv \sum_{i=1}^{N}
\f...
...em {model}}(\mbox{\boldmath$r$}_i;\mbox{\boldmath$a$})]^2\ \ .
\end{displaymath}

The theory of least squares states that the optimum value of the parameter vector a is obtained when $\chi^2(\mbox{\boldmath$a$})$ is minimized with respect to each parameter simultaneously. If the function $\chi^2(\mbox{\boldmath$a$})$ is thought of as a surface in the M-dimensional parameter space and if $\mbox{\boldmath$a$}_{\mbox{true}}$ is defined as the optimal parameter vector, then the absolute minimum of that surface is $\chi^2(\mbox{\boldmath$a$}_{\mbox{true}})$. The surface of the function $\chi^2(\mbox{\boldmath$a$})$ might be very complicated and the fitting algorithm must be able to converge to believable solutions (answers) even in ill-defined situations.

For some small correction parameter vector $\mbox{\boldmath$\delta$}$ we can approximate $\chi^2(\mbox{\boldmath$a$}+\mbox{\boldmath$\delta$})$ by its Taylor series expansion:

\begin{displaymath}
\chi^2(\mbox{\boldmath$a$}+\mbox{\boldmath$\delta$})
=
\sum_...
...delta$}\cdot \mbox{\boldmath$H$}\cdot \mbox{\boldmath$\delta$}
\end{displaymath}

where

\begin{displaymath}[H]_{jk} \equiv \frac{\partial^2\,\mbox{$\chi^2(\mbox{\boldmath$a$})$}}{\partial a_j \partial a_k}
\end{displaymath}

the jkth element of the M x M Hessian matrix H of $\chi^2(\mbox{\boldmath$a$})$ [see, e.g., Arfken 1970; Press et al. 1986]. If $\chi^2(\mbox{\boldmath$a$}+\mbox{\boldmath$\delta$})$ is a local minimum of $\chi^2$, then it can be shown that

\begin{displaymath}
\mbox{\boldmath$H$}\cdot \mbox{\boldmath$\delta$}= -\nabla\mbox{$\chi^2(\mbox{\boldmath$a$})$}\ .
\end{displaymath}

By solving this equation for the correction vector $\delta$ it is then possible to determine a better parameter vector $
\mbox{\boldmath$a$}_{\mbox{new}} = \mbox{\boldmath$a$}_{\mbox{old}} + \mbox{\boldmath$\delta$}\ .
$ When the parameter vector ( a) is redefined to be the better parameter ( $\mbox{\boldmath$a$}_{\mbox{new}}$), the Hessian matrix and the gradient of $\chi^2(\mbox{\boldmath$a$})$ can then be recalculated to determine a new correction vector ( $\delta$). This process repeats until $\delta$ is sufficiently small. The final parameter vector is called the optimal parameter vector, $\mbox{\boldmath$a$}_0$, and should be very close to the optimal parameter vector ( $\mbox{\boldmath$a$}_{\mbox{true}}$) if the fit is good.

This method will find the absolute minimum of $\mbox{$\chi^2(\mbox{\boldmath$a$})$}$ if the original guess $\mbox{\boldmath$a$}$ is near $\mbox{\boldmath$a$}_{\mbox{true}}$. Unfortunately, the original guess of the parameter vector $\mbox{\boldmath$a$}$ may not always be very good. For production stellar photometry software it is important that the search for the absolute minimum of $\chi^2$ be both robust and efficient.

3.1 Analytical Models

One can frequently create a realistic intensity model of a ground-based CCD observation of a total of K stars on a non-flat background with a combination of Moffat (1969) functions on a tilted plane:

\begin{eqnarray*}
\mbox{\em {model}}(\mbox{\boldmath$r$};\mbox{\boldmath$a$})
&\...
...&\equiv
&B_0 + B_Xx + B_Yy +
\sum_{k=1}^{K} I_k \mbox{$\Psi$}_k
\end{eqnarray*}


where the B0, BX, BY terms describe the tilted-plane model of the background and Ik is the intensity of the kth star which has a Point Spread Function of
\begin{displaymath}
\mbox{$\Psi$}_k
\equiv
\int_{x-0.5}^{x+0.5}
\int_{y-0.5}^{y+...
...\rho_k^2}
\right\}
\right]^{-\tau_k}
\mbox{d}y
\
\mbox{d}x
~.
\end{displaymath} (2)

The partial derivatives of $\mbox{\em {model}}(r;a)$ must now be determined:

\begin{displaymath}
\frac{\partial \mbox{\em {model}}}{\partial B_0}
= 1,~~~
\fr...
...{\partial \mbox{\em {model}}}{\partial I_k}
=
\mbox{$\Psi$}_k,
\end{displaymath}


\begin{displaymath}
\frac{\partial \mbox{\em {model}}}{\partial X_k}
=
I_k \frac...
...tau_k}
=
I_k \frac{\partial \mbox{$\Psi$}_k}{\partial \tau_k}.
\end{displaymath}

If the PSF, $\mbox{$\Psi$}_k$, does not have a general analytical solution, it can usually be approximated numerically. For example, if we subdivide a pixel into $\eta^2$ subpixels, then Equation (2) can be approximated as $\mbox{$\Psi$}_k \approx \mbox{$\gamma_k^{-\tau_k}$}$ where

\begin{displaymath}
\mbox{$\gamma_k^{-\tau_k}$}
\equiv
\!\!\frac{1}{\eta^2}
\sum...
...i+\eta^{-1}j]-\!Y_k)^2
}{\rho_k^2}
\right\}
\right]^{-\tau_k}
\end{displaymath}

and $\xi\equiv(\eta+1)/(2\eta)$. The new partial derivatives of $\mbox{\em {model}}(r;a)$ can now be determined as follows

\begin{displaymath}
\frac{\partial \mbox{\em {model}}}{\partial I_k}
\approx
\mb...
...c{2 I_k \tau_k}{\rho_k^2} \frac{(y-Y_k)}{\gamma_k^{\tau_k+1}},
\end{displaymath}


\begin{displaymath}
\frac{\partial \mbox{\em {model}}}{\partial \rho_k}
\approx
...
...au_k}
\approx
-I_k \mbox{$\gamma_k^{-\tau_k}$}\ln(\gamma_k)\,.
\end{displaymath}

Simulations and practical experience has demonstrated that the $\gamma_k^{-\tau_k}$ approximation with the Levenberg-Marquardt method of non-linear least squares (Levenberg 1944, Marquardt 1963) produces accurate and precise CCD stellar photometry (Mighell 1989, 1990).

3.2 Digital Models

One can create a realistic intensity model of a Hubble Space Telescope ( HST) CCD observation of K stars on a non-flat background with a combination of digital Point Spread Functions on a tilted plane:

\begin{eqnarray*}
\mbox{\em {model}}(\mbox{\boldmath$r$};\mbox{\boldmath$a$})
&\...
...v
&B_0 + B_Xx + B_Yy +
\sum_{k=1}^{K}
I_k \Psi_k (x-X_k,y-Y_k)
\end{eqnarray*}


where the B0, BX, BY terms describe the tilted-plane model of the background and Ik is the intensity of the kth star which has a digital PSF, $\Psi_k$, that is represented with matrix of numbers with a sum of one. If the kth star is isolated and on a flat background, then the digital PSF, $\mbox{$\Psi$}_k$, can be easily derived from the actual image. Alternatively, $\mbox{$\Psi$}_k$ could be a synthetic PSF computed by PSF modeling software [e.g., the TINY TIM package of Krist (1993, Krist & Hook 1997)]. With a digital PSF in hand, one can then easily determine the partial derivatives with respect to x and y using standard numerical differentiation techniques (see Fig. 2). The art of CCD photometry with digital Point Spread Functions lies in the implementation details. While the mathematics is the same as for the case of analytical PSFs, the software engineering issues are significantly more challenging.

Figure 2: Typical Point Spread Functions ($\Psi $) of the Hubble Space Telescope WFPC2 (left) and WF/PC (right) instruments (M star, F555W filter) and their respective partial derivatives with respect to the directions x $(\partial \Psi / \partial x)$ and y $(\partial \Psi / \partial y)$. These PSFs were synthesized using TINY TIM VERSION 4.0B (Krist 1994).
\begin{minipage}[h]{\textwidth} %5.25truein\}
\vspace*{-105truemm}
\hspace*{-5...
...textwidth}
{\bf {\hskip 26truemm}{WFPC2}{\hskip 49truemm}{WF/PC}}
\end{minipage}

3.3 Recent Advances

Traditional PSF-fitting CCD photometric reduction packages like DAOPHOT (Stetson 1987) use analytical functions to represent the Point Spread Function. All the major partial derivative computations are computed on the analytical model of the PSF. Any deviations of the real-world PSF from the analytical PSF are then generally stored in a residual matrix which is only used to determine the $\chi^2$ goodness-of-fit.

I have recently demonstrated the feasibility of doing accurate CCD stellar photometry with digital Point Spread Functions. I have developed a new digital PSF-fitting algorithm which does not require a residual matrix because all partial derivative computations are done on the digital PSF itself using standard numerical differentiation techniques. This algorithm has already passed the proof-of-principle stage with the successful reduction of simulated Next Generation Space Telescope (NGST) CCD stellar observations (see Fig. 3.).

I investigated the performance of CCD stellar photometry with a 1.5-micron diffraction-limited 8-m Next Generation Space Telescope. These simulations used artificial Point Spread Functions for three different 8-m NGST design concepts which were kindly provided by John Krist. Assuming that the 8-m NGST primary mirror has 1/13 wave RMS errors at 1.5 micron, I determined that 90% of the light from a star falls within an aperture radius of 0.1 arcsec -- the size of one WF pixel of the Hubble Space Telescope WFPC2 instrument. The three NGST design concepts have nearly identical V-band encircled-energy functions; the degradation caused by the differences between the three NGST design concepts is quite negligible when state-of-the-art digital-PSF photometric reduction software is used to analyze uncrowded stellar fields.

Figure 3: Measured signal-to-noise ratios of simulated NGST CCD stellar observations using the Levenberg-Marquardt method of non-linear least squares fitting with digital Point Spread Functions. Each mark represents the average signal-to-noise ratio of 100 stars that were simulated using using synthetic V-band PSFs on a flat background of 300 electrons per pixel. The filled circles show the result using a critically sampled V-band PSF for a perfect 0.5-$\mu$m diffraction-limited 8-m NGST mirror. The V-band results for three different 8-m NGST design concepts of a 1.5-$\mu$m diffraction-limited mirror.
\begin{figure}
\epsscale{0.4}
\plotone{mighellkj3.eps}
\end{figure}

All the V-band PSFs used in these simulations were sampled at 0.0064 arcsec pixel-1 which is the critical V-band sampling rate of a perfect 0.5-micron diffraction-limited 8-m NGST. Better photometric performance for all three NGST design concepts could be obtained by using larger CCD pixels: a pixel size of $\sim$0.013 ($\sim$ 2 x 0.0064) arcsec pixel-1 should provide an optimal pixel sampling when the 8-m primary mirror is diffraction-limited at 1.5 micron.

References

Arfken, G. 1970, Mathematical Methods for Physicists (2nd ed.), (New York: Academic Press)

Da Costa, G. S. 1992, in ASP Conf. Ser., Vol. 23, Astronomical CCD Observing and Reduction Techniques, ed. S. B. Howell, (San Francisco: ASP), 90

Davis, L. 1987, ``Specifications for the Aperture Photometry Package'', National Optical Astronomy Observatories, ftp://iraf.noao.edu/iraf/docs/apspec.ps.Z

Fischer, P. & Kochanski, G. P. 1994, AJ, 107, 802

Haldane, J. B. S. 1943, Biometrika, 32, 294

Howell, S. B. 1989, PASP, 101, 616

\ibid , 1992, in ASP Conf. Ser., Vol. 23, Astronomical CCD Observing and Reduction Techniques, ed. S. B. Howell, (San Francisco: ASP), 105

Howell, S. B., Koehn, B., Bowell, E., & Hoffman, M. 1996, AJ, 112, 1302

Irwin, M. J. 1985, MNRAS, 214, 575

Kendall, M. G. & Stuart A. 1958, The Advanced Theory of Statistics, Vol. I, (London: Charles Griffin & Co.) pp. 39 and 179

King, I. R. 1971, PASP, 83, 199

Krist, J. 1993, in ASP Conf. Ser., Vol. 52, Astronomical Data Analysis Software and Systems II, ed. R. J. Hanisch, R. J. V. Brissenden, & J. Barnes (San Francisco: ASP), 530

\ibid , 1994, Tiny Tim User's Manual (Ver. 4.0)

Krist, J. & Hook, R. 1997, Tiny Tim User's Guide (Ver. 4.4)

Levenberg, K. 1944, Quart. of Appl. Math., 2, 164

Malin, D. 1977, AAS Photo-Bull., No. 16, 10

\ibid , 1981, J. Phot. Sci., 29(5), 199

Marquardt, D. 1963, J. SIAM, 11, 431

Merline, W. J. & Howell, S. B. 1995, Exp. Astron., 6, 163

Mighell, K. J. 1989, MNRAS, 238, 807

\ibid , 1990, A&AS, 82, 1

Mighell, K. J. & Rich, R. M. 1995, AJ, 110, 1649

Moffat, A. F. J. 1969, A&A, 3, 455

Newberry, M. V. 1991, PASP, 103, 122

\ibid , 1992, in ASP Conf. Ser., Vol. 25, Astronomical Data Analysis Software and Systems I, ed. D. M.  Worrall, C. Biemesderfer, & J. Barnes (San Francisco: ASP), 307

Pearson, K. 1895, Phil. Trans., 186, 343

Press, W. H., Flannery, B. P., Teukolsky, S. A., & Vetterling, W. T. 1986, Numerical Recipes, (Cambridge: Cambridge Univ. Press)

Pritchet, C. & Kline, M. I. 1981, AJ, 86, 1859

Secker, J. 1995, PASP, 107, 496

Stetson, P. B. 1987, PASP, 99, 191

\ibid , 1990, PASP, 102, 932

Stone, R. C. 1989, AJ, 97, 1227

Wells, D.C. 1979, in SPIE Proc., Vol. 172, Instrumentation in Astronomy III, ed. D. L. Crawford, (Bellingham: SPIE), 418



Footnotes

... FWHM1
The Full-Width-at-Half-Maximum of a critically-sampled Gaussian distribution is $\sim$2.35 pixels.

© Copyright 1999 Astronomical Society of the Pacific, 390 Ashton Avenue, San Francisco, California 94112, USA
Next: Joint Processing of Radio Data Produced by the SSRT Together With Data of Other Spectral Ranges
Up: Data Analysis and Processing Techniques
Previous: The Pixon Method of Image Reconstruction
Table of Contents - Subject Index - Author Index - Search - PS reprint - PDF reprint

adass@ncsa.uiuc.edu