pro make_sigma_image_new,image_fits,mask_fits,wt_fits,sigma_fits,$ check=check,fudge0=fudge0,$ test=test,band0=band0,stop0=stop0,$ image_new_fits=image_new_fits,$ mask_new_fits=mask_new_fits,datadir=datadir,$ sky_subtract0=sky_subtract0 if(n_params() le 0 and not keyword_set(test)) then begin print,'-------------------------------------------------------------------' print,'make_sigma_image, image_fits, mask_fits, wt_fits, sigma_fits' print,'-------------------------------------------------------------------' print,'Make a sigma-image based on data, mask, and weight images' print,'-------------------------------------------------------------------' print,'Input:' print,' image_fits = name of the fits-file of data image (unit=Mjy/sr)' print,' mask_fits = -- " -- mask image (0-good pixel)' print,' wt_fits = -- " -- weight im (#frames/pixel*10)' print,'' print,'Output:' print,' sigma_fits = output sigma image is written (unit=MJy/sr)' print,'Output via keywords:' print,' image_new_fits = image cleaned of NaNs, also set exptime=1' print,' mask_new_fits = masks which exclude NaNs of image and sigma-image' print,'' print,' (output is also stored to IDL sav-file sigma_fits.save)' print,'--------------------------------------------------------------------' print,'Keywords:' print,' check=[xc,yc,wpix] ->' print,' compare stdev in the image with the calculated sigma-image' print,' using circular annulae of width wpix centered at xc,yc' print,' (note xc,yc are in idl-coordinates starting from 0,0)' print,'' print,' fudge = fudge-factor by which the calculated sigmas are multiplied' print,' FOR SOME UNKNOWN REASON fudge=2 SEEMS TO WORK BEST' print,'' print,'---------------------------------------------------------------------' print,'Examples:' print," make_sigma_image_new,'ngc4594.newcut.fits','ngc4594.newmask.fits',$" print," 'ngc4594.newwt.fits','ngc4594.newsigma.fits',$" print," fudg=2,check=[517.579,451.251]" print,'' print,'Normally called by make_sigma_image_driver.pro' print,'---------------------------------------------------------------------' print,'Method: using header information, convert image numbers to electrons,' print,'calculate poisson noise + readout-noise in electrons' print,'then convert sigma(electrons) back to image numbers (MJy/sr)' print,'---------------------------------------------------------------------' print,' heikki.salo@oulu.fi 280809/050110' print,'---------------------------------------------------------------------' print,' ADDITION 090610 HS:' print,' since new P1 images are not background-subtracted, ' print,' this can be done here, after calculation of sigma-images' print,' keywosrd sky_subtract=value --> ' print,' subtract value from the nonan-image stored by this procedure' return endif ;==================================================== ; factor to multiply the estimated sigma ; should be 1, but 2 seems to make more sense! ;==================================================== fudge=1. if(keyword_set(fudge0)) then fudge=fudge0 ;==================================================== ;read ngc1097 test-data ? ;==================================================== if(keyword_set(test)) then begin if(test eq 1) then begin ide='ngc1097_36_sigma' ima=readfits('ngc1097_v7.phot.1.fits',h1) wt=readfits('ngc1097_v7.phot.1_wt.fits') mask=readfits('mask_ngc1097.fits',h2) XC = 750.881 YC = 738.824 wpix=10. band=1 sigma_fits='ngc1097_36_sigma.fits' endif if(test eq 2) then begin ide='ngc1097_45_sigma' ima=readfits('ngc1097_v7.phot.2.fits',h1) wt=readfits('ngc1097_v7.phot.2_wt.fits') mask=readfits('ngc1097_mask.2.fits',h2) XC = 1001.98 YC = 1219.95 wpix=10. band=2 sigma_fits='ngc1097_45_sigma.fits' endif if(test ne 1 and test ne 2) then begin print,'use test=1 orr test=2' return endif endif ;==================================================== ; use input data ;==================================================== if(not keyword_set(test)) then begin usedir='' if(keyword_set(datadir)) then usedir=datadir ide=image_fits ima=readfits(usedir+image_fits,h1) wt=readfits(usedir+wt_fits,/silent) mask=readfits(usedir+mask_fits,h2,/silent) size_ima=size(ima) size_wt=size(wt) size_mask=size(mask) if(size_ima(4) ne size_wt(4)) then begin print,'sizes of data and weight images do not match -> check input' return endif if(size_ima(4) ne size_mask(4)) then begin print,'sizes of data and mask images do not match -> check input' return endif band=1 if(keyword_set(band0)) then band=band0 endif ;============================================================== ;read header data ;============================================================== ;example header (ngc1097_v7.phot.1.fits) : ;============================================================== ;FRAMTIME= 30. ;/ [sec] Time spent integrating (whole array) ;;COMMENT Photons in Well = Flux[photons/sec/pixel] * FRAMTIME ;EXPTIME = 1876. ;/ [sec] Effective integration time per pixel ;;COMMENT DN per pixel = Flux[photons/sec/pixel] / GAIN * EXPTIME ;BUNIT = 'MJy/sr ' ;/ Units of image data ;FLUXCONV= 0.1088 ;/ Flux Conv. factor (MJy/sr per DN/sec) ;GAIN = 3.3 ;/ e/DN conversion ;RONOISE = 7.8 ;/ [Electrons] Readout Noise from Array ;ZODY_EST= .06041035 ;/ [MJy/sr] Zodiacal Background Estimate ;BACKGRND= 0.0367242 ;============================================================== gain=sxpar(h1,'GAIN') ronoise=sxpar(h1,'RONOISE') fluxconv=sxpar(h1,'FLUXCONV') framtime=sxpar(h1,'FRAMTIME') backgrnd=sxpar(h1,'BACKGRND') ;============================================================== ;Find bad pixels in the image (=NaN) and set the mask to -1 for these pixels ;indbad contains list of pixels not to use in smoothing lnan=where(finite(ima) eq 0, cnan) if(cnan gt 0) then mask(lnan)=-1. indbad_ima=lnan cnan_ima=cnan lnan=where(finite(wt) eq 0, cnan) if(cnan gt 0) then mask(lnan)=-1. indbad=where(mask ne 0, cbad) indok=where(mask eq 0) ;number of frames/pixel weight=wt/10. wt=0. ;to save memory ;============================================================================= ;convert Mjy/sr -> orig. digital units/second -> digital units -> electrons ele=(ima+backgrnd)/fluxconv*framtime*weight*gain>1d-12 ;sigma of electrons: add Poisson noise and RON quadratically sigma_ele=sqrt(ele+weight*ronoise^2) ;convert to sigma of dnsec and sigma of ima sigma_dnsec=sigma_ele/(gain*framtime*weight) sigma_ima=sigma_dnsec*fluxconv * fudge ;NOTE FUDGE!! sigma_dnsec=0. ;to save memory ;sigma without RON (just for test) sigma_ele0=sqrt(ele) sigma_dnsec0=sigma_ele0/(gain*framtime*weight) sigma_ima0=sigma_dnsec0*fluxconv sigma_dnsec0=0. ;to save memory ;--------------------------------------------------------------------- ;make sure that sigma is not unnecessarily stored in double precisoin sigma_ima=float(sigma_ima) ;smooth sigma_ima by 20 pixels (arbitrary, refine later) ;before smoothing, replace bad values by median medi=median(sigma_ima(indok)) sigma_ima=smooth(sigma_ima,20,/nan,missing=medi) ;--------------------------------------------------------------------- ;write sigma-image to fits file ;mark bad values with -1 (they will be excluded by the mask_new) if(cnan_ima ge 1) then begin ima(indbad_ima)=-1. ;eliminate NaNs from the image sigma_ima(indbad_ima)=-1. endif h11=h1 sxaddhist,'sigma-image made by make_sigma_image.pro '+!stime,h1,/comment writefits,sigma_fits,sigma_ima,h1 ;write new mask? removes Nan in data and weight if(keyword_set(mask_new_fits)) then begin sxaddhist,'mask-image modified by make_sigma_image.pro '+$ !stime,h2,/comment writefits,mask_new_fits,mask,h2 endif ;write new data image? Does not have NaNs ;for galfit use texp=1sec sky_subtract=0. if(keyword_set(image_new_fits)) then begin sxaddhist,'data-image modified by make_sigma_image.pro '+$ !stime,h11,/comment sxaddhist,'NaNs removed',h11,/comment sxaddhist,'exptime set to 1 second',h11,/comment sxaddpar,h11,'exptime','1.' if(keyword_set(sky_subtract0)) then begin sky_subtract=sky_subtract0 ima=ima-sky_subtract sxaddhist,'subtracted sky value '+string(sky_subtract),h11,/comment endif writefits,image_new_fits,ima,h11 endif save,file=sigma_fits+'.save',ima,sigma_ima,mask,fudge,sky_subtract ;---------------------------------------------------------------------- ;compare actual and calculated sigma in circular annulae ;---------------------------------------------------------------------- if(keyword_set(test) or keyword_set(check)) then begin if(keyword_set(check)) then begin if(n_elements(check) ne 3) then begin print,'must specify check=[xc,yc,wpix]' return endif xc=check(0) yc=check(1) wpix=check(2) endif ;create circular distances via ASTRO ROUTINE nx=n_elements(ima(*,0)) ny=n_elements(ima(0,*)) dist_ellipse,rad,[nx,ny],xc,yc,1.,0. ;makes a window (works also with PS) nwin_hs ;calculate actual image stdev (=dyt) in circular zones, excluding bad values if(cbad ge 1) then ima(indbad)=-999. tabulate_f,rad,ima,0,1500.,wpix,xt,yt,dyt=dyt,myt=myt,bad=-999. plot,xt,dyt,psym=0,/ylog,yr=[0.0001,100],xtitle='pixels from center',$ ytitle='sigma(image) [MJy/sr]',title=ide,/nod,pos=[0.12,0.1,0.9,0.9] oplot,xt,dyt,col=2 ;tabulate theoretical mean sigma (=myt) in a similar fahion if(cbad ge 1) then sigma_ima(indbad)=-999. tabulate_f,rad,sigma_ima,0,1500.,wpix,xt,yt2,dyt=dyt2,myt=myt2,/finite,$ bad=-999. oplot,xt,myt2,col=3,lines=2 if(cbad ge 1) then sigma_ima0(indbad)=-999. tabulate_f,rad,sigma_ima0,0,1500.,wpix,xt,yt2,dyt=dyt2,myt=myt2,/finite,$ bad=-999. oplot,xt,myt2,col=5,lines=1 ff='(i2)' label_data,0.5,0.9,['STDEV(image)','MEAN(sigma_image) * '+$ string(fudge,ff),'without RON'],$ title='IN '+string(wpix,ff)+'pix ANNULAE:',lines=[0,2,1],$ col=[2,3,5] loadct,0,/sil !p.position=[0.25,0.6,.5,.88] xx=shrink(findgen(nx)-xc,f=4) yy=shrink(findgen(ny)-yc,f=4) tvplot2,alog(shrink(sigma_ima>0.001<1,f=4)),/nonew,/asp,/noe,xx=xx,yy=yy tek_color endif !p.position=0 if(keyword_set(stop0)) then stop end