Viewing contents of file '../idllib/uit/pro/circint.pro'
PRO circint, image, x0, y0, radius, tot, npix, totdif, npixdif, t8, $
mask=mask, verbose=verbose
;+
; NAME:
; CIRCINT
; PURPOSE:
; Integrate flux in circular apertures. User responsible for
; subtracting sky first. See notes (1) and (2) below.
; CALLING SEQUENCE:
; circint, image, x0, y0, radius, tot, npix, totdif, npixdif, t8, $
; mask=mask, /verbose
; INPUT PARAMETERS:
; image 2-D image array
; x0, y0 center of circular apertures (need not be integral)
; radius vector of aperture radii
; OUTPUT PARAMETERS (see note 2):
; tot total raw flux in each aperture
; npix number of pixels in each aperture
; totdif total raw flux for each annulus (difference of successive
; apertures)
; npixdif number of pixels in each annulus
; t8 like totdif, but each annulus is divided into eight sections
; (X dimension = 8, Y dimension = number of annuli). Allows
; assessment of error due to lumpiness of distribution, by
; method of Djorgovski & King
; OPTIONAL INPUT KEYWORDS:
; mask Must be same size as image. If a pixel in
; the mask is 1, the corresponding pixel in the image is
; counted. If 0, the corresponding pixel in the image is
; ignored (in all results). See note (3).
; verbose Use of /verbose will make the routine tell you what stage
; it's at. See note (1).
; COMMON BLOCKS:
; none
; NOTES:
; (1) If you aren't sure you've set up right, use the /verbose keyword,
; because the routine is fairly slow.
; (2) For a surface brightness profile, plot totdif/npixdif vs. radius.
; For an aperture growth curve, plot tot vs. radius.
; (3) Mask is intended to mask out stars or garbage. Depending on your
; application, you might be better off modifying the input image.
; SIDE EFFECTS: none
; PROCEDURE: Similar to IDL/UIT/DAOPHOT aperture routine.
; MODIFICATION HISTORY:
; Integrate in circular aperture.
; RSH - STX - 20 Aug 1990
; Modified to ignore some pixels according to mask image.
; RSH - STX - 19 Sept 90
; Totals done for annuli rather than discs.
; RSH - STX - 20 Sept 90
; Use of mask corrected. RSH - STX - 27 Sept 90
; Small change to conserve array space. RSH - STX - 5 Nov 90
; Fractional-pixel approximation adopted from Wayne Landsman's version
; of DAOPHOT APER. Annuli computed from discs. RSH - STX - 17 July 91
; Annuli divided into octants for subsequent error estimate due to random
; distribution of sources. RSH - STX - 22 July 91
; Several bugs fixed. RSH - STX - 3 Oct 91
; Spiffed up for UIT IDL library. RSH - Hughes STX - 11 June 92
; Speeded up initializaion. RSH - HSTX - 5 August 1993
;-
on_error,1
IF n_params(0) LT 1 THEN BEGIN
message,'Calling sequence: circint,image,x0,y0,' $
+'radius,tot,npix,totdif,npixdif,t8,' $
+'mask=mask,/verbose'
ENDIF
IF NOT keyword_set(verbose) THEN verbose=0
sz = size(image)
bell = string(07b)
IF keyword_set(mask) THEN BEGIN
szm = size(mask)
IF (szm(1) NE sz(1)) OR (szm(2) NE sz(2)) THEN BEGIN
message,'Mask image must be same size as flux image.'
ENDIF
ENDIF
IF verbose NE 0 THEN $
message,/inf,'Beginning routine to integrate in circular apertures ...'
nrad = n_elements(radius)
order_test = radius(1:nrad-1) - radius(0:nrad-2)
ww = where((order_test LE 0),count)
IF count GT 0 THEN $
message,'Radii must be specified in ascending order.'
ww = where((radius LT 0.5), count)
IF count GT 0 THEN $
message,'All radii must be 0.5 pixels or greater.'
maxrad = float(max(radius))
maxrad2 = maxrad^2
distance = replicate(2.0*maxrad2, sz(1), sz(2)) ; Like FLTARR with init
xdist2 = (findgen(sz(1)) - x0)^2
ydist2 = (findgen(sz(2)) - y0)^2
xmin = min(where(xdist2 LE maxrad2)) - 1
xmax = max(where(xdist2 LE maxrad2)) + 1
ymin = min(where(ydist2 LE maxrad2)) - 1
ymax = max(where(ydist2 LE maxrad2)) + 1
FOR ix=xmin,xmax DO FOR iy=ymin,ymax DO $
distance(ix,iy) = xdist2(ix) + ydist2(iy)
distance = sqrt(distance)-0.5
tot = fltarr(nrad)
totdif = tot
npix = fltarr(nrad)
npixdif = npix
; Sections for error estimates a la Djorgovski and King
t8 = fltarr(8,nrad)
octants = replicate(255b, sz(1), sz(2))
FOR ix=xmin,xmax DO BEGIN
FOR iy=ymin,ymax DO BEGIN
iyp = iy-y0 & ixp = ix-x0
IF (ixp GT 0) AND (iyp GE 0) AND (iyp LT ixp) THEN octants(ix,iy)=0b
IF (ixp GT 0) AND (iyp GE 0) AND (iyp GE ixp) THEN octants(ix,iy)=1b
IF (ixp LE 0) AND (iyp GT 0) AND (iyp GT -ixp) THEN octants(ix,iy)=2b
IF (ixp LE 0) AND (iyp GT 0) AND (iyp LE -ixp) THEN octants(ix,iy)=3b
IF (ixp LT 0) AND (iyp LE 0) AND (iyp GT ixp) THEN octants(ix,iy)=4b
IF (ixp LT 0) AND (iyp LE 0) AND (iyp LE ixp) THEN octants(ix,iy)=5b
IF (ixp GE 0) AND (iyp LT 0) AND (iyp LT -ixp) THEN octants(ix,iy)=6b
IF (ixp GE 0) AND (iyp LT 0) AND (iyp GE -ixp) THEN octants(ix,iy)=7b
octants(x0,y0)=99 ;center pixel not included an any octant
ENDFOR
ENDFOR
;
IF verbose NE 0 THEN message,/inf,'Initialization complete. Now integrating.'
FOR ir=0, nrad-1 DO BEGIN
IF verbose NE 0 THEN $
message,/inf,'Now doing aperture of ' $
+strtrim(radius(ir),2)+' pixel radius'
IF keyword_set(mask) THEN BEGIN
within = (distance LE radius(ir)) AND mask
ENDIF ELSE within = distance LE radius(ir)
thisapr = distance(where(within))
fractn = (radius(ir)-thisapr < 1.0 >0.0 ) ;Fraction of pixels to count
tot(ir) = total(image(where(within))*fractn)
npix(ir)= total(fractn)
FOR ioct = 0,7 DO BEGIN
section = where(within AND (octants EQ ioct))
thisoct = distance(section)
ofrac = (radius(ir)-thisoct < 1.0 > 0.0 )
IF verbose NE 0 THEN message,/inf, $
'Octant '+strtrim(ioct,2)+' with '+strtrim(n_elements(section),2) $
+' pixels'
t8(ioct,ir) = total(image(section)*ofrac)
ENDFOR
thisapr=0 & fractn=0 & within=0
ENDFOR
totdif (0) = tot(0)
npixdif(0) = npix(0)
FOR ir=1,nrad-1 DO BEGIN
npixdif(ir) = npix(ir) - npix(ir-1)
totdif(ir) = tot(ir) - tot(ir-1)
FOR ioct = 0,7 DO BEGIN
t8(ioct,nrad-ir) = t8(ioct,nrad-ir) - t8(ioct,nrad-ir-1)
ENDFOR
ENDFOR
IF verbose NE 0 THEN message,/inf,bell+'Done.'
RETURN
END