#include 'synthesis.g'
dowait:=T	# Tell the Object system to wait for a method to finish
#
# Read UVF into MS
#
# m:=fitstoms(msfile='1608.ms',fitsfile='F16156554.UVF');
# m.close() 
#
# Run imager on MS to make high-res large image
#
cell_siz := dq.quantity('0.05arcsec')
#
imgr:=imager('1608.ms')
imgr.setimage(cellx=cell_siz, celly=cell_siz, nx=4096, ny=4096, 
  stokes='I', spwid=1, fieldid=1)
imgr.setdata(spwid=1, mode="none")
imgr.weight('robust')
imgr.uvrange()	
bmaj:=F; bmin:=F; bpa:=F; imgr.image(type='psf',image='1608.large.psf')
imgr.fitpsf(psf='1608.large.psf', bmaj=bmaj, bmin=bmin, bpa=bpa)
system.print.precision := 5
print 'Bmaj = ', bmaj;	print 'Bmin = ', bmin;  print 'Bpa = ', bpa;
imgr.image(type='observed',image='1608.large.dirty')
imgr.close()
imgr.done()
#
# Use image constructor to access dirty image and find peak
#
im := image('1608.large.dirty')
im.statistics(statsout=dirty_large_stat,list=T)
dirty_large_max := dirty_large_stat.maxpos
direc_large_max := im.coordmeasures(pos=dirty_large_max).direction
print 'Peak at ',dm.time(direc_large_max.m0),dm.angle(direc_large_max.m1),' at ',dirty_large_max
im.close()
#
# Go back and make small image
#
imgr:=imager('1608.ms')
imgr.setimage(cellx=cell_siz, celly=cell_siz, nx=512, ny=512, 
  stokes='I', spwid=1, fieldid=1, doshift=T, phasecenter=direc_large_max)
imgr.setdata(spwid=1, mode="none")
imgr.weight('robust')
imgr.uvrange()	
bmaj:=F; bmin:=F; bpa:=F; imgr.image(type='psf',image='1608.sm.psf')
imgr.fitpsf(psf='1608.sm.psf', bmaj=bmaj, bmin=bmin, bpa=bpa)
system.print.precision := 5
print 'Bmaj = ', bmaj;	print 'Bmin = ', bmin;  print 'Bpa = ', bpa;
imgr.image(type='observed',image='1608.sm.dirty')
#
# Look at small dirty image
#
im := image('1608.sm.dirty')
dirty_ref := im.shape()/2; dirty_50mas_ref[3:4] := [1,1];
direc_ref := im.coordmeasures(pos=dirty_ref).direction
print 'Center = ',dm.time(direc_ref.m0),dm.angle(direc_ref.m1),' at ',dirty_ref
im.statistics(statsout=dirty_stat,list=T)
dirty_max := dirty_stat.maxpos
direc_max := im.coordmeasures(pos=dirty_max).direction
print 'Max value = ',dirty_stat.max,' at ',dirty_stat.maxpos
#
# Get statistics in an outer region
#
regoff := im.shape()/4; regoff[3:4] := [1,1];
reglen := im.shape()/6; reglen[3:4] := [0,0];
reg_blc := regoff - reglen; reg_trc := regoff + reglen;
r := drm.box(blc=reg_blc,trc=reg_trc)
im.statistics(statsout=noise_stat,list=T,region=r)
im.done()
#
# Make clean box (mask) around max
#
boxlen := as_integer( dq.div(bmaj,cell_siz).value + 0.5 )
box_siz := [ boxlen, boxlen, 0, 0 ]
blc_max := dirty_max - box_siz; trc_max := dirty_max + box_siz;
imgr.boxmask(mask='box_mask',blc=blc_max,trc=trc_max)
imgr.make('box.clean')
cln_lim := dq.quantity(noise_stat.rms,'Jy')
imgr.clean(mask='box_mask',model='box.clean',niter=500,gain=0.05,
   threshold=cln_lim)
#
# Make residual image and look at it
#
imgr.residual(model='box.clean',image='1608.sm.residual')
im := image('1608.sm.residual')
im.statistics(statsout=resid_stat,list=T)
resid_max := resid_stat.maxpos
resdr_max := im.coordmeasures(pos=resid_max).direction
print 'Max residual = ',resid_stat.max,' at ',resid_stat.maxpos
#
# Get statistics in an outer region
#
regoff := im.shape()/4; regoff[3:4] := [1,1];
reglen := im.shape()/6; reglen[3:4] := [0,0];
reg_blc := regoff - reglen; reg_trc := regoff + reglen;
r := drm.box(blc=reg_blc,trc=reg_trc)
im.statistics(statsout=noise_stat,list=T,region=r)
print 'Rms residual = ',noise_stat.rms,' Jy '
im.done()
#
# Next clean iteration
#
boxlen := as_integer( dq.div(bmaj,cell_siz).value + 0.5 )
box_siz := [ boxlen, boxlen, 0, 0 ]
blc_max2 := resid_max - box_siz; trc_max2 := resid_max + box_siz;
imgr.boxmask(mask='box_mask2',blc=blc_max2,trc=trc_max2)
cln_lim := dq.quantity(noise_stat.rms,'Jy')
imgr.clean(mask='box_mask2',model='box.clean',niter=500,gain=0.05,threshold=cln_lim)
#
#  End for now
#
imgr.done()