pro aper,image,xc,yc,mags,errap,sky,skyerr,phpadu,apr,skyradii,badpix, $ SETSKYVAL = setskyval,PRINT = print, SILENT = silent, FLUX=flux, $ EXACT = exact, Nan = nan, READNOISE = readnoise, MEANBACK = meanback, $ CLIPSIG=clipsig, MAXITER=maxiter,CONVERGE_NUM=converge_num, $ MINSKY = minsky ;+ ; 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, /NAN, /EXACT, /FLUX, PRINT = , /SILENT, ; /MEANBACK, MINSKY=, 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. Ignored if the SETSKYVAL ; keyword is set. ; BADPIX - Two element vector giving the minimum and maximum value ; of a good pixel. If badpix is not supplied or if BADPIX[0] is ; equal to BADPIX[1] then it is assumed that there are no bad ; pixels. Note that fluxes will not be computed for any star ; with a bad pixel within the aperture area, but that bad pixels ; will be simply ignored for the sky computation. The BADPIX ; parameter is ignored if the /NAN keyword is set. ; ; OPTIONAL KEYWORD INPUTS: ; CLIPSIG - if /MEANBACK is set, then this is the number of sigma at which ; to clip the background. Default=3 ; CONVERGE_NUM: if /MEANBACK is set then if the proportion of ; rejected pixels is less than this fraction, the iterations stop. ; Default=0.02, i.e., iteration stops if fewer than 2% of pixels ; excluded. ; /EXACT - By default, APER counts subpixels, but uses a polygon ; approximation for the intersection of a circular aperture with ; a square pixel (and normalizes 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. ; MAXITER if /MEANBACK is set then this is the ceiling on number of ; clipping iterations of the background. Default=5 ; /MEANBACK - if set, then the background is computed using the 3 sigma ; clipped mean (using meanclip.pro) rather than using the mode ; computed with mmm.pro. This keyword is useful for the Poisson ; count regime or where contamination is known to be minimal. ; MINSKY - Integer giving mininum number of sky values to be used with MMM ; APER will not compute a flux if fewer valid sky elements are ; within the sky annulus. Default = 20. ; /NAN - If set then APER will check for NAN values in the image. /NAN ; takes precedence over the BADPIX parameter. Note that fluxes ; will not be computed for any star with a NAN pixel within the ; aperture area, but that NAN pixels will be simply ignored for ; the sky computation. ; 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'. ; READNOISE - Scalar giving the read noise (or minimum noise for any ; pixel. This value is passed to the procedure mmm.pro when ; computing the sky, and is only need for images where ; the noise is low, and pixel values are quantized. ; /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). If the /FLUX keyword is not set, then ; a flux of 1 digital unit is assigned a zero point magnitude of ; 25. ; ERRAP - NAPER by NSTAR array giving error for each star. If a ; magnitude could not be determined then ERRAP = 9.99 (if in ; magnitudes) or ERRAP = !VALUES.F_NAN (if /FLUX is set). ; SKY - NSTAR element vector giving sky value for each star in ; flux units ; SKYERR - NSTAR element vector giving error in sky values ; ; EXAMPLE: ; Determine the flux and error for photometry radii of 3 and 5 pixels ; surrounding the position 234.2,344.3 on an image array, im. Compute ; the partial pixel area exactly. Assume that the flux units are in ; Poisson counts, so that PHPADU = 1, and the sky value is already known ; to be 1.3, and that the range [-32767,80000] for bad low and bad high ; pixels ; ; ; IDL> aper, im, 234.2, 344.3, flux, eflux, sky,skyerr, 1, [3,5], -1, $ ; [-32767,80000],/exact, /flux, setsky = 1.3 ; ; PROCEDURES USED: ; 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 ; (5) The total computed flux is negative. In this case the negative ; flux and error are returned. ; ; ; For the case where the source is fainter than the background, APER will ; return negative fluxes if /FLUX is set, but will otherwise give ; invalid data (since negative fluxes can't be converted to magnitudes) ; ; 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 ; 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 ; Don't abort for insufficient sky pixels W. Landsman May 2000 ; Added /EXACT keyword W. Landsman June 2000 ; Allow SETSKYVAL = 0 W. Landsman December 2000 ; Set BADPIX[0] = BADPIX[1] to ignore bad pixels W. L. January 2001 ; Fix chk_badpixel problem introduced Jan 01 C. Ishida/W.L. February 2001 ; Set bad fluxes and error to NAN if /FLUX is set W. Landsman Oct. 2001 ; Remove restrictions on maximum sky radius W. Landsman July 2003 ; Added /NAN keyword W. Landsman November 2004 ; Set badflux=0 if neither /NAN nor badpix is set M. Perrin December 2004 ; Added READNOISE keyword W. Landsman January 2005 ; Added MEANBACK keyword W. Landsman October 2005 ; Correct typo when /EXACT and multiple apertures used. W.L. Dec 2005 ; Remove VMS-specific code W.L. Sep 2006 ; Add additional keywords if /MEANBACK is set W.L Nov 2006 ; Allow negative fluxes if /FLUX is set W.L. Mar 2008 ; Previous update would crash if first star was out of range W.L. Mar 2008 ; Fix floating equality test for bad magnitudes W.L./J.van Eyken Jul 2009 ; Added MINSKY keyword W.L. Dec 2011 ; Don't ever modify input skyrad variable W. Landsman Aug 2013 ; Avoid integer overflow for very big images W. Landsman/R. Gutermuth Mar 2016 ; Eliminate limit on maximum number of sky pixels W. Landsman Dec 2016 ;- COMPILE_OPT IDL2 On_error,2 ; Set parameter limits ;Smallest number of pixels from which the sky may be determined if N_elements(minsky) EQ 0 then minsky = 20 ; 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=, ]' print,' /SILENT, /NAN, MINSKY=' 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 ~keyword_set(nan) then begin 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 [none]: ',ans if ans EQ '' then badpix = [0,0] 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 chk_badpix = badpix[0] LT badpix[1] ;Ignore bad pixel checks? 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 else begin if N_elements(skyradii) NE 2 then begin skyrad = fltarr(2) read,'Enter inner and outer sky radius (pixel units): ',skyrad endif else skyrad = float(skyradii) endelse 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,1x,3A,3(/,28x,4A,:))' ;Screen format fmt3 = '(I4,5F8.2,1x,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 if keyword_set(EXACT) then begin bigrad = apr + 0.5 smallrad = apr/sqrt(2) - 0.5 endif if N_elements(SETSKYVAL) EQ 0 then begin rinsq = (skyrad[0]> 0.)^2 routsq = skyrad[1]^2 endif if keyword_set(PRINT) then begin ;Open output file and write header info? if size(PRINT,/TNAME) NE 'STRING' 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'),' Host: ',getenv('HOST') for j = 0, Naper-1 do printf,lun, $ format='(a,i2,a,f4.1)','Radius of aperture ',j,' = ',apr[j] if N_elements(SETSKYVAL) EQ 0 then begin 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] endif else printf,lun,'Sky values fixed at ', strtrim(setskyval[0],2) 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 ~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 = long(xc-skyrad[1]) > 0 ;Lower limit X direction ux = long(xc+skyrad[1]) < (ncol-1) ;Upper limit X direction nx = ux-lx+1 ;Number of pixels X direction ly = long(yc-skyrad[1]) > 0 ;Lower limit Y direction uy = long(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' if keyword_set(flux) then begin badval = !VALUES.F_NAN baderr = badval endif else begin badval = 99.999 baderr = 9.999 endelse for i = 0L, Nstars-1 do begin ;Compute magnitudes for each star apmag = replicate(badval, Naper) & magerr = replicate(baderr, Naper) skymod = 0. & skysig = 0. & skyskw = 0. ;Sky mode sigma and skew if badstar[i] then goto, BADSTAR error1=apmag & error2 = apmag & error3 = apmag 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 N_elements(SETSKYVAL) EQ 0 then begin skypix = ( rsq GE rinsq ) and ( rsq LE routsq ) if keyword_set(nan) then skypix = skypix and finite(rotbuf) $ else if chk_badpix then skypix = skypix and ( rotbuf GT badpix[0] ) and $ (rotbuf LT badpix[1] ) sindex = where(skypix, Nsky) if ( nsky LT minsky ) then begin ;Sufficient sky pixels? if ~silent then $ message,'There aren''t enough valid pixels in the sky annulus.',/con goto, BADSTAR endif skybuf = rotbuf[ sindex[0:nsky-1] ] if keyword_set(meanback) then $ meanclip,skybuf,skymod,skysig, $ CLIPSIG=clipsig, MAXITER=maxiter, CONVERGE_NUM=converge_num else $ mmm, skybuf, skymod, skysig, skyskw, readnoise=readnoise,minsky=minsky ; 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 the modal sky value could not be determined, then all apertures for this ; star are bad if ( skysig LT 0.0 ) then goto, BADSTAR 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] GE apr[k] ) then begin ;Does aperture extend outside the image? 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[k] )) ;Fix 05-Dec-05 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 keyword_set(NaN) then $ badflux = min(finite(thisapd)) EQ 0 $ else if chk_badpix then begin minthisapd = min(thisapd, max = maxthisapd) badflux = (minthisapd LE badpix[0] ) or ( maxthisapd GE badpix[1]) endif else badflux = 0 if ~badflux then $ apmag[k] = total(thisapd*fractn) ;Total over irregular aperture endif endfor ;k if keyword_set(flux) then g = where(finite(apmag), Ng) else $ g = where(abs(apmag - badval) GT 0.01, Ng) if Ng GT 0 then begin apmag[g] = apmag[g] - skymod*area[g] ;Subtract sky from the integrated brightnesses ; Add in quadrature 3 sources of error: (1) random noise inside the star ; aperture, including readout noise and the degree of contamination by other ; stars in the neighborhood, as estimated by the scatter in the sky values ; (this standard error increases as the square root of the area of the ; aperture); (2) the Poisson statistics of the observed star brightness; ; (3) the uncertainty of the mean sky brightness (this standard error ; increases directly with the area of the aperture). error1[g] = area[g]*skyvar ;Scatter in sky values error2[g] = (apmag[g] > 0)/phpadu ;Random photon noise error3[g] = sigsq*area[g]^2 ;Uncertainty in mean sky brightness magerr[g] = sqrt(error1[g] + error2[g] + error3[g]) if ~keyword_set(FLUX) then begin good = where (apmag GT 0.0, Ngood) ;Are there any valid integrated fluxes? if ( Ngood GT 0 ) then begin ;If YES then compute errors magerr[good] = 1.0857*magerr[good]/apmag[good] ;1.0857 = log(10)/2.5 apmag[good] = 25.-2.5*alog10(apmag[good]) endif endif endif BADSTAR: ;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 ~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