;+ ; NAME: ; mbmclean ; PURPOSE: ; Remove a multi-beam set of PSFs from an image via the ``clean'' algorithm. ; DESCRIPTION: ; CATEGORY: ; ALFA data processing ; CALLING SEQUENCE: ; mbmclean,dcube,psf,xloc,yloc,maxdist,iters,new_dcube, $ ; CH1=ch1,CH2=ch2,DISPLAY=display,GAIN=gain,VERBOSE=verbose ; INPUTS: ; dcube - Original source data cube structure to be cleaned. ; dcube.d(nRA,nDec,ncz) are flux densities ; dcube.bm(nDec) records the ALFA pixel used to ; acquire each row ; psf - set of nBm PSF images at same sampling resolution as image. ; xloc - X location of "object" ; yloc - Y location of "object" ; maxdist - Maximum distance from xyloc to look for local max. ; iters - Number of cleaning iterations to perform. ; OPTIONAL INPUT PARAMETERS: ; KEYWORD INPUT PARAMETERS: ; CH1 - lowest velocity channel to be cleaned ; CH2 - highest velocity channel to be cleaned ; DISPLAY - Display intermediate results, if 0, this is suppressed, ; if non-zero, this is the interval for the display, that is, ; DISPLAY=10 would cause a display every 10th iteration. ; GAIN - "Gain" of the clean process, the default value is 0.05 and ; is the scaled amount of the psf removed at each step. ; VERBOSE - Verbose printout of intermediate steps to the screen. Just ; like display, VERBOSE=0 suppresses output, VERBOSE=n will ; print information every nth iteration. ; OUTPUTS: ; new_dcube - Clean-ed image result. Pol 0 has stokes I, pol 1 has ; remains of the original image after clean-ed image is removed. ; KEYWORD OUTPUT PARAMETERS: ; COMMON BLOCKS: ; SIDE EFFECTS: ; RESTRICTIONS: ; PROCEDURE: ; MODIFICATION HISTORY: ; 1/26/94, written by Marc W. Buie, Lowell Observatory, algorithmic insight ; graciously provided by Tod Lauer (NOAO, Tucson). ; 2/25/94, MWB, added GAIN and CROSS keywords. ; 5/19/94, MWB, changed psf normalization from peak to total. Should ; now be strictly flux conserving. ; 5/21/94, MWB, Changed CROSS to allow for only a portion of psf for ; convolution. ; 10/17/2005, GLH, begun modification for 7 beam ALFA data ; 7/16/2007, GLH, version ready for production runs ;- pro mbmclean,dcube,psf,xloc,yloc,maxdist,iters,new_dcube, $ CH1=ch1,CH2=ch2,DISPLAY=display,GAIN=gain,VERBOSE=verbose print,'Enter name for output save file: ' outname=string(64) read,outname ; Define 2D circular gaussian beam for convolution with result ; hardwired for 41X41 map, 3.9' FWHM gssbm=fltarr(41,41) sigmasq=(3.9^2)/(8*alog(2)) for ix=0,40 do begin for iy=0,40 do begin sepsq=((ix-20.0)^2 + (iy-20.0)^2)/4 ; squared separation from center in arcmin gssbm[ix,iy]=exp(-sepsq/(2*sigmasq)) endfor endfor bmtot=total(gssbm) ; Set up the new dcube new_dcube = dcube new_dcube.d[*,*,*,*] = 0.0 s=size(dcube.d) ncz = s[1] npol = s[2] imw = s[3] imh = s[4] if (not keyword_set(ch2)) then ch2=ncz-1 if (not keyword_set(gain)) then gain=0.05 ; Set up information on the psf array. npsf = psf psfnorm = fltarr(7) for ibm = 0,6 do begin npsf[*,*,ibm]= psf[*,*,ibm] / max(psf[*,*,ibm]) psfnorm[ibm] = total(npsf[*,*,ibm]) endfor s=size(psf) psfw = s[1] psfh = s[2] boxm,npsf,psfw/2,psfh/2,psfw/2,psfh/2,psfx,psfy for iz = ch1,ch2 do begin ; Set up information on the input image and make a working copy. work = reform((dcube.d[iz,0,*,*]*dcube.w[iz,0,*,*]+dcube.d[iz,1,*,*]*dcube.w[iz,1,*,*])$ /(dcube.w[iz,0,*,*]+dcube.w[iz,1,*,*])) print,'Velocity channel ',iz cleanit=1 for i=0L,iters do begin if (cleanit) then begin kernal=gssbm workimg = work workc = convol(workimg,kernal) boxm,workc,xloc,yloc,maxdist,maxdist,xpos,ypos,/ABSMAX peak = workimg[xpos,ypos] if(peak lt 0.0) then $ print,'Iter ',i,' Negative at ',xpos,ypos if(peak lt 0.002) then cleanit=0 ibm=dcube.bm[ypos] scale_fac = peak * gain ; Increment the output image using beam ibm appropriate to ypos. new_dcube.d[iz,0,xpos,ypos] = new_dcube.d[iz,0,xpos,ypos] + scale_fac*psfnorm[ibm] ; Determine the sub-array for subtraction il = max([0,xpos-psfx[ibm]]) pl = max([0,psfx[ibm]-xpos]) ir = min([imw-1,xpos+(psfw-1-psfx[ibm])]) pr = min([psfw-1,psfx[ibm]+(imw-1-xpos)]) ib = max([0,ypos-psfy[ibm]]) pb = max([0,psfy[ibm]-ypos]) it = min([imh-1,ypos+(psfh-1-psfy[ibm])]) pt = min([psfh-1,psfy[ibm]+(imh-1-ypos)]) ; Subtract scale_fac at the max location using appropriate beam at each dec for iy = ib,it do begin py = pb + iy - ib jbm=dcube.bm[iy] work[il:ir,iy] = work[il:ir,iy] - npsf[pl:pr,py,jbm]*scale_fac endfor if keyword_set(verbose) then begin if i mod verbose eq 0 then begin print,'#',i,' bm ',ibm,' x,y ',xpos,ypos,' sf ',scale_fac*1000,' (',il,':',ir,',', $ ib,':',it,') (',pl,':',pr,',',pb,',',pt,')', $ peak*1000,work[xpos,ypos]*1000, $ format='(a,i6.6,a,i3,a,i3,1x,i3,a,f7.4,a,i3.3,a,i3.3,a,' + $ 'i3.3,a,i3.3,a,i3.3,a,i3.3,a,i3.3,a,i3.3,a,1x,f7.4,1x,f7.4)' endif endif ; Optional display if display ne 0 then begin if i mod display eq 0 then begin mn = min(work[il:ir,ib:it]) mx = max(work[il:ir,ib:it]) sz=20 ; 35 ; 20 zf= 6 ; 4 ; 6 new_wid = (2*sz+1)*zf tmp1 = work[xloc-sz:xloc+sz,yloc-sz:yloc+sz] tmp2 = new_dcube.d[iz,0,xloc-sz:xloc+sz,yloc-sz:yloc+sz] if !d.x_size lt new_wid or !d.y_size lt new_wid then $ setwin,0,xsize=new_wid*2,ysize=new_wid ypos = (!d.y_size - new_wid)/2 xpos1 = (!d.x_size)/2 - new_wid xpos2 = (!d.x_size)/2 tvscl,rebin(tmp1,new_wid,new_wid,/sample),xpos1,ypos tvscl,rebin(tmp2,new_wid,new_wid,/sample),xpos2,ypos endif endif endif ; cleanit endfor ; iters loop workimg = reform(new_dcube.d[iz,0,*,*]) workc = convol(workimg,gssbm) new_dcube.d[iz,0,*,*] = workc/bmtot new_dcube.d[iz,1,*,*] = work endfor ; iz loop save,new_dcube,filename=outname print,'Saved ',outname end