pro make_sigma_image,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 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 read from fits-files' print,'-----------------------------------------------------------------------------' print,'Input:' print,' image_fits = name of the fits file containing data image (unit=Mjy/sr)' print,' mask_fits = -- " -- mask image (0-good pixel)' print,' wt_fits = -- " -- weight image (#frames/pixel*10)' print,'' print,'Output:' print,' sigma_fits = name of fits file where sigma image is written (unit=MJy/sr)' print,'Output via keywords:' print,' image_new_fits = fitsfile containing image cleaned of NaNs' print,' mask_new_fits = fitsfile containing masks which exclude also 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] -> 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,' test=1 --> use ngc1097 3.6micron uncut images as test data' print,' test=2 --> use ngc1097 4.5micron uncut images as test data' 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,test=1,fudg=2 ; ngc1097 test-data' print," make_sigma_image,'ngc4594.newcut.fits','ngc4594.newmask.fits',$" print," 'ngc4594.newwt.fits','ngc4594.newsigma.fits',fudg=2,check=[517.579,451.251]" print,'' 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,' HS 280809' print,'-----------------------------------------------------------------------------' 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 ide='sigma_fits' ima=readfits(image_fits,h1) wt=readfits(wt_fits) mask=readfits(mask_fits,h2) 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') ; skydrkzb=sxpar(h1,'skydrkzb') 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. 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. ;total magnitude f0=280.9 ;3.6 mu zeropoint in Jy from http://ssc.spitzer.caltech.edu/irac/calib/ if(band eq 2) then f0=179.7 pix=0.75 arcsec2_to_sr=(!dpi/180./3600.)^2 fluxtot=total(ima(indok)*1e6)*pix^2*arcsec2_to_sr mag=-2.5*alog10(fluxtot/f0) mu0=-2.5*alog10(1e6*pix^2*arcsec2_to_sr/f0) ; zeropoint 1MJy/sr corresponds to mu0 mag/arcsec2 ; mag=-2.5*alog10(total(ima(indok)))+mu0 ; yields the same print,'total magnitude ',mag,' in band ',band ;============================================================================= ;convert Mjy/sr -> orig. digital units/second -> digital units -> electrons dnsec=(ima+backgrnd)/fluxconv dntot=dnsec*framtime*weight ele=dntot*gain>1d-12 ;sigma of electrons: add Poisson noise and RON quadratically sigma2_ele=ele+weight*ronoise^2 sigma_ele=sqrt(sigma2_ele) ;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 without Ron (just for test) sigma2_ele0=ele sigma_ele0=sqrt(sigma2_ele0) sigma_dnsec0=sigma_ele0/(gain*framtime*weight) sigma_ima0=sigma_dnsec0*fluxconv ;--------------------------------------------------------------------- ;make sure 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(cbad ge 1) then begin ima(indbad)=-1. sigma_ima(indbad)=-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 ne data image? Does not have NaNs if(keyword_set(ima_new_fits)) then begin sxaddhist,'data-image modified by make_sigma_image.pro '+!stime,h11,/comment writefits,ima_new_fits,mask,h11 endif save,file=sigma_fits+'.save',ima,sigma_ima,mask,mu0,fudge ;---------------------------------------------------------------------- ;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] endif if(keyword_set(stop0)) then stop end