idlastro / DAOPHOT-Type Photometry Procedures: GETPSF

[Source code]

NAME
GETPSF
PURPOSE
To generate a point-spread function (PSF) from observed stars. 
EXPLANATION
The PSF is represented as a 2-dimensional Gaussian
(integrated over each pixel) and a lookup table of residuals.
The lookup table and Gaussian parameters are output in a FITS
image file.   The PSF FITS file created by GETPSF can be
read with the procedure RDPSF.      Adapted from the 1986 STSDAS 
version of DAOPHOT
CALLING SEQUENCE
GETPSF, image, xc, yc, apmag, sky, [ronois, phpadu, gauss, psf, 
                idpsf, psfrad, fitrad, psfname, /DEBUG ]
INPUTS
IMAGE  - input image array
XC     - input vector of x coordinates (from FIND), these should be
        IDL (first pixel is (0,0)) convention.
YC     - input vector of y coordinates (from FIND)
APMAG  - vector of magnitudes (from APER), used for initial estimate
        of gaussian intensity.  If APMAG is multidimensional, (more
        than 1 aperture was used in APER) then the first aperture
        is used.
SKY    - vector of sky values (from APER)                
OPTIONAL INPUTS
The user will be prompted for the following parameters if not supplied.
RONOIS - readout noise per pixel, (in electrons, or equivalent photons)
PHPADU - photons per analog digital unit, used to scale the data
        numbers in IMAGE into photon units
IDPSF  - subscripts of the list of stars created by 
        APER which will be used to define the PSF.   Stars whose
        centroid does not fall within PSFRAD of the edge of the frame,
        or for which a Gaussian fit requires more than 25 iterations,
        will be ignored when creating the final PSF.
PSFRAD - the scalar radius, in pixels, of the circular area within
        which the PSF will be defined.   This should be slightly larger
        than the radius of the brightest star that one will be
        interested in.
FITRAD - the scalar radius, in pixels of the circular area used in the
        least-square star fits.  Stetson suggest that FITRAD should
        approximately equal to the FWHM, slightly less for crowded
        fields.  (FITRAD must be smaller than PSFRAD.)
PSFNAME- Name of the FITS file that will contain the table of residuals,
        and the best-fit Gaussian parameters.    This file is 
        subsequently required for use by NSTAR.  
OPTIONAL OUTPUTS
GAUSS  - 5 element vector giving parameters of gaussian fit to the 
        first PSF star
        GAUSS(0) - height of the gaussian (above sky)
        GAUSS(1) - the offset (in pixels) of the best fitting gaussian
                and the original X centroid
        GAUSS(2) - similiar offset from the Y centroid 
        GAUSS(3) - Gaussian sigma in X
        GAUSS(4) - Gaussian sigma in Y
PSF    - 2-d array of PSF residuals after a Gaussian fit.
PROCEDURE
GETPSF fits a Gaussian profile to the core of the first PSF star 
and generates a look-up table of the residuals of the
actual image data from the Gaussian fit.  If desired, it will then
fit this PSF to another star (using PKFIT) to determine its precise 
centroid, scale the same Gaussian to the new star's core, and add the
differences between the actual data and the scaled Gaussian to the
table of residuals.   (In other words, the Gaussian fit is performed
nly on the first star.)
OPTIONAL KEYWORD INPUT
DEBUG - if this keyword is set and non-zero, then the result of each
        fitting iteration will be displayed.
PROCEDURES CALLED
DAOERF, MAKE_2D, MKHDR, RINTER(), PKFIT, STRNUMBER(), STRN(), WRITEFITS
REVISON HISTORY
Adapted from the 1986 version of DAOPHOT in STSDAS
IDL Version 2  W Landsman           November 1988
Use DEBUG keyword instead of !DEBUG  W. Landsman       May 1996
Converted to IDL V5.0   W. Landsman   September 1997
eturn to caller
ave time in RINTER
nitialize the common blocks
nough parameters passed?
et number of rows and columns in image
otal # of stars identified in image
rray of PSF id's defined?
id User hit the [RETURN] key
o stellar ID's supplied
of stars used to create the PSF
s APMAG multidimensional?
Odd) width of box that contains PSF circle
ookup table has half pixel interpolation
Even) Width of subarray to be extracted from image
Find the first PSF star in the star list.
ounter for number of stars used to create PSF
D number of first PSF star
Now a subarray F will be read in from the big image, given by 
IXCEN-NBOX/2+1 <= x <= IXCEN+NBOX/2, IYCEN-NBOX/2+1 <= y <= IYCEN+NBOX/2.  
(NBOX is an even number.)  In the subarray, the coordinates of the centroid
of the star will lie between NBOX/2 and NBOX/2+1 in each coordinate.
pper & lower bounds in X
tar too close to edge?
ead in subarray, subtract off sky
An integrated Gaussian function will be fit to the central part of the
stellar profile. Initially, a 5x5 box centered on the centroid of the
star is used, but if the sigma in one coordinate drops to less than
1 pixel, then the box width of 3 will be used in that coordinate.
If the sigma increases to over 3 pixels, then a box width of 7 will be
used in that coordinate
coordinate of stellar centroid in subarray F
coordinate of stellar centroid in subarray F
ndex of pixel containing centroid
;Begin least squares
nitial guess for peak intensity
eginning of big iteration loop
Print the current star
egin the iterative loop
o convergence after 100 iterations?
default box width 
The T's are the first derivatives of the model profile with respect
to the five fitting parameters H, DXCEN, DYCEN, SIGX, and SIGY.
Note that the center of the best-fitting Gaussian profile is
expressed as an offset from the centroid of the star.  In the case of
a general, asymmetric stellar profile, the center of symmetry of the
best-fitting Gaussian profile will not necessarily coincide with the
centroid determined by any arbitrary centroiding algorithm.  
ubtract best fit Gaussian from subarray
DL version assumes INVERT is successful
ultiply by vector of residuals
orrect the fitting parameters
est for convergence
Now that the solution has converged, we can generate an
array containing the differences between the actual stellar profile
and the best-fitting Gaussian analytic profile.
arameters for Gaussian fit
ompute Gaussian
esiduals (Real profile - Gaussian)
The look-up table is obtained by interpolation within the array of
fitting residuals. We need to interpolate because we want the look-up
table to be centered accurately on the centroid of the star, which of
course is at some fractional-pixel position in the original data.
ndex function for PSF array
nterpolate residuals onto current star
ndex of first good star
For each additional star, determine the precise coordinates of the
centroid and the relative brightness of the star
by least-squares fitting to the current version of the point-spread
function. Then subtract off the appropriately scaled integral under
the analytic Gaussian function and add the departures of the actual
data from the analytic Gaussian function to the look-up table.
oop for additional PSF stars begins here
ave all the stars been done?
Fit the current version of the point-spread function to the data for
this star.
tar too close to edge?
onvergence in less than 25 iterations?
arameters of successful fit
ompute array of residuals
Values of the array of residuals are now interpolated to an NPSF by
NPSF (NPSF is an odd number) array centered on the centroid of the
star, and added to the existing look-up table of corrections to the
analytic profile
Now correct both the height of the analytic Gaussian, and the value
of the aperture-magnitude of the point-spread function for the
inclusion of the additional star.
Create FITS file containing the PSF created.
reate a minimal FITS header