nstar.pro

basic median filter simulation

      xold = dblarr(nterm) 
      clamp = replicate(1.,nterm)               ;Release all clamps
      clip = 0b
      niter = niter-1                           ;Back up iteration counter
      goto, RESTART 
      endif
   endfor
   endif
   j = j+1
 endwhile

 xpsfmin = (fix (xg - psfrad+1)) > 0
 xpsfmax = (fix (xg + psfrad  )) < (nx-1)
 ypsfmin = (fix (yg - psfrad+1)) > 0
 ypsfmax = (fix (yg + psfrad  )) < (ny-1)
 npsfx = xpsfmax-xpsfmin+1 & npsfy = ypsfmax-ypsfmin+1
 wt = fltarr(nx,ny)
 mask = bytarr(nx,ny)
 nterm = 3*nstr + varsky
 chi = fltarr(nstr) & sumwt = chi & numer = chi & denom = chi
 c = fltarr(nterm,nterm) & v = fltarr(nterm)

 for j = 0,nstr-1 do begin   ;Mask of pixels within fitting radius of any star
     x1 = xfitmin[j]  &  y1 = yfitmin[j]
     x2 = xfitmax[j]  &  y2 = yfitmax[j]
     rpixsq = fltarr(nfitx[j],nfity[j])
     xfitgen2 = (findgen(nfitx[j]) + x1 - xg[j])^2
     yfitgen2 = (findgen(nfity[j]) + y1 - yg[j])^2
     for k=0,nfity[j]-1 do rpixsq[0,k] = xfitgen2 + yfitgen2[k]
     temp = (rpixsq LE 0.999998*radsq)
     mask[x1,y1] = mask[x1:x2,y1:y2] or temp
     good = where(temp)
     rsq = rpixsq[good]/radsq
     temp1 = wt[x1:x2,y1:y2] 
     temp1[good] = temp1[good] > (5./(5.+rsq/(1.-rsq)) )
     wt[x1,y1] = temp1
 endfor

 igood = where(mask, ngoodpix)
 x = dblarr(ngoodpix,nterm)
 if varsky then x[0, nterm-1] = replicate(-1.0d, ngoodpix)

 psfmask = bytarr(ngoodpix,nstr)
 d = dimage[igood] - skybar
 for j = 0,nstr-1 do begin ;Masks of pixels within PSF radius of each star
     x1 = xpsfmin[j]   &    y1 = ypsfmin[j]
     x2 = xpsfmax[j]   &    y2 = ypsfmax[j]
     xgen = lindgen(npsfx[j]) + x1 - xg[j]
     ygen = lindgen(npsfy[j]) + y1 - yg[j]
     xgen2 = xgen^2 & ygen2 = ygen^2
     rpxsq = fltarr( npsfx[j],npsfy[j] )
     for k = 0,npsfy[j]-1 do rpxsq[0,k] = xgen2 + ygen2[k]
     temp =  mask[x1:x2,y1:y2] and (rpxsq LT psfrsq)
     temp1 = bytarr(nx,ny)
     temp1[x1,y1] = temp 
     goodfit = where(temp1[igood])
     psfmask[goodfit+ngoodpix*j] = 1b
     good = where(temp)
     xgood = xgen[good mod npsfx[j]] & ygood = ygen[good/npsfx[j]]
     model = dao_value(xgood,ygood,gauss,psf,dvdx,dvdy)
     d[goodfit] = d[goodfit] - magg[j]*model
     x[goodfit + 3*j*ngoodpix] = -model
     x[goodfit + (3*j+1)*ngoodpix] = magg[j]*dvdx
     x[goodfit + (3*j+2)*ngoodpix] = magg[j]*dvdy
 endfor

 wt = wt[igood] & idimage = dimage[igood]
 dpos = (idimage-d) > 0
 sigsq = dpos/phpadu + ronois + (0.0075*dpos)^2 + (pkerr*(dpos-skybar))^2

 relerr = abs(d)/sqrt(sigsq)
 if clip then begin   ;Reject pixels with 20 sigma errors (after 1st iteration)
        bigpix = where(relerr GT 20.*chiold, nbigpix)
        if ( nbigpix GT 0 ) then begin
              keep = indgen(ngoodpix)
              for i = 0,nbigpix-1 do keep = keep[ where( keep NE bigpix[i]) ]
              wt= wt[keep] & d = d[keep] & idimage = idimage[keep] 
              igood= igood[keep]  & relerr = relerr[keep]
              psfmask = psfmask[keep,*]   &  x = x[keep,*]
        endif
 endif

 sumres = total(relerr*wt)
 grpwt = total(wt)

 dpos = ((idimage-skybar) > 0) + skybar
 sig = dpos/phpadu + ronois + (0.0075*dpos)^2 + (pkerr*(dpos-skybar))^2
 for j = 0,nstr-1 do begin
     goodfit = where(psfmask[*,j])
     chi[j] = total(relerr[goodfit]*wt[goodfit])
     sumwt[j] = total(wt[goodfit])
     xgood = igood[goodfit] mod nx & ygood = igood[goodfit]/nx
     rhosq = ((xg[j] - xgood)/gauss[3])^2  +  ((yg[j] - ygood)/gauss[4])^2
     goodsig = where(rhosq LT 36)     ;Include in sharpness index only
     rhosq = 0.5*rhosq[goodsig]       ;pixels within 6 sigma of centroid
     dfdsig = exp(-rhosq)*(rhosq-1.)
     sigpsf = sig[goodfit[goodsig]] & dsig = d[goodfit[goodsig]]
     numer[j] = total(dfdsig*dsig/sigpsf)
     denom[j] = total(dfdsig^2/sigpsf)
 endfor

 wt = wt/sigsq
 if clip then $  ;After 1st iteration, reduce weight of a bad pixel
   wt = wt/(1.+(0.4*relerr/chiold)^8) 

 v = d * wt # x
 c = transpose(x) # ( ( wt # replicate(1.,nterm) ) * x )

 if grpwt GT 3 then begin
        chiold = 1.2533*sumres*sqrt(1./(grpwt*(grpwt-3.)))
        chiold = ((grpwt-3.)*chiold+3.)/grpwt
 endif

 i = where(sumwt GT 3, ngood)
 if ngood GT 0 then begin 
     chi[i] = 1.2533*chi[i]*sqrt(1./((sumwt[i]-3.)*sumwt[i]))
     chi[i] = ((sumwt[i]-3.)*chi[i]+3.)/sumwt[i]
 endif

