aper.pro

basic median filter simulation

pro aper,image,xc,yc,mags,errap,sky,skyerr,phpadu,apr,skyrad,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
;+
; 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, 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.
;     /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
;
;
;       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
;-
 COMPILE_OPT IDL2
 On_error,2
;             Set parameter limits
 minsky = 20   ;Smallest number of pixels from which the sky may be determined
 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=, ]'
  print,'             /SILENT, /NAN'
  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 not 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

;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