Date - 2000.02.08
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

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

1. Test goal: read in MS 1608.ms, produce dirty image, find peak of
   dirty image (continued)

   Start over, remake ms from UVF
#
# Read UVF into MS
#
m:=fitstoms(msfile='1608.ms',fitsfile='F16156554.UVF');
m.close() 
#
# Run imager on MS to make low-res dirty image
#
imgr:=imager('1608.ms')
imgr.setimage(cellx='0.2arcsec', celly='0.2arcsec', nx=512, ny=512, 
  stokes='I', spwid=1, fieldid=1)	
imgr.setdata(spwid=1, mode="none")	
imgr.uvrange(uvmin=0.0, uvmax=258000.0)	
imgr.weight('robust')
bmaj:=F; bmin:=F; bpa:=F;	
imgr.image(type='psf',image='1608.200mas.psf')	
imgr.fitpsf(psf='1608.200mas.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.200mas.dirty')	
#
# Use image constructor to access dirty image and find peak
#
im := image('1608.200mas.dirty')
dirty_200mas_ref := im.shape()/2; dirty_200mas_ref[3:4] := [1,1];
im.coordmeasures(pos=dirty_200mas_ref,direction=direc_200mas_ref)
print 'Center of 200mas map at ',dm.time(direc_200mas_ref.m0),dm.angle(direc_200mas_ref.m1)
im.statistics(statsout=dirty_200mas_stat,list=T) 
print 'Max value = ',dirty_200mas_stat.max,' at ',dirty_200mas_stat.maxpos
#
# Store pixel position of max, and print in sensible units
#
dirty_200mas_max := dirty_200mas_stat.maxpos
direc_200mas_max := im.coordmeasures(pos=dirty_200mas_max).direction
print 'Peak at ',dm.time(direc_200mas_max.m0),dm.angle(direc_200mas_max.m1),' at ',dirty_200mas_max
im.close()
#
# Back to imager, but shift phasecenter and remake dirty image
#
imgr.setimage(cellx='0.05arcsec', celly='0.05arcsec', nx=512, ny=512, 
  stokes='I', spwid=1, fieldid=1, doshift=T, phasecenter=direc_200mas_max)
imgr.uvrange()
imgr.weight('robust')
imgr.setdata(spwid=1, mode="none")
bmaj:=F; bmin:=F; bpa:=F; imgr.image(type='psf',image='1608.50mas.psf')
imgr.fitpsf(psf='1608.50mas.psf', bmaj=bmaj, bmin=bmin, bpa=bpa) 
print 'Bmaj = ', bmaj;	print 'Bmin = ', bmin;  print 'Bpa = ', bpa;
imgr.image(type='observed',image='1608.50mas.dirty')
#
# Look at high-res dirty image
#
im := image('1608.50mas.dirty')
dirty_50mas_ref := im.shape()/2; dirty_50mas_ref[3:4] := [1,1];
im.coordmeasures(pos=dirty_50mas_ref,direction=direc_50mas_ref)
print 'Center of 50mas map at ',dm.time(direc_50mas_ref.m0),dm.angle(direc_50mas_ref.m1)
im.statistics(statsout=dirty_stat,list=T)
print 'Max value = ',dirty_stat.max,' at ',dirty_stat.maxpos

   Output:
Bmaj =  [value=0.5, unit=arcsec]
Bmin =  [value=0.25, unit=arcsec]
Bpa =  [value=57.296, unit=deg]
Center of 200mas map at  16:09:10.732 +065.32.45.800
Max value =  0.03593  at  [156 172 1 1] 
Peak at  16:09:13.953 +065.32.28.998  at  [156 172 1 1] 
Bmaj =  [value=0.30155, unit=arcsec]
Bmin =  [value=0.19959, unit=arcsec]
Bpa =  [value=296.38, unit=deg]
Center of 50mas map at  16:09:13.961 +065.32.28.948
Max value =  0.04142  at  [257 256 1 1] 

>>>BUG: THIS DOESNT WORK WITHOUT THE EXPLICT RESETTING OF THE WEIGHTS
        WITH IMAGER.WEIGHT IN SPITE OF THE IMAGER.UVRANGE.  SETTING
        UVRANGE SHOULD RESET THE WEIGHTS CORRECTLY (PERHAPS FORCE A
        WEIGHT?), BUT I SEE THAT THERE IS A PROBLEM WITH ONLY ONE 
        WEIGHT COLUMN (IT WOULD HAVE TO BACK OUT THE PREVIOUS
        UVRANGE).  SOLUTION - HAVE UVRANGE SET A SEPARATE UV MASK
        (TABLE) THAT GETS APPLIED WHEN THE UVDATA IS ACCESSED?

   This task now seems to work.

2. Test goal: compare efficiency of doing low-resolution big image for
   object finding versus a very large full resolution map.  The
   critical thing will be imaging speed.

   Aconfig/Xband shapshots 0.05" cells, cover primary beam (5.3' FWHM)
   Should reasonably to out to a radius of 120" from the center:

             240"/0.05" = 4800 pix => 4096 image @ 0.05" = 204.8" sq

imgr:=imager('1608.ms')
imgr.setimage(cellx='0.05arcsec', celly='0.05arcsec', nx=4096, ny=4096, 
  stokes='I', spwid=1, fieldid=1)
imgr.setdata(spwid=1, mode="none")
imgr.weight('robust')

   This takes 37sec to do weights.  (On 128MB RAM Linux box kernow)
   Disk usage: 1.26MB for 1608.ms cf. 313KB for F16156554.UVF