chibad = where(sumwt LE 3, ngood)
if ngood GT 0 then chi[chibad] = chiold

 c = invert(c)
 x = c # transpose(v)

 if (not clip) or (niter LE 1) then redo = 1b else redo = 0b 
 if varsky then begin
        skybar = skybar - x[nterm-1]
        if abs(x[nterm-1]) GT  0.01 then redo = 1b
 endif
 clip = 1b

 j = 3*indgen(nstr) & k = j+1 & l=j+2
 sharp = sharpnrm*numer/(magg*denom)
 if not redo then begin
   redo = max(abs(x[j]) GT ( (0.05*chi*sqrt(c[j+nterm*j])) > 0.001*magg) )
   if redo EQ 0 then redo = max( abs([x[k], x[l]]) GT 0.01)
 endif

 sgn = where( xold[j]*x[j]/magg^2 LT -1.E-37, Nclamp )  
 if Nclamp GT 0 then clamp[j[sgn]] = 0.5*clamp[j[sgn]]
 sgn = where( xold[k]*x[k]        LT -1.E-37, Nclamp )
 if Nclamp GT 0 then clamp[k[sgn]] = 0.5*clamp[k[sgn]]
 sgn = where( xold[l]*x[l]        LT -1.E-37, Nclamp )
 if Nclamp GT 0 then clamp[l[sgn]] = 0.5*clamp[l[sgn]]

 magg = magg-x[j] / (1.+ ( (x[j]/(0.84*magg)) > (-x[j]/(5.25*magg)) )/ clamp[j] )
 xg = xg - x[k]   /(1.+abs(x[k])/( clamp[k]*0.5))
 yg = yg - x[l]   /(1.+abs(x[l])/( clamp[l]*0.5))
 xold = x

 magerr = c[j+nterm*j]*(nstr*chi^2 + (nstr-1)*chiold^2)/(2.*nstr-1.)

 dx = (-xg) > ( (xg - nx) > 0.) ;Find stars outside subarray
 dy = (-yg) > ( (yg-  ny) > 0.)
 badcen = where(    $                     ;Remove stars with bad centroids
     (dx GT 0.001) or (dy GT 0.001) or ( (dx+1)^2 + (dy+1)^2 GE radsq ), nbad)
 if nbad GT 0 then begin
        nstr = nstr - nbad
        print,strn(nbad),' stars eliminated by centroid criteria'
        if nstr LE 0 then goto, DONE_GROUP 
        remove, badcen, idg, xg, yg, magg, skyg, magerr
        nterm = nstr*3 + varsky
        redo = 1b
 endif

 faint = 1        
 toofaint =  where (magg LE 1.e-5,nfaint) 
                              ;Number of stars 12.5 mags fainter than PSF star
 if nfaint GT 0 then begin        
         faint = min( magg[toofaint], min_pos )
         ifaint = toofaint[ min_pos ]
         magg[toofaint] = 1.e-5
         goto, REM_FAINT                ;Remove faintest star
 endif else begin
         faint = 0.
         ifaint = -1
         if (not redo) or (niter GE 4) then $
            faint = max(magerr/magg^2, ifaint) else $ 
            goto,START_IT 
 endelse

 if keyword_set(DEBUG) then begin 
   err = 1.085736*sqrt(magerr)/magg
    for i=0,nstr-1 do  print,format=fmt,idg[i],xg[i]+ixmin,yg[i]+iymin, $
          psfmag-1.085736*alog(magg[i]),err[i],skyg[i],niter,chi[i],sharp[i]
 endif

 if redo and (niter LE 50) and (faint LT wcrit) then goto,START_IT   
