nstar.pro

basic median filter simulation

pro nstar,image,id,xc,yc,mags,sky,group,phpadu,readns,psfname,DEBUG=debug, $  
          errmag,iter,chisq,peak,PRINT=print,SILENT=silent, VARSKY = varsky
;+
; NAME:
;       NSTAR
; PURPOSE:
;       Simultaneous point spread function fitting (adapted from DAOPHOT)
; EXPLANATION:
;       This PSF fitting algorithm is based on a very old (~1987) version of 
;       DAOPHOT, and much better algorithms (e.g. ALLSTAR) are now available
;       -- though not in IDL.
;       
; CALLING SEQUENCE:
;       NSTAR, image, id, xc, yc, mags, sky, group, [ phpadu, readns, psfname,
;               magerr, iter, chisq, peak, /PRINT , /SILENT, /VARSKY, /DEBUG ]
;
; INPUTS:
;       image - image array
;       id    - vector of stellar ID numbers given by FIND
;       xc    - vector containing X position centroids of stars (e.g. as found
;               by FIND)
;       yc    - vector of Y position centroids
;       mags  - vector of aperture magnitudes (e.g. as found by APER)
;               If 9 or more parameters are supplied then, upon output
;               ID,XC,YC, and MAGS will be modified to contain the new
;               values of these parameters as determined by NSTAR.
;               Note that the number of output stars may be less than 
;               the number of input stars since stars may converge, or 
;               "disappear" because they are too faint.
;       sky   - vector of sky background values (e.g. as found by APER)
;       group - vector containing group id's of stars as found by GROUP
;
; OPTIONAL INPUT:
;       phpadu - numeric scalar giving number of photons per digital unit.  
;               Needed for computing Poisson error statistics.   
;       readns - readout noise per pixel, numeric scalar.   If not supplied, 
;               NSTAR will try to read the values of READNS and PHPADU from
;               the PSF header.  If still not found, user will be prompted.
;       psfname - name of FITS image file containing the point spread
;               function residuals as determined by GETPSF, scalar string.  
;               If omitted, then NSTAR will prompt for this parameter.
;
; OPTIONAL OUTPUTS:
;       MAGERR - vector of errors in the magnitudes found by NSTAR
;       ITER - vector containing the number of iterations required for
;               each output star.  
;       CHISQ- vector containing the chi square of the PSF fit for each
;               output star.
;       PEAK - vector containing the difference of the mean residual of
;               the pixels in the outer half of the fitting circle and
;               the mean residual of pixels in the inner half of the
;               fitting circle
;
; OPTIONAL KEYWORD INPUTS:
;       /SILENT - if set and non-zero, then NSTAR will not display its results
;               at the terminal
;       /PRINT - if set and non-zero then NSTAR will also write its results to
;               a file nstar.prt.   One also can specify the output file name
;               by setting PRINT = 'filename'.
;       /VARSKY - if this keyword is set and non-zero, then the sky level of
;               each group is set as a free parameter.
;       /DEBUG - if this keyword is set and non-zero, then the result of each
;               fitting iteration will be displayed.
;
; PROCEDURES USED:
;       DAO_VALUE(), READFITS(), REMOVE, SPEC_DIR(), STRN(), SXPAR()
;
; COMMON BLOCK:
;       RINTER - contains pre-tabulated values for cubic interpolation
; REVISION HISTORY
;       W. Landsman                 ST Systems Co.       May, 1988
;       Adapted for IDL Version 2, J. Isensee, September, 1990
;       Minor fixes so that PRINT='filename' really prints to 'filename', and
;       it really silent if SILENT is set.  J.Wm.Parker HSTX 1995-Oct-31
;       Added /VARSKY option   W. Landsman   HSTX      May 1996
;       Converted to IDL V5.0   W. Landsman   September 1997
;       Replace DATATYPE() with size(/TNAME)  W. Landsman November 2001
;       Assume since V5.5, remove VMS calls W. Landsman September 2006
;-
 compile_opt idl2
 common rinter,c1,c2,c3,init            ;Save time in RINTER()
 npar = N_params()
 if npar LT 7 then begin
   print,'Syntax - NSTAR, image, id, xc, yc, mags, sky, group, [phpadu, '
   print, $
     '   [readns, psfname, magerr, iter, chisq, peak, /SILENT, /PRINT, /VARSKY]'
   return
 endif                                                        

 if ( N_elements(psfname) EQ 0 ) then begin
   psfname=''
   read,'Enter name of FITS file containing PSF: ',psfname
 endif else zparcheck,'PSFNAME',psfname,10,7,0,'PSF disk file name'

 psf_file = file_search( psfname, COUNT = n)
 if n EQ 0 then message, $
      'ERROR - Unable to locate PSF file ' + spec_dir(psfname)

 if npar LT 9 then begin                                                     
   ans = ''
   read, $
     'Do you want to update the input vectors with the results of NSTAR? ',ans
   if strmid(strupcase(ans),0,1) EQ 'Y' then npar = 9
 endif

 if npar LT 9 then $
    message,'Input vectors ID,XC,YC and MAGS will not be updated by NSTAR',/INF