imgr.uvrange()	
bmaj:=F; bmin:=F; bpa:=F; imgr.image(type='psf',image='1608.large.psf')

   Uses tiled gridding (8192 tiles, 814 accesses, 175R/175W).
   Took 454.65 sec = 7.6 min for imager.image
   Disk usage: 65.8MB for 1608.large.psf

imgr.fitpsf(psf='1608.large.psf', bmaj=bmaj, bmin=bmin, bpa=bpa)

   Took 56sec for fitpsf.
   Output: Beam fit: 0.322863 by 0.206411 (arcsec) at pa 125.169 (deg) 

imgr.image(type='observed',image='1608.large.dirty')

   Took 452 sec for imager.image
   Hey, in the middle of this some "Progress Meter" popped up by
   itself and told me all tools were quiet.  Huh?

im := image('1608.large.dirty')
dirty_large_ref := im.shape()/2; dirty_large_ref[3:4] := [1,1];
direc_large_ref := im.coordmeasures(pos=dirty_large_ref).direction
print 'Center = ',dm.time(direc_large_ref.m0),dm.angle(direc_large_ref.m1),' at ',dirty_large_ref

   Output: Center =  16:09:10.708 +065.32.45.950  at  [2048 2048 1 1]

im.statistics(statsout=dirty_large_stat,list=T)
print 'Max value = ',dirty_large_stat.max,' at ',dirty_large_stat.maxpos

   Output: Max value =  0.034342  at  [1644 1708 1 1]

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

   Output: Peak at  16:09:13.961 +065.32.28.948  at  [1644 1708 1 1]

im.close()

   The time on this was pretty long, especially if I am to process the
   CLASS data for example!  Lets try 2048 sq images...

imgr.setimage(cellx='0.05arcsec', celly='0.05arcsec', nx=2048, ny=2048, 
  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')

   Took 106sec for imager.image

imgr.fitpsf(psf='1608.large.psf', bmaj=bmaj, bmin=bmin, bpa=bpa)

   Took 6 sec for fitpsf
   Output: Beam fit: 0.318515 by 0.206957 (arcsec) at pa 124.022 (deg)

imgr.image(type='observed',image='1608.large.dirty')

   Took 107sec

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

   Output: Peak at  16:09:13.961 +065.32.28.948  at  [620 684 1 1]

im.close()
imgr.close()
imgr.done()

   Go back and make small image.
>>>QUERY: YOU MUST CLOSE AND KILL THE OLD IMAGER AND CONSTRUCT A NEW
          ONE, OR THE NEW MAP WILL BE AT THE OLD PHASECENTER. IS THIS
          ON PURPOSE, OR A BUG?

imgr:=imager('1608.ms')
imgr.setimage(cellx='0.05arcsec', celly='0.05arcsec', 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)

   Output: Beam fit: 0.301553 by 0.199593 (arcsec) at pa 296.384 (deg)

imgr.image(type='observed',image='1608.sm.dirty')
im := image('1608.sm.dirty')
dirty_ref := im.shape()/2; dirty_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

   Output: Center =  16:09:13.969 +065.32.28.898  at  [256 256 1 1]

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

   Output: Max value =  0.035457  at  [258 258 1 1]

   CONCLUSION: THIS IS COMPETITIVE WITH MAKING A LOW-RES IMAGE FIRST.

3. Test Goal: put a clean box around the peak, and clean to some
   appropriate level.

   Get some statistics on an outer region of the dirty map

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)

   Output: Number points = 2.924100e+04  
   Sum           = -3.088177e+00        Mean     = -1.056112e-04 
   Variance      = 4.362766e-06         Sigma    = 2.088723e-03  
   Rms           = 2.091356e-03  

   Minimum value is -6.385271e-03  at [187, 105, 1, 1]
   Maximum value is 1.623204e-02   at [190, 110, 1, 1]

   Note: rms noise 2.1 mJy looks OK.

im.done()

   Use the bmaj as the box size

cell_siz := dq.quantity('0.05arcsec')
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
print blc_max, trc_max

   Output: [252 252 1 1]  [264 264 1 1

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=100,gain=0.05,
   threshold=cln_lim)

   Output: Clean used 90 iterations to get to a max residual of 0.00208398

   Now make residual image.  Note that the Getting Results implies
   that for difference imaging (a la difmap) you want to go back to
   the original UVDATA and thus should use the residual function.
   MORE OR CLEARER EXPLANATION SHOULD BE AVAILABLE.

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

   Output: Max residual  =  0.014734  at  [243 219 1 1]

   There is more to clean

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 '

   Output: Rms residual =  0.0012656  Jy

im.done()
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')

   Note that I want to clean only at the new mask location (keeping
   the lens "components" separate.  But Im not sure how to make the
   new model go into box.clean2 instead of box.clean. From the 
   documentation, this is not possible!
>>>QUERY: IS IT POSSIBLE (DESIRABLE) TO HAVE A SEPARATION OF FIXED AND
          VARIABLE MODELS IN THE CLEAN SO THAT THE RESULTS OF CLEAN
          CAN BE KEPT TRACK OF EASIER! OR ALLOW IT TO WORK FROM THE
          PREVIOUS MODEL_DATA COLUMN.

imgr.clean(mask='box_mask2',model='box.clean',niter=500,gain=0.05,threshold=cln_lim)

   Output: Clean used 143 iterations to get to a max residual of 0.00125784

   It is clear that this mostly works.  The next step is to put these
   cycles into functions and build into an automapper for this data.

imgr.done()

   BELOW IS A GLISH SCRIPT TO CARRY OUT THE ABOVE:

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

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()