REM_FAINT: 
 if (faint GE 0.25) or (nfaint GT 0) then begin
         if not SILENT then $
               message,'Star '+ strn(idg[ifaint]) + ' is too faint',/INF
         nstr = nstr-1
         if nstr LE 0 then goto,DONE_GROUP  
         remove,ifaint,idg,xg,yg,magg,skyg,magerr
         nterm = nstr*3 + varsky
         xold = dblarr(nterm)
         clamp = replicate(1.,nterm)
         clip = 0b
         niter = niter-1
         goto,RESTART  
 endif

 err = 1.085736*sqrt(magerr)/magg
 magg = psfmag - 1.085736*alog(magg)
 sharp = sharp > (-99.999) < 99.999
 xg = xg+ixmin & yg = yg+iymin

; Print results to terminal and/or file

 if not SILENT then for i = 0,nstr-1 do print,format=fmt, $
     idg[i],xg[i],yg[i],magg[i],err[i],skyg[i],niter,chi[i],sharp[i]
 if PRINT then for i = 0,nstr-1 do printf,lun,format=fmt, $
     idg[i],xg[i],yg[i],magg[i],err[i],skyg[i],niter,chi[i],sharp[i]

 if ( npar GE 9 ) then begin                  ;Create output vectors?
   if ( N_elements(newid) EQ 0 ) then begin  ;Initialize output vectors?
       newid = idg &  newx = xg  &  newy = yg & newmag = magg
       iter = replicate(niter,nstr) & peak = sharp & chisq = chi
       errmag = err
   endif else begin           ;Append current group to output vector
       newid = [newid,idg] & newx = [newx ,xg] & newy = [newy,yg]
       newmag = [newmag,magg] & iter = [iter,replicate(niter,nstr)]
       peak = [peak,sharp]     & chisq = [chisq,chi] & errmag = [errmag,err]
   endelse
 endif

DONE_GROUP: 
 endfor

 if  ( npar GE 9 ) then begin
    if N_elements(newid) GT 0 then begin
       id = newid &  xc = newx &  yc = newy  & mags = newmag
    endif else $
     message,'ERROR - There are no valid stars left, variables not updated',/CON
 endif

 if PRINT then free_lun,lun

 return
 end