; Read in the FITS file containing the PSF

 s = size(image)
 icol = s[1]-1 & irow = s[2]-1  ;Index of last row and column
 psf = readfits(psfname, hpsf)
 if  N_elements(phpadu) EQ 0 then begin 
    par = sxpar(hpsf,'PHPADU', Count = N_phpadu)
    if N_phpadu eq 0 $
        then read, 'Enter photons per analog digital unit: ',phpadu $
        else phpadu = par
endif

 if ( N_elements(readns) EQ 0 ) then begin 
    par = sxpar(hpsf,'RONOIS', Count = N_ronois)
    if N_ronois EQ 0 $
        then read, 'Enter the readout noise per pixel: ',readns $
        else readns = par
 endif

 gauss = sxpar(hpsf,'GAUSS*')
 psfmag = sxpar(hpsf,'PSFMAG')
 psfrad = sxpar(hpsf,'PSFRAD')
 fitrad = sxpar(hpsf,'FITRAD')
 npsf = sxpar(hpsf,'NAXIS1')
;                               Compute RINTER common block arrays
 p_1 = shift(psf,1,0) & p1 = shift(psf,-1,0) & p2 = shift(psf,-2,0)
 c1 = 0.5*(p1 - p_1)
 c2 = 2.*p1 + p_1 - 0.5*(5.*psf + p2)
 c3 = 0.5*(3.*(psf-p1) + p2 - p_1)
 init = 1

 ronois = readns^2
 radsq = fitrad^2   &  psfrsq = psfrad^2
 sepmin = 2.773*(gauss[3]^2+gauss[4]^2)

