Viewing contents of file '../idllib/astron/pro/aper.pro'
pro aper,image,xc,yc,mags,errap,sky,skyerr,phpadu,apr,skyrad,badpix, $
SETSKYVAL = setskyval,PRINT = print, SILENT = silent, FLUX=flux, $
EXACT = exact
;+
; NAME:
; APER
; PURPOSE:
; Compute concentric aperture photometry (adapted from DAOPHOT)
; EXPLANATION:
; APER can compute photometry in several user-specified aperture radii.
; A separate sky value is computed for each source using specified inner
; and outer sky radii.
;
; CALLING SEQUENCE:
; APER, image, xc, yc, [ mags, errap, sky, skyerr, phpadu, apr, skyrad,
; badpix, /EXACT, /FLUX, PRINT = , /SILENT, SETSKYVAL = ]
; INPUTS:
; IMAGE - input image array
; XC - vector of x coordinates.
; YC - vector of y coordinates
;
; OPTIONAL INPUTS:
; PHPADU - Photons per Analog Digital Units, numeric scalar. Converts
; the data numbers in IMAGE to photon units. (APER assumes
; Poisson statistics.)
; APR - Vector of up to 12 REAL photometry aperture radii.
; SKYRAD - Two element vector giving the inner and outer radii
; to be used for the sky annulus
; BADPIX - Two element vector giving the minimum and maximum value
; of a good pix (Default [-32765,32767])
;
; OPTIONAL KEYWORD INPUTS:
; /EXACT - By default, APER counts subpixels, but uses a polygon
; approximation for the intersection of a circular aperture with
; a square pixel (and normalize the total area of the sum of the
; pixels to exactly match the circular area). If the /EXACT
; keyword, then the intersection of the circular aperture with a
; square pixel is computed exactly. The /EXACT keyword is much
; slower and is only needed when small (~2 pixels) apertures are
; used with very undersampled data.
; /FLUX - By default, APER uses a magnitude system where a magnitude of
; 25 corresponds to 1 flux unit. If set, then APER will keep
; results in flux units instead of magnitudes.
; PRINT - if set and non-zero then APER will also write its results to
; a file aper.prt. One can specify the output file name by
; setting PRINT = 'filename'.
; /SILENT - If supplied and non-zero then no output is displayed to the
; terminal.
; SETSKYVAL - Use this keyword to force the sky to a specified value
; rather than have APER compute a sky value. SETSKYVAL
; can either be a scalar specifying the sky value to use for
; all sources, or a 3 element vector specifying the sky value,
; the sigma of the sky value, and the number of elements used
; to compute a sky value. The 3 element form of SETSKYVAL
; is needed for accurate error budgeting.
;
; OUTPUTS:
; MAGS - NAPER by NSTAR array giving the magnitude for each star in
; each aperture. (NAPER is the number of apertures, and NSTAR
; is the number of stars). A flux of 1 digital unit is assigned
; a zero point magnitude of 25.
; ERRAP - NAPER by NSTAR array giving error in magnitude
; for each star. If a magnitude could not be deter-
; mined then ERRAP = 9.99.
; SKY - NSTAR element vector giving sky value for each star
; SKYERR - NSTAR element vector giving error in sky values
;
; PROCEDURES USED:
; DATATYPE(), GETOPT, MMM, PIXWT(), STRN(), STRNUMBER()
; NOTES:
; Reasons that a valid magnitude cannot be computed include the following:
; (1) Star position is too close (within 0.5 pixels) to edge of the frame
; (2) Less than 20 valid pixels available for computing sky
; (3) Modal value of sky could not be computed by the procedure MMM
; (4) *Any* pixel within the aperture radius is a "bad" pixel
;
; APER was modified in June 2000 in two ways: (1) the /EXACT keyword was
; added (2) the approximation of the intersection of a circular aperture
; with square pixels was improved (i.e. when /EXACT is not used)
; REVISON HISTORY:
; Adapted to IDL from DAOPHOT June, 1989 B. Pfarr, STX
; Adapted for IDL Version 2, J. Isensee, July, 1990
; Code, documentation spiffed up W. Landsman August 1991
; TEXTOUT may be a string W. Landsman September 1995
; FLUX keyword added J. E. Hollis, February, 1996
; SETSKYVAL keyword, increase maxsky W. Landsman, May 1997
; Work for more than 32767 stars W. Landsman, August 1997
; Converted to IDL V5.0 W. Landsman September 1997
; Don't abort for insufficient sky pixels W. Landsman May 2000
; Added /EXACT keyword W. Landsman June 2000
;-
On_error,2
; Set parameter limits
minsky = 20 ;Smallest number of pixels from which the sky may be determined
maxrad = 100. ;Maximum outer radius permitted for the sky annulus.
maxsky = 10000 ;Maximum number of pixels allowed in the sky annulus.
;
if N_params() LT 3 then begin ;Enough parameters supplied?
print, $
'Syntax - APER, image, xc, yc, [ mags, errap, sky, skyerr, phpadu, apr, '
print,' skyrad, badpix, /EXACT, /FLUX, SETSKYVAL = ,PRINT =, ]'
return
endif
s = size(image)
if ( s[0] NE 2 ) then message, $
'ERROR - Image array (first parameter) must be 2 dimensional'
ncol = s[1] & nrow = s[2] ;Number of columns and rows in image array
silent = keyword_set(SILENT)
if N_elements(badpix) NE 2 then begin ;Bad pixel values supplied
GET_BADPIX:
ans = ''
print,'Enter low and high bad pixel values, [RETURN] for defaults'
read,'Low and high bad pixel values [-32765,32767]: ',ans
if ans EQ '' then badpix = [-32765,32767] else begin
badpix = getopt(ans,'F')
if ( N_elements(badpix) NE 2 ) then begin
message,'Expecting 2 scalar values',/continue
goto,GET_BADPIX
endif
endelse
endif
if ( N_elements(apr) LT 1 ) then begin ;Read in aperture sizes?
apr = fltarr(10)
read, 'Enter first aperture radius: ',ap
apr[0] = ap
ap = 'aper'
for i = 1,9 do begin
GETAP:
read,'Enter another aperture radius, [RETURN to terminate]: ',ap
if ap EQ '' then goto,DONE
result = strnumber(ap,val)
if result EQ 1 then apr[i] = val else goto, GETAP
endfor
DONE:
apr = apr[0:i-1]
endif
if N_elements(SETSKYVAL) GT 0 then begin
if N_elements( SETSKYVAL ) EQ 1 then setskyval = [setskyval,0.,1.]
if N_elements( SETSKYVAL ) NE 3 then message, $
'ERROR - Keyword SETSKYVAL must contain 1 or 3 elements'
skyrad = [ 0., max(apr) + 1]
endif
;Get radii of sky annulii
if N_elements(skyrad) NE 2 then begin
skyrad = fltarr(2)
read,'Enter inner and outer sky radius (pixel units): ',skyrad
endif else skyrad = float(skyrad)
if ( N_elements(phpadu) LT 1 ) then $
read,'Enter scale factor in Photons per Analog per Digital Unit: ',phpadu
Naper = N_elements( apr ) ;Number of apertures
Nstars = min([ N_elements(xc), N_elements(yc) ]) ;Number of stars to measure
ms = strarr( Naper ) ;String array to display mag for each aperture
if keyword_set(flux) then $
fmt = '(F8.1,1x,A,F7.1)' else $ ;Flux format
fmt = '(F9.3,A,F5.3)' ;Magnitude format
fmt2 = '(I5,2F8.2,F7.2,3A,3(/,28x,4A,:))' ;Screen format
fmt3 = '(I4,5F8.2,6A,2(/,44x,9A,:))' ;Print format
mags = fltarr( Naper, Nstars) & errap = mags ;Declare arrays
sky = fltarr( Nstars ) & skyerr = sky
area = !PI*apr*apr ;Area of each aperture
apmag= fltarr(Naper) & magerr = apmag
error1=apmag & error2 = apmag & error3 = apmag
if keyword_set(EXACT) then begin
bigrad = apr + 0.5
smallrad = apr/sqrt(2) - 0.5
endif
if not keyword_set(SETSKYVAL) then begin
if skyrad[1] GT maxrad then begin
message,/INF,'WARNING - Outer sky radius being reset to '+ strn(maxrad)
skyrad[1] = (maxrad-0.001) ;Outer sky radius less than MAXRAD?
endif
rinsq = (skyrad[0]> 0.)^2
routsq = skyrad[1]^2
endif
if keyword_set(PRINT) then begin ;Open output file and write header info?
if datatype(PRINT) NE 'STR' then file = 'aper.prt' $
else file = print
message,'Results will be written to a file ' + file,/INF
openw,lun,file,/GET_LUN
printf,lun,' Program: APER '+ systime(), ' User: ', $
getenv('USER'),' Node: ',getenv('NODE')
for j = 0, Naper-1 do printf,lun, $
format='(a,i2,a,f4.1)','Radius of aperture ',j,' = ',apr[j]
printf,lun,f='(/a,f4.1)','Inner radius for sky annulus = ',skyrad[0]
printf,lun,f='(a,f4.1)', 'Outer radius for sky annulus = ',skyrad[1]
if keyword_set(FLUX) then begin
printf,lun,f='(/a)', $
'STAR X Y SKY SKYSIG SKYSKW FLUXES'
endif else printf,lun,f='(/a)', $
'STAR X Y SKY SKYSIG SKYSKW MAGNITUDES'
endif
print = keyword_set(PRINT)
; Print header
if not SILENT then begin
if (KEYWORD_SET(FLUX)) then begin
print, format="(/1X,'STAR',5X,'X',7X,'Y',6X,'SKY',8X,'FLUXES')"
endif else print, $
format="(/1X,'STAR',5X,'X',7X,'Y',6X,'SKY',8X,'MAGNITUDES')"
endif
; Compute the limits of the submatrix. Do all stars in vector notation.
lx = fix(xc-skyrad[1]) > 0 ;Lower limit X direction
ux = fix(xc+skyrad[1]) < (ncol-1) ;Upper limit X direction
nx = ux-lx+1 ;Number of pixels X direction
ly = fix(yc-skyrad[1]) > 0 ;Lower limit Y direction
uy = fix(yc+skyrad[1]) < (nrow-1); ;Upper limit Y direction
ny = uy-ly +1 ;Number of pixels Y direction
dx = xc-lx ;X coordinate of star's centroid in subarray
dy = yc-ly ;Y coordinate of star's centroid in subarray
edge = (dx-0.5) < (nx+0.5-dx) < (dy-0.5) < (ny+0.5-dy) ;Closest edge to array
badstar = ((xc LT 0.5) or (xc GT ncol-1.5) $ ;Stars too close to the edge
or (yc LT 0.5) or (yc GT nrow-1.5))
;
badindex = where( badstar, Nbad) ;Any stars outside image
if ( Nbad GT 0 ) then message, /INF, $
'WARNING - ' + strn(nbad) + ' star positions outside image'
for i = 0L, Nstars-1 do begin ;Compute magnitudes for each star
skymod = 0. & skysig = 0. & skyskw = 0. ;Sky mode sigma and skew
if badstar[i] then begin ;
apmag[*] = -1.0E-36
goto, BADSTAR
endif
rotbuf = image[ lx[i]:ux[i], ly[i]:uy[i] ] ;Extract subarray from image
; RSQ will be an array, the same size as ROTBUF containing the square of
; the distance of each pixel to the center pixel.
dxsq = ( findgen( nx[i] ) - dx[i] )^2
rsq = fltarr( nx[i], ny[i], /NOZERO )
for ii = 0, ny[i]-1 do rsq[0,ii] = dxsq + (ii-dy[i])^2
if keyword_set(exact) then begin
nbox = lindgen(nx[i]*ny[i])
xx = reform( (nbox mod nx[i]), nx[i], ny[i])
yy = reform( (nbox/nx[i]),nx[i],ny[i])
x1 = abs(xx-dx[i])
y1 = abs(yy-dy[i])
endif else begin
r = sqrt(rsq) - 0.5 ;2-d array of the radius of each pixel in the subarray
endelse
; Select pixels within sky annulus, and eliminate pixels falling
; below BADLO threshold. SKYBUF will be 1-d array of sky pixels
if not keyword_set (SETSKYVAL) then begin
sindex = where( ( rsq GE rinsq ) and ( rsq LE routsq ) and $
( rotbuf GT badpix[0] ) )
nsky = N_elements(sindex) < maxsky ;Must be less than MAXSKY pixels
skybuf = rotbuf[ sindex[0:nsky-1] ]
if ( nsky LT minsky ) then begin ;Sufficient sky pixels?
if not silent then $
message,'There aren''t enough valid pixels in the sky annulus.',/con
mags[*,i] = -1.0e36
goto, BADSTAR
endif else mmm, skybuf, skymod, skysig, skyskw
; Obtain the mode, standard deviation, and skewness of the peak in the
; sky histogram, by calling MMM.
skyvar = skysig^2 ;Variance of the sky brightness
sigsq = skyvar/nsky ;Square of standard error of mean sky brightness
if ( skysig LT 0.0 ) then begin ;If the modal sky value could not be
apmag[*] = -99.999 ;determined, then all apertures for
goto, BADSTAR ;this star are bad.
endif
skysig = skysig < 999.99 ;Don't overload output formats
skyskw = skyskw >(-99)<999.9
endif else begin
skymod = setskyval[0]
skysig = setskyval[1]
nsky = setskyval[2]
skyvar = skysig^2
sigsq = skyvar/nsky
skyskw = 0
endelse
for k = 0,Naper-1 do begin ;Find pixels within each aperture
if ( edge[i] LT apr[k] ) then $ ;Does aperture extend outside the image?
apmag[k] = -1.0E36 $
else begin
if keyword_set(EXACT) then begin
mask = fltarr(nx[i],ny[i])
good = where( ( x1 LT smallrad[k] ) and (y1 LT smallrad[k] ), Ngood)
if Ngood GT 0 then mask(good) = 1.0
bad = where( (x1 GT bigrad[k]) or (y1 GT bigrad ))
mask(bad) = -1
gfract = where(mask EQ 0.0, Nfract)
if Nfract GT 0 then mask[gfract] = $
PIXWT(dx[i],dy[i],apr[k],xx[gfract],yy[gfract]) > 0.0
thisap = where(mask GT 0.0)
thisapd = rotbuf[thisap]
fractn = mask[thisap]
endif else begin
;
thisap = where( r LT apr[k] ) ;Select pixels within radius
thisapd = rotbuf[thisap]
thisapr = r[thisap]
fractn = (apr[k]-thisapr < 1.0 >0.0 ) ;Fraction of pixels to count
full = fractn EQ 1.0
gfull = where(full, Nfull)
gfract = where(1 - full)
factor = (area[k] - Nfull ) / total(fractn[gfract])
fractn[gfract] = fractn[gfract]*factor
endelse
; If the pixel is bad, set the total counts in this aperture to a large
; negative number
;
if ( (min(thisapd) LE badpix[0] ) or ( max(thisapd) GE badpix[1]) ) then $
apmag[k] = -1.0E36 $
else begin
apmag[k] = total(thisapd*fractn) ;Total over irregular aperture
endelse
endelse
endfor ;k
apmag = apmag - skymod*area ;Subtract sky from the integrated brightnesses
good = where (apmag GT 0.0, Ngood) ;Are there any valid integrated fluxes?
if ( Ngood GT 0 ) then begin ;If YES then compute errors
error1[good] = area[good]*skyvar ;Scatter in sky values
error2[good] = apmag[good]/phpadu ;Random photon noise
error3[good] = sigsq*area[good]^2 ;Uncertainty in mean sky brightness
magerr[good] = sqrt(error1[good] + error2[good] + error3[good])
if not keyword_set(FLUX) then begin
magerr[good] = 1.0857*magerr[good]/apmag[good] ;1.0857 = log(10)/2.5
apmag[good] = 25.-2.5*alog10(apmag[good])
endif
endif
BADSTAR:
;Assign fluxes to bad stars
nogood = where (apmag LE 0.0, Nbad)
if ( nbad GT 0 ) then begin
apmag[nogood] = 99.999
magerr[nogood] = 9.999
endif
;Print out magnitudes for this star
for ii = 0,Naper-1 do $ ;Concatenate mags into a string
ms[ii] = string( apmag[ii],'+-',magerr[ii], FORM = fmt)
if PRINT then printf,lun, $ ;Write results to file?
form = fmt3, i, xc[i], yc[i], skymod, skysig, skyskw, ms
if not SILENT then print,form = fmt2, $ ;Write results to terminal?
i,xc[i],yc[i],skymod,ms
sky[i] = skymod & skyerr[i] = skysig ;Store in output variable
mags[0,i] = apmag & errap[0,i]= magerr
endfor ;i
if PRINT then free_lun, lun ;Close output file
return
end