Date - 2000.02.10
Tester - S.T. Myers (AOC)
Platform - Linux (kernow)
Version - Weekly

Note: bugs are denoted by lines with prefix:
>>>BUG
and queries are denoted by lines with prefix:
>>>QUERY

-------------------------------------------------------------------------------

GOAL: Make a set of functions to perform automapping operations

1. Function to make large map, and find peak and rms

   Read the data into a measurement set

include 'synthesis.g'
dowait:=T 
m:=fitstoms(msfile='1608.ms',fitsfile='F16156554.UVF');
m.close() 

   Set up some preliminary parameters

cell_siz := dq.quantity('0.05arcsec')
imag_bg := 2048
imag_sm := 512
prefix := '1608'

   Define function to make dirty map

#
# Function to make dirty map of measurement set with prefix
# Returns name of dirty image
#
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()
   return img_name
}

   Run function on data

img_name := makedirty(prefix,'.bg',cell_siz,imag_bg)

   This took 4m13s, which is much longer than expected given previous
   timing.  Will probably have to go back to low-resolution initial
   image.

   Try this (after deleting 1608.bg.*)

img_name := makedirty(prefix='1608',qual='.bg',cell_siz='0.1arcsec',imag_siz=1024)

   Much faster.  Ran in 34s!

print img_name

   Output: 1608.bg.dirty

   Now make function to reture peak value, direction, and rms of image
   (do the rms with an additional function):

#
# 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)
   im.statistics(statsout=img_stat,list=T)
   img_max := img_stat.maxpos
   img_dir := im.coordmeasures(pos=img_max).direction
   img_flx := img_stat.max
   # Get rms from an outer region
   img_rms := rmsbox( im )
   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 return rms in an outer region of an image
# Pass IMAGE object img
#
function rmsbox( img )
{
   regoff := img.shape()/4; regoff[3:4] := [1,1];
   reglen := img.shape()/6; reglen[3:4] := [0,0];
   reg_blc := regoff - reglen
   reg_trc := regoff + reglen;
   r := drm.box(blc=reg_blc,trc=reg_trc)
   img.statistics(statsout=noise_stat,list=T,region=r)
   rms := noise_stat.rms
   return rms
}

   Run this on the large dirty image
   
stat_rec := imagestats( img_name )
print stat_rec

   Output: [flux=[value=0.0351882838, unit=Jy], maxpos=[311 343 1 1] , 
            direction=[type=direction, refer=J2000, m1=[value=1.14391302, 
            unit=rad], m0=[unit=rad, value=-2.05411054]], 
            rms=[value=0.00212037354, unit=Jy]]

   These can now easily be accessed, eg. stat_rec.rms returns the rms
   quantity

ms_name := spaste(prefix,".ms")
imgr := imager(ms_name)

   Now a function to clean the residual image. This reuses the clean
   model from a previous run, if possible, but destroys the previous
   mask (clean box).  Eventually keep a list of masks used.

#
# Function to clean at peak of a residual image
# Works on existing imager object imgr
#
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
}
 
   Now run cleanresid but centered on the peak found from the big
   dirty image.  Keep running cleanresid until done.

phdir := stat_rec.direction
resid_rec := cleanresid(ref imgr,prefix='1608',qual='.sm',cell_siz='0.05arcsec',imag_siz=512,phase_dir=phdir)

   Peak flux  [value=0.0354564562, unit=Jy]  Jy at  16:09:13.953 +065.32.28.998  at  [257 257 1 1] 
   Map rms  [value=0.00208935607, unit=Jy]
   Making clean mask blc,trc =  [251 251 1 1]  [263 263 1 1] 
   Max residual  [value=0.0147339571, unit=Jy]  Jy at  16:09:14.073 +065.32.27.048  at  [242 218 1 1] 
   Map rms  [value=0.00126423652, unit=Jy]

