#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() 
#
# Load functions:
#
#
# Function to get image stats
# Returns glish records with max, maxpos, maxdirec, rms
#
function imagestats( img_name  )
{
   # Use image constructor to access dirty image and find peak 
   im := image(img_name)
   # Search inner quarter
   regoff := im.shape()/2; regoff[3:4] := [1,1];
   reglen := im.shape()/4; 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=img_stat,list=T,region=r)
   img_max := img_stat.maxpos
   img_dir := im.coordmeasures(pos=img_max).direction
   img_flx := img_stat.max
   # Get rms noise level from 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)
   img_rms := noise_stat.rms
   im.close()
   im.done()
   # Return glish record, with fluxes as quanta
   return [ flux=dq.quantity(img_flx,'Jy'), maxpos=img_max, direction=img_dir, rms=dq.quantity(img_rms,'Jy') ]
}
#
# Function to make dirty map of measurement set with prefix
# This will create and destroy its own imager object.
# Returns glish record of map stats
#
function makedirty(prefix,qual,cell_siz,imag_siz)
{	
   # Use imager to make large dirty image
   ms_name := spaste(prefix,".ms")
   imgr := imager(ms_name)
   imgr.setimage(cellx=cell_siz, celly=cell_siz, nx=imag_siz, ny=imag_siz, stokes='I', spwid=1, fieldid=1)
   imgr.setdata(spwid=1, mode="none")
   imgr.weight('robust')
   imgr.uvrange()
   pre_name := spaste(prefix,qual)
   psf_name := spaste(pre_name,".psf")
   bmaj:=F; bmin:=F; bpa:=F; imgr.image(type='psf',image=psf_name)
   imgr.fitpsf(psf=psf_name, bmaj=bmaj, bmin=bmin, bpa=bpa)
   img_name := spaste(pre_name,".dirty")
   imgr.image(type='observed',image=img_name)
   imgr.close()
   imgr.done()
   # Get stats of dirty image
   img_stats := imagestats( img_name )
   peak_flx := img_stats.flux
   peak_pix := img_stats.maxpos
   peak_dir := img_stats.direction
   img_rms := img_stats.rms
   print 'Peak in dirty image ',peak_flx,' Jy at ',dm.time(peak_dir.m0),dm.angle(peak_dir.m1),' at ',peak_pix
   print 'Map rms ',img_rms
   # return glish record of map stats
   return img_stats
}
#
# Function to clean at peak of a residual image
# Works on existing imager object imgr
# Returns glish record of last residual map stats
#
function cleanresid(ref imgr,prefix,qual,cell_siz,imag_siz,phase_dir)
{	
   # Some local parameters
   cln_niter := 500
   cln_gain := 0.05
   # Center image on input phasecenter
   imgr.setimage(cellx=cell_siz, celly=cell_siz, nx=imag_siz, ny=imag_siz, stokes='I', spwid=1, fieldid=1, doshift=T, phasecenter=phase_dir)
   imgr.setdata(spwid=1, mode="none")
   imgr.weight('robust')
   imgr.uvrange()	
   pre_name := spaste(prefix,qual)
   psf_name := spaste(pre_name,".psf")
   bmaj:=F; bmin:=F; bpa:=F; imgr.image(type='psf',image=psf_name)
   imgr.fitpsf(psf=psf_name, bmaj=bmaj, bmin=bmin, bpa=bpa)
   img_name := spaste(pre_name,".dirty")
   cln_name := spaste(pre_name,".clean")
   # check if clean model already exists, else make a blank one
   if ( tableexists(cln_name) )
      {
      print 'Using existing clean model table ',cln_name
      imgr.residual(model=cln_name,image=img_name)
      } 
   else
      {
      imgr.image(type='observed',image=img_name)
      imgr.make(cln_name)
      }
   # Get image stats
   img_stats := imagestats( img_name )
   peak_flx := img_stats.flux
   peak_pix := img_stats.maxpos
   peak_dir := img_stats.direction
   img_rms := img_stats.rms
   print 'Peak flux ',peak_flx,' Jy at ',dm.time(peak_dir.m0),dm.angle(peak_dir.m1),' at ',peak_pix
   print 'Map rms ',img_rms
   # Make clean box (mask) around max
   cln_lim := img_rms
   boxlen := as_integer( dq.div(bmaj,cell_siz).value + 0.5 )
   box_siz := [ boxlen, boxlen, 0, 0 ]
   blc_max := peak_pix - box_siz
   trc_max := peak_pix + box_siz;
   mask_name := spaste(pre_name,".mask")
   if ( tableexists(mask_name) ) tabledelete(mask_name)
   print 'Making clean mask blc,trc = ',blc_max,trc_max
   imgr.boxmask(mask=mask_name,blc=blc_max,trc=trc_max)
   imgr.clean(mask=mask_name,model=cln_name,niter=cln_niter,gain=cln_gain,threshold=cln_lim)
   res_name := spaste(pre_name,".residual")
   imgr.residual(model=cln_name,image=res_name)
   # Get stats of new residual image
   img_stats := imagestats( res_name )
   peak_flx := img_stats.flux
   peak_pix := img_stats.maxpos
   peak_dir := img_stats.direction
   img_rms := img_stats.rms
   print 'Max residual ',peak_flx,' Jy at ',dm.time(peak_dir.m0),dm.angle(peak_dir.m1),' at ',peak_pix
   print 'Map rms ',img_rms
   # note that imgr now points to the imager object
   # return last glish record of residual map stats
   return img_stats
}
#
# Run aips++ on 1608 test data for automapping emulation
#
# Make big low-res dirty map and find peak
#
stat_rec := makedirty(prefix='1608',qual='.bg',cell_siz='0.1arcsec',imag_siz=1024)
#
# Check if there is something above 6 sigma
#
pk_flx := stat_rec.flux
rms_flx := stat_rec.rms
if ( pk_flx.value > 6.0*rms_flx.value ) 
   {
   # Start new imager and make small hi-res field and clean it
   imgr := imager('1608.ms')
   phdir := stat_rec.direction
   resid_rec := cleanresid(ref imgr,prefix='1608',qual='.sm',cell_siz='0.05arcsec',imag_siz=512,phase_dir=phdir)
   #
   pk_flx := resid_rec.flux
   rms_flx := resid_rec.rms
   # 
   # Iterate on this position while there is signal
   #
   while ( pk_flx.value > 6.0*rms_flx.value ) 
      {
      resid_rec := cleanresid(ref imgr,prefix='1608',qual='.sm',cell_siz='0.05arcsec',imag_siz=512,phase_dir=phdir)
      #
      pk_flx := resid_rec.flux
      rms_flx := resid_rec.rms
      }
   #
   # Done with this field
   #
   imgr.close()
   imgr.done()
}
# 
# Done with this data
#