function thz_clean,map,beam

   ; Display map for drawing clean box
   tvscl,map
   ; Create copies of map and beam to preserve original
   dmap = map
   dbeam = beam/max(beam)   ; Normalize beam to its own maximum
   n = n_elements(map[*,0]) ; Size of map (number of cells in one dimension)
   hn = n/2
;   clbox = [40,40,59,59]
   ; Fit Gaussian to beam center
   clparm = clean_beam(dbeam)  ; Fit gaussian to beam and return major/minor axes and orientation
   niter = 100    ; Number of components to clean
   clgain = 0.1   ; Clean gain parameter
   clstop = 10.   ; Percentage of peak at which to stop
   xyint = 1.0    ; Cell size [arcsec]
   inroi=defroi(n,n,xv,yv,/nofill) ;to draw any shape ROI (INROI is interior cells)
;   nroi=[1]               ;to keep track of the number of ROI
;   nel_roi=[n_elements(inroi)]
;   xvi=xv mod n
;   yvi=yv/n
;   verts=[[xv],[yv]]  ;to save the initial roi vertices
;   cl={a:clparm[0], b:clparm[1], p:clparm[2], inroi:inroi,nroi:nroi,box:verts,$
;       niter:niter, gain:clgain, stop:clstop, xyint:xyint}
;   clean_only, dmap, dbeam, cl, cmap, cbeam, pmap, maxitr ;, /show

   ; Ellipse parameters, for creating clean beam
   eparm =[1.0, 0.,0., clparm[0], clparm[1], -clparm[2], 0.0]
   ; Make arrays of indexes, convenient for creating clean beam
   m = 2*n-1  ; Beam size
   ii = (jj = fltarr(m,m))
   for i = 0, m-1 do ii(i,*) = float(i) - m/2
   jj = transpose(ii)
   ; Create clean beam
   cbeam = cvd_ellipse(ii,jj,eparm)
   ; Start with map
   resid = dmap
   ; Create arrays to hold clean component information (x, y, and peak)
   ix = (iy = intarr(niter))
   comp = fltarr(niter)
   mapmax = max(resid[inroi])
   mincomp = mapmax*clstop/100.  ; Minimum component at which to stop cleaning
   for i = 0, niter-1 do begin
      comp[i] = clgain*max(resid[inroi],imax) ; peak value to be cleaned
      ix[i] = inroi[imax] mod n  ; x location of component, in map
      iy[i] = inroi[imax]/n      ; y location of component, in map
      sx = ix[i] - hn            ; Amount to shift beam, in x
      sy = iy[i] - hn            ; Amount to shift beam, in y
      resid = resid - (comp[i]*shift(dbeam,sx,sy))[hn:3*hn,hn:3*hn]
      tv,bytscl(resid,max=max(dmap[inroi])),110,0 & wait,0.1
      if (max(resid[inroi]) le mincomp) then break
   endfor
   
   print,'Number of iterations = ',i
   if (i ne niter) then niter = i
   ; Restore clean components
   for i = 0, niter-1 do begin
      sx = ix[i] - hn            ; Amount to shift beam, in x
      sy = iy[i] - hn            ; Amount to shift beam, in y
      resid = resid + (comp[i]*shift(cbeam,sx,sy))[hn:3*hn,hn:3*hn]
      tv,bytscl(resid,max=max(dmap[inroi])),220,0 & wait,0.1
   endfor
   cmap = resid
   tvscl,cmap,220,0
   
return,cmap
end