aper.pro

basic median filter simulation

 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

 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 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'
      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) 
 Nsky =   Nsky < maxsky   ;Must be less than MAXSKY pixels
 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
    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
           
 

;  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 not 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(apmag NE badval, Ng)
   if Ng GT 0 then begin 
  apmag[g] = apmag[g] - skymod*area[g]  ;Subtract sky from the integrated brightnesses

   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 not 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 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