;      PKERR will be used to estimate the error due to interpolating PSF
;      Factor of 0.027 is estimated from good-seeing CTIO frames

 pkerr = 0.027/(gauss[3]*gauss[4])^2     
 sharpnrm = 2.*gauss[3]*gauss[4]/gauss[0]
 if (N_elements(group) EQ 1) then groupid = group[0] else $
     groupid = where(histogram(group,min=0))    ;Vector of distinct group id's

 mag = mags                        ;Save original magnitude vector
 bad = where( mag GT 99, nbad )     ;Undefined magnitudes assigned 99.9
 if nbad GT 0 then mag[bad] = psfmag + 7.5
 mag = 10.^(-0.4*(mag-psfmag)) ;Convert magnitude to brightness, scaled to PSF
 fmt = '(I6,2F9.2,3F9.3,I4,F9.2,F9.3)'

 SILENT = keyword_set(SILENT)
 VARSKY = keyword_set(VARSKY)

 if keyword_set(PRINT) then begin
     if ( size(print,/TNAME) NE 'STRING' ) then file = 'nstar.prt' $
                                           else file = print
    message,'Results will be written to a file '+ file,/INF
     openw,lun,file,/GET_LUN
     printf,lun,'NSTAR:    '+ getenv('USER') + ' '+ systime()
     printf,lun,'PSF File:',psfname
 endif 
 PRINT = keyword_set(PRINT)

 hdr='   ID      X       Y       MAG     MAGERR   SKY   NITER     CHI     SHARP'
 if not(SILENT) then print,hdr
 if PRINT then printf,lun,hdr

 for igroup = 0, N_elements(groupid)-1 do begin

 index = where(group EQ groupid[igroup],nstr) 
 if not SILENT then print,'Processing group ', $
               strtrim(groupid[igroup],2),'    ',strtrim(nstr,2),' stars'
 if nstr EQ 0 then stop
 magerr = fltarr(nstr)
 chiold = 1.0
 niter = 0
 clip = 0b
 nterm = nstr*3 + varsky
 xold = dblarr(nterm)
 clamp = replicate(1.,nterm)
 xb = double(xc[index])  &   yb = double(yc[index])
 magg = double(mag[index]) & skyg = double(sky[index])
 idg = id[index]
 skybar = total(skyg)/nstr
 reset = 0b
;
START_IT : 
   niter = niter+1
RESTART: 
 case 1 of              ;Set up critical error for star rejection
   niter GE 4 : wcrit = 1
   niter GE 8 : wcrit = 0.4444444
   niter GE 12: wcrit = 0.25             
   else       : wcrit = 400                                                   
 endcase

 if reset EQ 1b then begin
     xb = xg + ixmin & yb = yg + iymin
 endif

 reset = 1b
 xfitmin = fix(xb - fitrad) > 0
 xfitmax = fix(xb + fitrad)+1 < (icol-1)
 yfitmin = fix(yb - fitrad) > 0
 yfitmax = fix(yb + fitrad)+1 < (irow-1)
 nfitx = xfitmax - xfitmin + 1
 nfity = yfitmax - yfitmin + 1
 ixmin = min(xfitmin)& iymin = min(yfitmin)
 ixmax = max(xfitmax)& iymax = max(yfitmax)
 nx = ixmax-ixmin+1 & ny = iymax-iymin+1
 dimage = image[ixmin:ixmax,iymin:iymax]
 xfitmin = xfitmin -ixmin & yfitmin = yfitmin-iymin
 xfitmax = xfitmax -ixmin & yfitmax = yfitmax-iymin
;                                        Offset to the subarray
 xg = xb-ixmin & yg = yb-iymin
 j = 0

 while (j LT nstr-1) do begin
   sep = (xg[j] - xg[j+1:*])^2 + (yg[j] - yg[j+1:*])^2
   bad = where(sep LT sepmin,nbad)
   if nbad GT 0 then begin      ;Do any star overlap?
      for l = 0,nbad-1 do begin
      k = bad[l] + j + 1
      if magg[k] LT magg[j] then imin = k else imin = j ;Identify fainter star
      if ( sep[l] LT 0.14*sepmin) or  $
         ( magerr[imin]/magg[imin]^2 GT wcrit ) then begin
      if  imin EQ j then imerge = k else imerge = j
      nstr = nstr - 1
      if not SILENT then print, $
       'Star ',strn(idg[imin]),' has merged with star ',strn(idg[imerge])
      totmag = magg[imerge] + magg[imin]
      xg[imerge] = (xg[imerge]*magg[imerge] + xg[imin]*magg[imin])/totmag
      yg[imerge] = (yg[imerge]*magg[imerge] + yg[imin]*magg[imin])/totmag
      magg[imerge] = totmag     
      remove,imin,idg,xg,yg,magg,skyg,magerr    ;Remove fainter star from group
      nterm = nstr*3 + varsky                   ;Update matrix size