resid_rec := cleanresid(ref imgr,prefix='1608',qual='.sm',cell_siz='0.05arcsec',imag_siz=512,phase_dir=phdir)

   Using existing clean model table  1608.sm.clean
   Peak flux  [value=0.0147339571, unit=Jy]  Jy at  16:09:14.073 +065.32.27.048  at  [242 218 1 1] 
   Map rms  [value=0.00126423652, unit=Jy]
   Making clean mask blc,trc =  [236 212 1 1]  [248 224 1 1] 
   Max residual  [value=0.0130313775, unit=Jy]  Jy at  16:09:14.081 +065.32.28.498  at  [241 247 1 1] 
   Map rms  [value=0.00100806146, unit=Jy]

resid_rec := cleanresid(ref imgr,prefix='1608',qual='.sm',cell_siz='0.05arcsec',imag_siz=512,phase_dir=phdir)

   Using existing clean model table  1608.sm.clean
   Peak flux  [value=0.0130313775, unit=Jy]  Jy at  16:09:14.081 +065.32.28.498  at  [241 247 1 1] 
   Map rms  [value=0.00100806146, unit=Jy]
   Making clean mask blc,trc =  [235 241 1 1]  [247 253 1 1] 
   Max residual  [value=0.00542213256, unit=Jy]  Jy at  16:09:11.908 +065.32.16.397  at  [511 5 1 1] 
   Map rms  [value=0.000784952019, unit=Jy]

   Note that peak residual is now at map edge, so we are done.  There
   is still one lens image to be found, however.  Note that we should
   have restricted imagestat to the inner quarter of the image.  Get
   rid of that bothersome extra function getrms in the process:

#
# 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') ]
}

   In the meantime, also modify makedirty to work similarly to
   cleanresid

#
# 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
}

   We can put this all together into the following script:

#
# Run aips++ on 1608 test data for automapping emulation
#
# Read data into measurements set
#
include 'synthesis.g'
dowait:=T 
m:=fitstoms(msfile='1608.ms',fitsfile='F16156554.UVF');
m.close() 
#
# Read in functions to use
#
include 'makedirty.g'
include 'imagestats.g'
include 'cleanresid.g'
#
# 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
#

   This is what comes out:

Peak in dirty image  [value=0.0351882838, unit=Jy]  Jy at  16:09:13.953 +065.32.28.998  at  [311 343 1 1] 
Map rms  [value=0.00212037354, unit=Jy]
Peak flux  [value=0.0354564562, unit=Jy]  Jy at  16:09:13.953 +065.32.28.998  at  [257 257 1 1] 
Map rms  [value=0.00208935607, unit=Jy]
Making clean mask blc,trc =  [251 251 1 1]  [263 263 1 1] 
Max residual  [value=0.0147339571, unit=Jy]  Jy at  16:09:14.073 +065.32.27.048  at  [242 218 1 1] 
Map rms  [value=0.00126423652, unit=Jy]
Using existing clean model table  1608.sm.clean
Peak flux  [value=0.0147339571, unit=Jy]  Jy at  16:09:14.073 +065.32.27.048  at  [242 218 1 1] 
Map rms  [value=0.00126423652, unit=Jy]
Making clean mask blc,trc =  [236 212 1 1]  [248 224 1 1] 
Max residual  [value=0.0130313775, unit=Jy]  Jy at  16:09:14.081 +065.32.28.498  at  [241 247 1 1] 
Map rms  [value=0.00100806146, unit=Jy]
Using existing clean model table  1608.sm.clean
Peak flux  [value=0.0130313775, unit=Jy]  Jy at  16:09:14.081 +065.32.28.498  at  [241 247 1 1] 
Map rms  [value=0.00100806146, unit=Jy]
Making clean mask blc,trc =  [235 241 1 1]  [247 253 1 1] 
Max residual  [value=0.00460882764, unit=Jy]  Jy at  16:09:13.687 +065.32.32.748  at  [290 332 1 1] 
Map rms  [value=0.000784952019, unit=Jy]
T 

   This works pretty well!  Next step will be to add selfcalibration,
   and multiple big-small field iterations.  To do that I need to work
   out whether to do it as multiple fields.

   Below are the final scripts to load the functions to carry out the
   above sequence:

-------------------------------------------------------------------------------

#
# 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
}

-------------------------------------------------------------------------------