Viewing contents of file '../idllib/uit/pro/aper.pro'
pro aper,image,xc,yc,mags,errap,sky,skyerr,phpadu,apr,skyrad,badpix, $
TEXTOUT=textout, SILENT = silent
;+
; NAME:
; APER
; PURPOSE:
; Procedure (adapted from DAOPHOT) to compute concentric aperture
; photometry. A separate sky value is determined for each source
; using specified inner and outer sky radii.
;
; CALLING SEQUENCE:
; aper, image, xc, yc, [ mags, errap, sky, skyerr, phpadu, apr, skyrad,
; badpix, TEXTOUT = , SILENT = ]
; INPUTS:
; IMAGE - input image array
; XC - vector of x coordinates. COORDINATE VALUES
; MUST BE IN IDL CONVENTION! (First pixel is 0 not 1!)
; 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:
; TEXTOUT - Determines output device (see TEXTOPEN)
; textout=1 TERMINAL using /more option
; textout=2 TERMINAL without /more option
; textout=3 <program>.prt
; textout=4 laser.tmp
; textout=5 user must open file
; textout = filename (default extension of .prt)
; SILENT - If supplied and non-zero then no output is displayed to the
; terminal.
;
; 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
;
; SYSTEM VARIABLES:
; If keyword textout is not used, the non-standard system variable
; !TEXTOUT will be used to determine the output device.
; PROCEDURES USED:
; MMM, TEXTOPEN, STRN
; 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
;-
On_error,2
; Set parameter limits
minsky = 20 ;Smallest number of pixels from which the sky may be determined
maxrad = 35. ;Maximum outer radius permitted for the sky annulus.
maxsky = 4000 ;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, TEXTOUT =, SILENT = ]'
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
if not keyword_set(TEXTOUT) then textout = !TEXTOUT
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
;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(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(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
fmt = '(F7.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,6A,:))' ;Print format
mags = fltarr( Naper, Nstars) & errap = mags ;Declare arrays
sky = fltarr( Nstars ) & skyerr = sky
area =fltarr( Naper ) & apmag=area & magerr = area
error1=area & error2 = area & error3 = area
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
if textout NE 1 then begin ;Open output file and write header info?
textopen, 'APER', TEXTOUT=textout
printf,!TEXTUNIT,' Program: APER '+strmid(systime(),0,20), ' User: ', $
getenv('USER'),' Node: ',getenv('NODE')
for j = 0, Naper-1 do printf,!TEXTUNIT, $
format='(a,i2,a,f4.1)','Radius of aperture ',j,' = ',apr(j)
printf,!TEXTUNIT,f='(/a,f4.1)','Inner radius for sky annulus = ',skyrad(0)
printf,!TEXTUNIT,f='(a,f4.1)', 'Outer radius for sky annulus = ',skyrad(1)
printf,!TEXTUNIT,f='(/a)', $
'STAR X Y SKY SKYSIG SKYSKW MAGNITUDES'
endif
; Print header
if not SILENT then print, $
format="(/1X,'STAR',5X,'X',7X,'Y',6X,'SKY',8X,'MAGNITUDES')"
rinsq = (skyrad(0)> 0.)^2
routsq = skyrad(1)^2
; 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 = 0, 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) )
for ii = 0, ny(i)-1 do rsq(0,ii) = dxsq + (ii-dy(i))^2
; Select pixels within sky annulus, and eliminate pixels falling
; below BADLO threshold. SKYBUF will be 1-d array of sky pixels
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?
print,'APER: There aren''t enough pixels in the sky annulus.'
print,' Are you sure your bad pixel threshold is all right?'
print,' If so, then you need a larger outer sky radius.'
return
endif
; Obtain the mode, standard deviation, and skewness of the peak in the
; sky histogram, by calling MMM.
mmm, skybuf, skymod, skysig, skyskw
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
r = sqrt(rsq) - 0.5 ;2-d array of the radius of each pixel in the subarray
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
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
; 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
area(k) = total(fractn)
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))
magerr(good) = 1.0857*magerr(good)/apmag(good) ;1.0857 = log(10)/2.5
apmag(good) = 25.-2.5*alog10(apmag(good))
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 ( textout GT 2 ) then printf,!TEXTUNIT, $ ;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 textout NE 1 then textclose, TEXTOUT=textout ;Close output file
return
end