pro dao, img0, imgexp, fwhm, gain, outname, $ px, py, pmags, perrmags, $ apx, apy, apmags, aperrmags, $ imgsub, $ model = model, $ zpt = zpt, $ sigcut = sigcut, $ medfilt = medfilt0, $ psfrad = psfrad0, fitrad = fitrad0, $ aprad = aprad0, skyrad = skyrad0 ;+ ; convenient wrapper for running DAOPHOT routines on an image ; to do aperture & psf fitting photometry ; ; script will ; 1. median-filter the image to remove large-scale structure ; 2. identify peaks ; 3. do aperture phot, ; 4. use criteria (S/N, shape, FWHM, and isolation) to select psf stars ; 5. extract the psf stars from the image ; 6. do psf-fitting phot ; 7. write aperture & psf-fitting phot results to a file ; 8. check for blends (disabled for now) ; ; INPUTS ; img input image ; imgexp exposue map for input image ; fwhm FWHM - used for finding psf so should be fairly accurate ; gain instrument gain (for imgexp=1) ; needs to be a reasonable value ; if (much) too large or small, will oversubtract ; also affects error estimates (which i think are bogus anyways) ; outname string prefix for output files - see below ; ; OUTPUTS ; px,py,pmags,perrmags coords & mags from psf-fitting phot ; apx,apy,apmags,aperrmags coords & mags from aperture phot ; imgsub image after point-source subtraction ; ; also writes a bunch of output files with the prefix 'outname': ; *.hhh, *.hhd DAOPHOT psf information ; *.aper aperture phot results ; *.nstar psf fitting results ; ; KEYWORD INPUTS ; model model image to be subtracted (i.e. elliptical galaxy), ; will not do median filtering ; ; KEYWORD PARMS ; zpt effective zeropoint for image ; sigcut sigma threshold for object detection ; medfilt box size for median filtering (units of FWHM) ; psfrad extraction radius for psf (units of FWHM) ; fitrad radius for psf fitting (units of FWHM) ; aprad radii for aperture phot (units of FWHM) ; skyrad inner and outer radii for ; sky annulus for aperture phot (units of FWHM) ; ; unresolved issues: ; - what is proper S/N value of the peaks from FIND.PRO? ; need to know if order to determine psf-selection cutoff ; though value I am using now works ok ; - specified S/N level ('sigcut') is not the actual one ; to do right, need to scale by the 'relative error' ; printed by FIND.PRO (10-20% effect) which accounts for ; the weighting of the detection scheme -> not a big deal ; - I suspect all the computed 'errors' are nonsense... ; - what does the printed output of NSTAR mean? ; many objects have the same ID # and almost the same position ... ; - need to apply *aperture correction* to aperture & psf mags ; - psf extraction & aperture phot is done on *filtered* image, ; not the unsharp masked one (I think this is right) ; ; started as DAO.IDL ; ; 05/13/99 MCL ; ; can pass elliptical galaxy model through 'model' keyword ; if this exists, it will be subtracted from the original image ; before any photometric processing, w/no unsharp masking ; if no model is passed, then unsharp masking is still done ; 05/25/99 MCL ;- ; default zeropoint and approx S/N threshold for detection ; to do right, use 'relative error' from FIND.PRO (10-20% effect) if not(keyword_set(zpt)) then zpt = 25.0 if not(keyword_set(sigcut)) then sigcut = 4.0 ; default sizes (in units of FWHM) if not(keyword_set(medfilt0)) then medfilt0 = 5.0 if not(keyword_set(psfrad0)) then psfrad0 = 4.0 if not(keyword_set(fitrad0)) then fitrad0 = 1.0 if not(keyword_set(aprad0)) then aprad0 = 2.5 if not(keyword_set(skyrad0)) then skyrad0 = [4, 6] ; user info if n_params() lt 5 then begin print, 'pro dao, img, imgexp, fwhm, gain, outname, {in}' print, ' px, py, pmags, perrmags, apx, apy, apmags, aperrmags, {out}' print, ' [sigcut='+strc(sigcut)+'],' print, ' [medfilt='+strc(medfilt0)+'], [psfrad='+strc(psfrad0)+ $ '], [fitrad='+strc(fitrad0)+'], [aprad='+strc(aprad0)+ $ '], [skyrad='+commalist(skyrad0)+']' print, ' note: all parms on last line are in units of the FWHM, not pixels' retall endif ;------------------------------------------------------------ ; initialization & hidden parms ;------------------------------------------------------------ BADVAL = -1e6 win, sz = 350, 0 win, sz = 350, 1 win, sz = 350, 2 psfname = outname ; convert from units of FWHM to pixel units medfilt = medfilt0 * fwhm psfrad = psfrad0 * fwhm fitrad = fitrad0 * fwhm aprad = aprad0 * fwhm skyrad = skyrad0 * fwhm ; psf star selection criteria: ; DWFHM: must be within this amt of user-specified FWHM ; PSFS2N: psf S/N threshold (faked S/N calculation) ; PSFROUND: max magnitude for psf 'roundness' PSF_DFWHM = 0.5 PSF_S2N = 20 PSF_ROUND = 0.2 ; inner sky radius must be larger than psf extraction radius if (skyrad(0) lt psfrad) then begin message, 'inner sky radius < psf radius == bad!', /info stop endif ; for Goddard routines defsysv,'!TEXTOUT',1 defsysv,'!TEXTUNIT',0 defsysv,'!DEBUG',0 ; DAOPHOT parms - don't change roundlim = [-1.0, 1.0] ; roundness stat cutoff - change for elongated stars sharplim = [0.2, 1.0] ; sharpness stat - change only if the stars have ; much larger or or smaller concentration ; than a Gaussian ; sanity check on images sz = size(img0) szexp = size(imgexp) if (sz(0) ne 2) or (szexp(0) ne 2) then begin message, 'input image is not 2d!' retall endif if total(sz ne szexp) ne 0 then begin message, 'image & exposure map are not the same size!' retall endif ; check if psf & phot files already exist repeat begin done = 1 if exist(psfname+'.hhh') then begin print, 'DAOPHOT files with the prefix of "', psfname, '" already exist!' if getyn('do you wish to change prefix name?', 0) then begin read, 'enter new prefix: ', psfname done = 0 endif else $ print, 'ok - overwriting existing files ...' endif endrep until (done) ;-------------------------------------------------- ; prepare the input image and noise-normalize ;-------------------------------------------------- ; subtract model image if passed ; otherwise do median filtering if keyword_set(model) then begin szmod = size(model) if (total(sz-szmod) ne 0) then begin message, 'input and model images are diff sizes!', /info stop endif message, 'subtracting model image passed by user', /info img = imsub(img0, model) endif else begin message, 'median filtering image, box size = '+ $ strc(medfilt)+' pixels ('+strc(medfilt0)+' x FWHM)', /info imgmed = filter_image(img, med = medfilt, /all) img = imsub(img0, imgmed) endelse ; create bad pixel fixed image fixpix, img, img ne BADVAL, imgfix, /silent ; create noise-normalized image ; this means the effective gain is non-uniform over the image ; as is the limiting magnitude/completeness limits imgnorm = img * sqrt(imgexp) ; use iterstat sigma as readnoise in DN, convert to e- iterstat, imgnorm, istat, /silent rn = istat(3)*gain ;------------------------------------------------------------ ; run DAOPHOT on the object image ;------------------------------------------------------------ ; determine sky level in normalized image sky, imgnorm, nskymode, nskysig, /silent uscl = 3.0*nskysig ; find positive perturbations (i.e stars) in noise-normalized image ; and compute centroids and shape parameters, ; use sky noise to estimate sigma in the image, rather ; than computing directly from gain & RN calc - ; this should be more accurate since images are sky subtracted hmin = sigcut*nskysig ; threshold intensity for a point source find, imgnorm, x, y, flux, sharp, roundness, $ hmin, fwhm, roundlim, sharplim, /silent nfind = n_elements(x) print, 'Final number of sources = ', strc(nfind) ; display found objects wset, 0 display2, /silent, imgnorm < uscl, tit = 'noise-normalized image' wset, 1 display2, /silent, imgnorm < uscl, $ tit = 'found objects!Ccircle = phot aperture' ;oplot, x, y, ps = 1, sym = 5, col = 0 for i = 0, nfind-1 do $ tvcircle, /data, aprad, x(i), y(i), col = 0 ; Do photometry with FWHM-sized aperture on noise-normalized image ; to determine approx S/N of each detected object ; this assumes uniform weighting, which FIND.PRO doesn't use - ; to do right, use 'relative error' printed by FIND.PRO: 10-20% effect) aper, imgnorm, x, y, s2nraw, /flux, $ errap, sky, skyerr, $ gain, fwhm/2.0, skyrad, [(min(imgnorm)-1) < BADVAL, max(imgnorm+1)], /SILENT s2n = s2nraw/(sqrt(!pi)*fwhm/2.0*nskysig) ; do aperture phot on original image to get real fluxes ; a separate sky value is determined for each source aper, imgfix, x, y, mags, errmags, sky, skyerr, $ gain, aprad, skyrad, [(min(imgfix)-1) < BADVAL, max(imgfix+1)], $ print = psfname+'.aper', /SILENT ; print = psfname+'.aper', /SILENT, /nosky ;print, '/nosky is set!!!' ; group the nearby stars rmax = psfrad+fitrad group, x, y, rmax, group ;-------------------------------------------------- ; extract psf from image ;-------------------------------------------------- ; use my RADPLOTF routine to get FWHM's for objects objfwhm = fltarr(nfind) for i = 0, nfind-1 do begin radplotf, imgfix, x(i), y(i), radout, ff, $ outrad = skyrad(0), insky = skyrad(0), outsky = skyrad(1), $ /silent objfwhm(i) = ff endfor ; choose candidate stellar objects as those with ; FWHM within +/- 0.5 pixels of the user-specified FWHM wpsf0 = where(abs(objfwhm-fwhm) le PSF_DFWHM, npsf0) if (npsf0 eq 0) then begin message, 'no objects pass FWHM criteria for psf selection!', /info stop endif message, 'number of PSF candidates = '+strc(npsf0), /info ; from the candidate psf stars, choose only those which: ; - are isolated, as designated by GROUP routine ; - have abs(ROUNDNESS) <= PSF_ROUND ; - have S/N>=20 for aperture phot npsf = 0 for i = 0, npsf0-1 do begin ii = wpsf0(i) g0 = group(ii) ; group # of candidate psf w = where(group eq g0, nw) if (nw eq 1) and (abs(roundness(ii)) le PSF_ROUND) $ and (s2n(ii) ge PSF_S2N) then begin print, 'psf '+strc(npsf), x(ii), y(ii), $ objfwhm(ii), s2n(ii), roundness(ii) if (npsf eq 0) then wpsf = ii $ else wpsf = [wpsf, ii] npsf = npsf+1 endif else $ print, ' bad', x(ii), y(ii), objfwhm(ii), s2n(ii), roundness(ii) endfor message, 'number of PSF stars = '+strc(npsf), /info if (npsf eq 0) then begin message, 'no psf stars!!', /info stop endif ; identify psf stars wset, 2 display2, /silent, imgnorm < uscl, tit = strc(npsf)+' psf stars!Ccircle = psf size' for i = 0, npsf-1 do $ tvcircle, /data, psfrad, x(wpsf(i)), y(wpsf(i)), col = 0, thick = 2 ; extract psf from the ORIGINAL image ; in case portions of psf star have diff exposure times sky, imgfix, skymode, skysig, /silent idpsf = wpsf getpsf, imgfix, x, y, mags, sky, rn, $ gain, gauss, psf, idpsf, psfrad, fitrad, psfname win, sz = 350, 3 display2, /silent, psf, tit = 'psf: non-gaussian residuals' ; write list of psf stars to file openw, unit0, psfname+'.psf', /get_lun printf, unit0, '# DAO.PRO: ', systime(), userid() printf, unit0, '# list of psf stars selected' printf, unit0, '#' printf, unit0, '# input parms:' printf, unit0, '# fwhm = ', strc(fwhm) printf, unit0, '# gain = ', strc(gain) printf, unit0, '# aprad = ', strc(aprad) printf, unit0, '# skyrad = ', commalist(skyrad) printf, unit0, '#' printf, unit0, '# psf selection criteria:' printf, unit0, '# PSF_DFWHM = ', strc(PSF_DFWHM) printf, unit0, '# PSF_S2N = ', strc(PSF_S2N) printf, unit0, '# PSF_ROUND =', strc(PSF_ROUND) printf, unit0, '#' printf, unit0, '# ID X Y MAG SKY' free_lun, unit0 printarray, [[idpsf], [x(idpsf)], [y(idpsf)], [mags(idpsf)], [sky(idpsf)]], $ outfile = psfname+'.psf', /silent ;-------------------------------------------------- ; do psf-fitting photometry ;-------------------------------------------------- ; first save results from aperture phot apx = x apy = y apmags = mags aperrmag = errmags ; run NSTAR on all the detected objects ; be sure to save aperture phot mags & positions ; note that NSTAR may return a different # of stars than inputted ; b/c some may be too faint (<0 net flux) or rejected by the fit id = indgen((size(x))(1)) nstar, imgfix, id, x, y, mags, sky, group, gain, rn, psfname, $ errmags, iter, chisq, peak, print = psfname+'.nstar', /silent ; rename the psf variable for clarity px = x py = y pmags = mags perrmags = errmags ; show which stars had good fit win, sz = 350, 4 display2, imgnorm < uscl, tit = 'stars removed by PSF fitting' w = where(pmags lt 99, nw) for i = 0, nw-1 do $ tvcircle, /data, aprad, px(i), py(i), col = 0 ; subtract stars imgsub = imgfix substar, imgsub, x, y, pmags, -1, psfname ; display subtracted image win, sz = 350, 5 display2, /silent, imgsub < uscl/3.0, tit = 'subtracted image' ; correct mags (aperture & psf-fitted) to correct zpt ; don't do this until the very end since mags used by ; the other routines (GETPSF, NSTAR, SUBSTAR) w = where(pmags ne 99.999) pmags(w) = pmags(w)-25+zpt w = where(apmags ne 99.999) apmags(w) = apmags(w)-25+zpt ;;------------------------------------------------------------ ;; do 2nd pass FIND to see if there are any companions ;; use same threshold (though maybe should increase it a bit) ;; -> CURRENTLY DISABLED ;;------------------------------------------------------------ ;find, imgsub, sx, sy, sflux, ssharp, sroundness,$ ; hmin, fwhm, roundlim, sharplim, /silent ;ws = where(sflux gt 0, ns) ; ;; if there are companions, append them to object list ;; and re-do APER & NSTAR to get phot of complete list ;if (ns gt 0) then begin ; ; print, '** number of companions = ', strc(ns) ; ; ; overplot on residual image ; oplot, sx(ws), sy(ws), ps = 1, sym = 13, col = 0, thick = 2 ; for i = 0, ns-1 do $ ; tvcircle, /data, 2, sx(ws(i)), sy(ws(i)), col = 0 ; ; ; append new companions to old list ; x = [x, sx(ws)] ; y = [y, sy(ws)] ; nfind = nfind + n_elements(ns) ; ; ; re-do APER ; aprad = fwhm ; aper, img, x, y, mags, errap, sky, skyerr, $ ; gain, aprad, skyrad, [min(img)-1, max(img)+1], /SILENT ; ; ; group the nearby stars (trivial step for this work) ; group, x, y, rmax, group ; ; ; do all the stars from find ; id = indgen(nfind) ; nstar,img,id,x,y,mags,sky,group,gain,rn,psfname, $ ; errmag,iter,chisq,peak ; ; ; subtract stars ; imgsub2 = img ; substar,imgsub2,x,y,mags,-1,psfname ; ; ; correct mags (aperture & psf-fitted) to correct zpt ; w = where(mags ne 99.999) ; mags(w) = mags(w)-25+zpt ; w = where(apmags ne 99.999) ; apmags(w) = apmags(w)-25+zpt ; ; ; display center of subtracted image ; win, sz = 350, 3, tit = 'IDL 3: 2nd-pass residuals' ; display2, /silent, imgsub2, /tv ; ;endif else $ ; print, '** no companions found ' ;-------------------------------------------------- ; Postscript plots: (1) number counts ;-------------------------------------------------- ps_open, psfname+'.daophot', thick = 4, /port !p.multi = [0, 1, 2] ; get good mags and set range wg_ap = where(apmags lt 99, ng_ap) wg_psf = where(pmags lt 99, ng_psf) datap = apmags(wg_ap) datp = pmags(wg_psf) dat = [datap, datp] xr = [floor(min(dat)), ceil(max(dat))] ; construct histograms xhist = makex(xr(0), xr(1), 0.5) aphy = histogram(datap, bin = 0.5, min = xr(0)-0.25, max = xr(1)+0.25) phy = histogram(pmags, bin = 0.5, min = xr(0)-0.25, max = xr(1)+0.25) yr = [0, max([aphy, phy])+0.5] ; plot histogram for aperture mags ; (slightly shift histogram for PSF phot for greater clarity) plot, xhist, aphy, ps = 10, $ xr = xr, xs = 1, yr = yr, $ tit = psfname+': number counts', $ xtit = 'mag', ytit = 'number', $ chars = 1.5 oplot, xhist-0.05, phy+0.05, line = 1, ps = 10 ; label number stats xyouts, /data, !x.crange(0)+0.5, !y.crange(1)*0.9, chars = 1, $ textoidl('FIND threshold = '+strc(sigcut)+' \sigma'+ $ '!C!CAPER PHOT (solid):'+ $ '!C N_{tot}='+strc(n_elements(apx))+', N_{good}='+strc(ng_ap)+ $ '!C aprad ='+string(aprad, '(F4.1)')+' pix'+ $ '!C!CPSF PHOT (dotted):'+ $ '!C N_{tot}='+strc(n_elements(px))+', N_{good}='+strc(ng_psf)+ $ '!C # of psf stars = '+strc(npsf)) ;-------------------------------------------------- ; Postscript plot: (2) aperture vs psf-fitted mags ;-------------------------------------------------- ; loop through psf-fitted stars, and find ; the corresponding star from the aperture phot list it originated from ; (this may screw up in crowded fields) n = n_elements(px) dist = fltarr(n) dmag = fltarr(n) for i = 0, n-1 do begin mindist = min((px(i)-apx)^2.0 + (py(i)-apy)^2.0, wmin) dmag(i) = pmags(i)-apmags(wmin) print, printcoo(px(i), py(i)), printcoo(apx(wmin), apy(wmin)), $ mags(i), apmags(wmin) endfor ; plot difference plot, pmags(wg_psf), dmag(wg_psf), ps = 8, sym = 1.5, $ xtit = 'PSF mag', ytit = '(PSF - APERTURE) mag', $ xr = xr, xs = 1, $ tit = psfname+': PSF versus APERTURE mags', $ chars = 1.5 plots, !x.crange, [0, 0], line = 1 datestamp, 'DAO.PRO' ps_close gv, psfname+'.daophot.ps' stop end