pro clean,uv_ij,vis_ij,inp,vsl,cmap,cln,dirty=dirty,show=show, clbox=clbox, $
dmap = dirty_map
m = inp.imsiz
niter = inp.niter
clgain = inp.clgain
clstop = inp.clstop
f_GHz = vsl.f_GHz
t = vsl.t
;--------------------------
print,'CLEAN_FTmap: uv-gridding, dirty map/beam'
fft_uv2xy, f_GHz, m, vis_ij, uv_ij, xyint, dbeam, dmap
dirty_map=dmap ; keep for later plotting
; NB: while dmap is m x m arrary, dbeam has double size, i.e., 2mx2m
; this is so to do better cleaning (JL).
;--------------------------
print,'CLEAN_Cbijm: get clean beam'
;case 1 of
;dprd_m ge 48. : cbm = clean_beam(dbeam(m/2:3*m/2-1, m/2:3*m/2-1))
;dprd_m lt 48. : cbm = clean_bijm(uv_ij,dbeam(m/2:3*m/2-1, m/2:3*m/2-1))
;endcase
;cbm = clean_bijm(uv_ij,dbeam(m/2:3*m/2-1, m/2:3*m/2-1))
cbm = clean_beam(dbeam[m/2:3*m/2-1, m/2:3*m/2-1])
Bmin = min(cbm[0:1])*xyint
Bmaj = max(cbm[0:1])*xyint
Pang = cbm[2]
if keyword_set(dirty) then begin
cmap = dmap/(!pi*cbm[0]*cbm[1]) ; its dmap but just scaled to TB [MK]
TBmax = max(cmap)
window,2,xs=512,ys=512
clev = exp(alog(.05)+alog(.99/.05)*findgen(7)/6.)
x = (findgen(m)-m/2)*xyint
;---------plot--------
loadct,5
contour,cmap,x,x,xst=1,yst=1,nlev=22,/fill
contour,cmap,x,x,xst=1,yst=1,lev=clev[3:*],/noera, $
xtit='x [arcsec]',ytit='y [arcsec]',tit='Dirty Map', c_col=[0]
;------organize-------
cboxstr='No CLEAN BOX for a DIRTY map'
cln={TBmax:TBmax, f_GHz:f_GHz, maxitr:0,gain:0.0,$
stop:100,box:[0,0,0,0],cboxstr:cboxstr, $
xyint:xyint, Bmaj:Bmaj, Bmin:Bmin, Pa:Pang}
print,' Clean box [pixels] : ',cln.cboxstr
print,' Cbeam: Bmaj [arcsec] : ',cln.Bmaj
print,' Bmin [arcsec] : ',cln.Bmin
print,' PA [degree] : ',cln.PA
print,' Max iteration : ',cln.maxitr
print,' Gain/Clean : ',cln.gain
print,' Clean stop at [%] : ',cln.stop
print,' Maximum TB [M K] : ',cln.TBmax
print,' Centr freq [GHz] : ',cln.f_GHz
print,' Map Pxlsiz [arcsec] : ',cln.xyint
print,' '
return
endif
;-------------------------------------
print,'CLEAN: Choose clean box(es)'
;-------------------------------------
; if needed, plot dirty/clean beams to assist choosing clbox
; get the cbeam from clean_only.pro
;window,1,xs=512,ys=512
;contour,dbeam,xst=5,yst=5,nlev=10,/fill
;contour,cbeem,xst=1,yst=1,lev=[0.27],/noera, $
; tit='Dirty beam & Clean beam'
if NOT keyword_set(clbox) then begin
start:
ikeywd =0
window,0,xs=512,ys=512
!p.multi=0
!x.margin=0
!y.margin=0
device,decompo=0
loadct,0
dimage=bytscl(congrid(dmap,512,512,/cubic))
tv,dimage
contour,dimage,xst=5,yst=5,/nodata,/noerase,/follow
xyouts,51,470,'Dirty map: Choose clbox',col=255,/data
print,' '
print,'To set a CLEAN box, '
print,'-Left mouse button: move the box by dragging'
print,'-Middle mouse button: resize the box by dragging'
print,'-Right mouse button: exit'
box_cursor,lbx,lby,ncx,ncy
clbox=round(([lbx,lby,lbx+ncx,lby+ncy]-256.)*float(m)/512.)
again:
oplot,[lbx,lbx+ncx],[lby,lby]
oplot,[lbx+ncx,lbx+ncx],[lby,lby+ncy]
oplot,[lbx+ncx,lbx],[lby+ncy,lby+ncy]
oplot,[lbx,lbx],[lby+ncy,lby]
print,'Choose: 1. Select another'
print,' 2. Well-done'
print,' 3. Start all over'
read,ikeywd,prompt='Which option [1,2,3]?: '
if (ikeywd lt 1) or (ikeywd gt 3) then begin
print,'Wrong option no. Try again'
goto, again
endif
if ikeywd eq 2 then goto, gonow
if ikeywd eq 3 then goto, start
; otherwise ikeywd eq 1 so that we choose another
print,'To try another box,'
print,'-Left mouse button: move the box by dragging'
print,'-Middle mouse button: resize the box by dragging'
print,'-Right mouse button: exit'
box_cursor,lbx,lby,ncx,ncy
clbox_ =round(([lbx,lby,lbx+ncx,lby+ncy]-256.)*float(m)/512.)
clbox =[clbox,clbox_]
goto,again
gonow: ; CLBOX(es) are set. Proceed . . .
endif
;--------------
; Do CLEAN
;--------------
print,' '
print,' CLEAN: initial mapping by cleaning'
print,' =================================='
cl={a:cbm[0], b:cbm[1], p:cbm[2], box:clbox, $
niter:niter, gain:clgain, stop:clstop, xyint:xyint}
if NOT keyword_set(show) then begin
clean_only, dmap, dbeam, cl, cmap, cbeam, pmap, maxitr
endif else begin
clean_only, dmap, dbeam, cl, cmap, cbeam, pmap, maxitr, show=show
endelse
cbimm = fltarr(m,m) & ncb=(n_elements(cbeam(*,0))-1)/2
cbimm(m/2-ncb:m/2+ncb,m/2-ncb:m/2+ncb)=cbeam
TBmax = max(cmap)
cmap = cmap/TBmax
;--------------------------------------
; P l o t f i n a l r e s u l t s
;--------------------------------------
big_tit ='I M A G R / C L E A N'
tit3=['DIRTY MAP','DIRTY/CLEAN BEAM','CLEAN MAP']
map3=fltarr(m,m,3)
map3[*,*,0]=dirty_map
map3[*,*,1]=dbeam(m/2:3*m/2-1, m/2:3*m/2-1)
map3[*,*,2]=cmap
winsiz =max(abs(clbox*xyint))*2.5 < m/2 * xyint ; [arcsec]
;wdelete,0
plot_3maps, map3, tit3, big_tit, xyint, winsiz
; if you want to compare obs vs cleaned visibilities
if keyword_set(show) then begin
TB2flux = 3.6E-5 * (xyint * f_ghz)^2
uvint = 206265. / (xyint * m)
iu = float(uv_ij) / uvint + m / 2
iv = imaginary(uv_ij) / uvint + m / 2
cvis_ij = fft_xy2uv(cmap * tbmax, TB2flux, iu, iv)
plot_vis2vis_simpl, vis_ij, cvis_ij
endif
;----------------------------
; Organize outputs
;----------------------------
cboxstr= strcompress(string(cl.box(0)),/remove_all) + $
' '+strcompress(string(cl.box(1)),/remove_all) + $
' '+strcompress(string(cl.box(2)),/remove_all) + $
' '+strcompress(string(cl.box(3)),/remove_all)
Bmin = min([cl.a,cl.b])*xyint
Bmaj = max([cl.a,cl.b])*xyint
Pang = cl.p
cln={TBmax:TBmax, f_GHz:f_GHz, maxitr:maxitr,gain:clgain,$
stop:clstop,box:clbox,cboxstr:cboxstr, $
xyint:xyint, Bmaj:Bmaj, Bmin:Bmin, Pa:Pang, t:t}
print,' Clean box [pixels] : ',cln.cboxstr
print,' Cbeam: Bmaj [arcsec] : ',cln.Bmaj
print,' Bmin [arcsec] : ',cln.Bmin
print,' PA [degree] : ',cln.PA
print,' Max iteration : ',cln.maxitr
print,' Gain/Clean : ',cln.gain
print,' Clean stop at [%] : ',cln.stop
print,' Maximum TB [M K] : ',cln.TBmax
print,' Centr freq [GHz] : ',cln.f_GHz
print,' Map Pxlsiz [arcsec] : ',cln.xyint
print,' '
return
end