pro clean_only, dmap, dbeam, clparm, cmap, cbeam, pmap, maxitr, show=show

Bmaj  = max([clparm.a, clparm.b]); in clean-main program you need to multiply
Bmin  = min([clparm.a, clparm.b]); these by xyint, but here you don't
Bpa   = clparm.p
cbx   = clparm.box
ni    = clparm.niter
gain  = clparm.gain*max(dmap)
rmin  = clparm.stop/100.
m     = n_elements(dmap(*,0))
xyint = clparm.xyint

;>>>>>>>>>>>>>>>>>>>>>> Set-up for graphic display (JL)

;if keyword_set(show) then plot_set_xyscan,dmap,sz,sx,sy,pk,trx,wran,cord,gost
if keyword_set(show) then ssmem_xyscan, winid, dmap, /init
;>>>>>>>>>>>>>>>>>>>>>> Now the normal pgm begins

;--------------- double the map array to match beam array
dmap0=dmap ; just for keeping record, dmap will be changed.

;--------------- locate points within the clean box

ncbox=n_elements(cbx)/4
; Cleaning
ii=fltarr(m,m)
for i=0,m-1 do ii(i,*)=float(i)-m/2
jj=transpose(ii)
box=where(ii ge cbx(0) and jj ge cbx(1)   $
       and ii le cbx(2) and jj le cbx(3))
if ncbox gt 1 then begin
 for i=1,ncbox-1 do begin
  i4=i*4
  box_= where(ii ge cbx(i4) and jj ge cbx(i4+1)   $
       and ii le cbx(i4+2) and jj le cbx(i4+3))
  box=[box,box_]
 endfor
endif

;--------------- make a clean beam
    eparm =[1.0, 0.,0., clparm.a, clparm.b, -clparm.p, 0.0]
; Watch out sign convention
    cbeam = cvd_ellipse(ii,jj,eparm) ; /(!pi*Bmaj*Bmin)*xyint^2
; Want to reduce the array size since a gaussian falls off enough after 4*aa
    nmux=fix(Bmaj*4.) < (m/2 - 1)
    cbeam=cbeam(m/2-nmux:m/2+nmux,m/2-nmux:m/2+nmux)

;--------------- subtract clean components
ipx =intarr(ni) & jpx=intarr(ni)
nolessthan =max(dmap[box])*rmin
in1=m/2 & in2=3*m/2-1 & m_1=m-1 ; frequenctly used indexes

for l=0,ni-1 do begin
 part=dmap[box]
 peak=max(part)
 there=where(part eq peak)
 ipx(l)=box(there(0)) mod m
 jpx(l)=box(there(0)) / m
 ip = ipx[l]+m/2   & jp = jpx[l] + m/2 ; [0,0] of dmap -> [m/2,m/2] of dbeam
 shfted_dbeam=shift(dbeam,ip-m,jp-m) ; shift debam rel to its center, [m,m]
 dmap=dmap-gain*shfted_dbeam(in1:in2, in1:in2);take only inner part
 if (peak lt nolessthan) then goto,sk8
;if keyword_set(show) then plot_xyscan,dmap,sx,sy,pk,trx,wran,cord=cord,gost=gost ;
if keyword_set(show) then ssmem_xyscan, winid, dmap
endfor

print, '  CLEAN: Reached maximum iteration...'
goto,sk9
sk8:
print, '  CLEAN: Reached target value...'
sk9:
maxitr=l

;--------------- construct a clean map (with background added)
cmap=dmap       ; start from the residual map
pmap=fltarr(m,m)    ; start from zero background
; NB. cbeam has size (2*numx+1, 2*nmux+1)

for l=0,maxitr-1 do begin
 i0=ipx(l) & i1=(i0-nmux)>0 & i2=(i0+nmux)<m_1 ; safe-guard againts points
 j0=jpx(l) & j1=(j0-nmux)>0 & j2=(j0+nmux)<m_1 ;  lying outside of the map
 ri1=nmux+(i1-i0) & ri2=nmux+(i2-i0) & rj1=nmux+(j1-j0) & rj2=nmux+(j2-j0)
 cmap[i1:i2,j1:j2] = cmap[i1:i2,j1:j2] + gain*cbeam[ri1:ri2,rj1:rj2]
 pmap[i0,j0] = pmap[i0,j0] + gain
;if keyword_set(show) then plot_xyscan,cmap,sx,sy,pk,trx,wran,cord=cord,gost=gost
if keyword_set(show) then $
 ssmem_xyscan, winid, cmap
endfor

; Need to raise pmap to the orig. value (recall that last 10% ignored)
ok=where(pmap gt 0.)
;pmap(ok) = pmap(ok) + nolessthan ; this didnot work.

; this uses only the cln components and so cleaner than cmap
;cmap1=fltarr(m,m)
;for l=0,n_elements(ok)-1 do begin
; i0=ok(l) mod m & i1=i0-nmux & i2=i0+nmux
; j0=ok(l) / m & j1=j0-nmux & j2=j0+nmux
; cmap1[i1:i2, j1:j2] = cmap1[i1:i2, j1:j2] + pmap[i0,j0]*cbeam
;endfor
;cmap1=cmap1/max(cmap1)

; conclusion cmap & cmap1 are identical
; for cmap, divide it by the cln beam since it is convolved with cbeam
cmap = cmap /(!pi*Bmaj*Bmin)
;cmap = cmap * xyint^2

;pmap = cmap1*max(cmap)
if keyword_set(show) then wdelete, winid
return
end