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

   Base first part of operations on Recipe 3 ss433.g

include 'synthesis.g'
dowait:=T	# Tell the Object system to wait for a method to finish
imgr:=imager('1608.ms')	# make imager object

   Now have 1608.ms as imager object.  The data is a 30sec snapshot of
   the quad lens 1608+656 from the CLASS Survey (Spring 1994), in A
   config at Xband.

   Set image parameters.  This field 25.6" x 25.6" should be
   sufficient.

imgr.setimage(cellx='0.05arcsec', celly='0.05arcsec', nx=512, ny=512, 
  stokes='I', spwid=1, fieldid=1)	
imgr.setdata(spwid=1, mode="none")
 
   Make beam, report size

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

   Here we find that the beam position angle (bpa) is 660 deg.
   Further investigation:

system.print.precision := 5
print 'Bmaj = ', bmaj;print 'Bmin = ', bmin;  print 'Bpa = ', bpa;

   Output:
Bmaj =  [value=0.30635, unit=arcsec]
Bmin =  [value=0.19783, unit=arcsec]
Bpa =  [value=660.17, unit=deg]

>>>BUG - WHY IS BPA NOT WITHIN -180 to 180 deg OR 0 - 360 deg?

bug() - bug reported

Note: difmap gives:
! Estimated beam: bmin=0.2137 arcsec, bmaj=0.3038 arcsec, bpa=-62.29 degrees
for the same mapsize and cellsize, natural weighting

In the course of investigating, difmap finds the source with an offset:
! Total accumulated eastward shift  = -20.2 (arcsec).
! Total accumulated northward shift = 17 (arcsec).
so we need to make a larger trial map of 512x512 with 0.2" cells:

imgr.setimage(cellx='0.2arcsec', celly='0.2arcsec', nx=512, ny=512, 
  stokes='I', spwid=1, fieldid=1)
imgr.setdata(spwid=1, mode="none")
bmaj:=F; bmin:=F; bpa:=F; imgr.image(type='psf',image='1608.psf')
imgr.fitpsf(psf='1608.psf', bmaj=bmaj, bmin=bmin, bpa=bpa) 

  Beam fit: 0.5 by 0.25 (arcsec) at pa 57.2958 (deg) 

  Why is the beam 0.5 x 0.25 with such heavy uvcutoff?

  difmap gives:
! Map grid = 512x512 pixels with 0.200x0.200 arcsec cellsize.
! Only data in the UV range: 0 -> 258 (kilo-wavelengths) will be gridded.
! Inverting map and beam 
! Your chosen uvrange limits excluded 63% of the data.
! Estimated beam: bmin=0.5997 arcsec, bmaj=0.6992 arcsec, bpa=-40.68 degrees

  Try again after deleting 1608.psf (it should have overwritten OK)
  No change noted.
  Maybe it requires an explicit uvrange (it should be able to handle
  this).

imgr.uvrange(uvmin=0.0, uvmax=258000.0)
bmaj:=F; bmin:=F; bpa:=F; imgr.image(type='psf',image='1608.psf')
imgr.fitpsf(psf='1608.psf', bmaj=bmaj, bmin=bmin, bpa=bpa) 

  Beam fit: 0.666492 by 0.618333 (arcsec) at pa -33.5142 (deg)
  
  This seems to work reasonably.

>>>BUG: IMAGER SHOULD BE ABLE TO SELF-SELECT UVDATA TO GRID GIVEN
   NX,NY,CELLX,CELLY

bug() - bug reported

   Continue, now make dirty image at low resolution

imgr.image(type='observed',image='1608.200mas.dirty')

   Done with imager for now, will use image package on dirty map.

imgr.close()

   Note that doing this from the GUI seems to crash glish still!
   At least it works in glish...

im := image('1608.200mas.dirty')
im.statistics(statsout=dirty_stat,list=T)

   The desired stats should now be in glish list dirty_stat
   Get the position of the min,max

print 'Max value = ',dirty_stat.max,' at ',dirty_stat.maxpos

   Output: Max value =  0.041276  at  [156 172 1 1]

im.coordmeasures(pos=dirty_stat.maxpos,direction=direc_max)

   Print out position of max in sensible angle units

print dm.time(direc_max.m0),dm.angle(direc_max.m1),' at ',dirty_stat.maxpos

   Output: 16:09:13.953 +065.32.28.998  at  [156 172 1 1]

Note that doing the last bit of stuff took several hours, as I had to hop
back and forth between the IMAGER, MEASURES, and QUANTA documentation.
And the fact that you can use QUANTITIES angle and time as MEASURES funtions
seems obscure (it isnt listed in measures for example) but I found it in the
example for measures.measure of all things!  Seems useful though...

>>>QUERY: IS THERE A MORE DIRECT WAY TO DO THIS?

   Done with image

im.close()

   Now we want to move the phase center and make a high res image there

imgr.open('1608.ms')
imgr.setimage(cellx='0.05arcsec', celly='0.05arcsec', nx=512, ny=512, 
  stokes='I', spwid=1, fieldid=1, phasecenter=direc_max)	
imgr.image(type='observed',image='1608.50mas.dirty')

   Do the stat thing again

im := image('1608.50mas.dirty')
im.statistics(statsout=dirty_stat,list=T) 

   Get nonsense.  Did it really shift?  I forgot the setdata crap...

imgr.close()
imgr.open('1608.ms')
imgr.setimage(cellx='0.05arcsec', celly='0.05arcsec', nx=512, ny=512, 
  stokes='I', spwid=1, fieldid=1, doshift=T, phasecenter=direc_max)	
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;

   Output
   Beam fit: 0.671522 by 0.621957 (arcsec) at pa -33.639 (deg) 
   This is the beam from the 0.2" resolution data.  I wonder if the
   uvrange is still in effect? Check summary...

imgr.summary()

   No mention of the uvrange.  Must still be in the IMAGING_WEIGHT
   column.  Will uvrange() default to all the data?

imgr.uvrange()

   Looks like it (Allowed uv range: 0 to 1e+10 wavelengths)
   OK, try again (after deleting old 1608.50mas.psf and .dirty)

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) 

   BAH! Its still wrong.  Try killing the imager and restarting

imgr.close()
imgr.done()
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_max)	
imgr.uvrange()
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;

   STILL NO GOOD.  The 1608.ms must be corrupted.  Redo from FITS
   again after deleting old 1608.ms.

imgr.close()
imgr.done()
m:=fitstoms(msfile='1608.ms',fitsfile='F16156554.UVF');
m.close() 
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_max)	
imgr.uvrange()
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;

   Output: Beam fit: 0.320926 by 0.228607 (arcsec) at pa -961.384 (deg)
   THIS WORKED. NOTE THE BPA IS EVEN WEIRDER THOUGH!

>>> QUERY: IS WHAT HAPPENED EXPECTED?  WAS THE 1608.ms CORRUPTED?
           I DONT KNOW IF THIS IS SOMETHING TO BUG() OR JUST ANOTHER
           USER ERROR...
  
   Continuing to make dirty map, get statistics,

imgr.image(type='observed',image='1608.50mas.dirty')
imgr.close()
im := image('1608.50mas.dirty')
im.statistics(statsout=dirty_stat,list=T)
print 'Max value = ',dirty_stat.max,' at ',dirty_stat.maxpos

   Output: Max value =  0.019444  at  [38 326 1 1]

   If the shift worked, I would have expected the max to now be at
  the center [255 256 1 1] or thereabouts!  Something is wrong.

im.coordmeasures(pos=dirty_stat.maxpos,direction=direc_max)
print 'Peak at ',dm.time(direc_max.m0),dm.angle(direc_max.m1),' at ',dirty_stat.maxpos

   Output: Peak at  16:09:12.464 +065.32.49.449  at  [38 326 1 1]

   Just what is the reference pixel in this map?

ref_pix := im.shape()/2; ref_pix[3:4] := [1,1];
im.coordmeasures(pos=ref_pix,direction=direc_ref)
print 'Ref at ',dm.time(direc_ref.m0),dm.angle(direc_ref.m1),' at ',ref_pix

   Output: Ref at  16:09:10.708 +065.32.45.950  at  [256 256 1 1]

   Compared to peak in low-resolution map:
                   16:09:13.953 +065.32.28.998  at  [156 172 1 1]

   IT IS CLEAR THAT JUST SETTING THE PHASE CENTER OF THE IMAGE DOES NOT
   SHIFT THE IMAGE.

>>>QUERY: IS IT SUPPOSED TO?  RECIPE 6 SUGGESTS THAT IT DOES

   Try just shifting (using shifts in imagr.setimage)
   Go back to the low-resolution dirty map (which still exists)

im.open('1608.200mas.dirty')
im.statistics(statsout=dirty_stat,list=T)
dirty_200mas_max := dirty_stat.maxpos; 
dirty_200mas_ref := im.shape()/2; dirty_200mas_ref[3:4] := [1,1];
dirty_200mas_shift := dirty_200mas_max - dirty_200mas_ref
print 'Shift in 200mas map is ',dirty_200mas_shift

   Output: Shift in 200mas map is  [-100 -84 0 0]
   Or, even better, in angular units

im.coordmeasures(pos=dirty_200mas_max,direction=direc_200mas_max)
im.coordmeasures(pos=dirty_200mas_ref,direction=direc_200mas_ref)
direc_200mas_shift := direc_200mas_max - direc_200mas_ref

   This fails - evidently the system isnt smart enough to know how
   to subtract measures.  (Note that subtracting the quanta .m0 and
   .m1 doesnt even work.  They should.)  Do manually...

shift_200mas_ra := direc_200mas_max.m0; shift_200mas_dec := direc_200mas_max.m1
shift_200mas_ra.value := direc_200mas_max.m0.value - direc_200mas_ref.m0.value
shift_200mas_dec.value := direc_200mas_max.m1.value - direc_200mas_ref.m1.value
print 'Peak in 200mas map is ',dm.time(direc_200mas_max.m0),dm.angle(direc_200mas_max.m1),' at ',dirty_200mas_max
print 'Center of 200mas map is ',dm.time(direc_200mas_ref.m0),dm.angle(direc_200mas_ref.m1),' at ',dirty_200mas_ref
print 'Shift in 200mas map is ',dm.time(shift_200mas_ra),dm.angle(shift_200mas_dec)

   Output: Peak in 200mas map is 16:09:13.953 +065.32.28.998 at [156 172 1 1]
   Output: Center of 200mas map is 16:09:10.732 +065.32.45.800 at [256 256 1 1]
   Output: Shift in 200mas map is  00:00:03.220 -000.00.16.802

   What is the bloody sense of this?  The OFFSET of the peak from the
   center seems to be W and S, so the shift should be reversed E and N

shift_200mas_ra.value := direc_200mas_ref.m0.value - direc_200mas_max.m0.value
shift_200mas_dec.value := direc_200mas_ref.m1.value - direc_200mas_max.m1.value
print 'Shift in 200mas map is ',dm.time(shift_200mas_ra),dm.angle(shift_200mas_dec)

   Output: Shift in 200mas map is  23:59:56.780 +000.00.16.802

print shift_200mas_ra, shift_200mas_dec

   Output: [unit=rad, value=-0.00023419] [value=8.1459e-05, unit=rad]

   Looks OK.  Back to imagr.

im.close()
imgr.open('1608.ms')

   Get default phase center from advise()

imgr.advise(takeadvice=T,phasecenter=phase_center)
print 'Using phase center ',dm.time(phase_center.m0),dm.angle(phase_center.m1)

   Output: Using phase center  16:09:10.700 +065.32.46.000
   This is the same as the center of the 200mas map 16:09:10.732 +065.32.45.800
   just quantized, so OK.
   But imagr uses strings shiftx,shifty, so we will have to visit
   quanta again :-(

shift_dec := dq.tos( dq.convert(shift_200mas_dec, 'arcsec') )
shift_ra := dq.tos( dq.convert(shift_200mas_ra, 'arcsec') )
print 'Shift x,y = ',shift_ra, shift_dec

   Output: Shift x,y =  -48.305 arcsec 16.802 arcsec

   To compare to difmap:

shift_rasec := dq.tos( dq.convert( dq.mul( shift_200mas_ra,dq.cos(phase_center.m1)), 'arcsec') )
print 'Shift x,y = ',shift_rasec, shift_dec

   Output: Shift x,y =  -19.996 arcsec 16.802 arcsec

   This compares well with difmap shift -20.2 (arcsec), 17 (arcsec).

   NOTE: I COULD HAVE USED THE DQ.SUB TO CALCULATE THE SHIFT ABOVE

shift_200mas_ra := dq.sub( direc_200mas_ref.m0, direc_200mas_max.m0 )
shift_200mas_dec := dq.sub( direc_200mas_ref.m1, direc_200mas_max.m1 )

   Anyway, continue...

imgr.setimage(cellx='0.05arcsec', celly='0.05arcsec', nx=512, ny=512, 
  stokes='I', spwid=1, fieldid=1, doshift=T, phasecenter=phase_center,
  shiftx=shift_ra, shifty=shift_dec)	
imgr.uvrange()
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) 

  Output: Beam fit: 0.321 by 0.228605 (arcsec) at pa -961.359 (deg)

  Approximately same as before (why is it slightly different?)

imgr.image(type='observed',image='1608.50mas.dirty')
imgr.close()
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 'Ref 50mas at ',dm.time(direc_50mas_ref.m0),dm.angle(direc_50mas_ref.m1)

  Output: Ref 50mas at  16:09:10.708 +065.32.45.950

  cf.     Peak in 200mas map is 16:09:13.953 +065.32.28.998 at [156 172 1 1]
        Center of 200mas map is 16:09:10.732 +065.32.45.800 at [256 256 1 1]

  This clearly hasnt worked, try going back to remake measurement set
  again.

m:=fitstoms(msfile='1608.ms',fitsfile='F16156554.UVF');
m.close() 
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=phase_center,
  shiftx=shift_ra, shifty=shift_dec)	
imgr.uvrange()
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) 
imgr.image(type='observed',image='1608.50mas.dirty')
imgr.close()
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 'Ref 50mas at ',dm.time(direc_50mas_ref.m0),dm.angle(direc_50mas_ref.m1)

   Output: Ref 50mas at  16:09:10.708 +065.32.45.950

>>>QUERY: DOES THE SHIFT IN SETIMAGE WORK?  

   It is late at night now.  Try again tomorrow.

   NOTE: This took all day.  There is a steep learning curve...

>>>QUERY: Below is the script I made from this session.  This is what
          I think should work:

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

include 'synthesis.g'
dowait:=T	# Tell the Object system to wait for a method to finish
imgr:=imager('1608.ms')	# make imager object
imgr.setimage(cellx='0.2arcsec', celly='0.2arcsec', nx=512, ny=512, 
  stokes='I', spwid=1, fieldid=1)	# Set the image properties. 
imgr.setdata(spwid=1, mode="none")	# Select the data 
imgr.uvrange(uvmin=0.0, uvmax=258000.0)	# Explicit uvrange for this cell/nx
bmaj:=F; bmin:=F; bpa:=F;	# Set up some return variables
imgr.image(type='psf',image='1608.200mas.psf')	# Make the PSF
imgr.fitpsf(psf='1608.200mas.psf', bmaj=bmaj, bmin=bmin, bpa=bpa) # Fit beam
system.print.precision := 5	# set global print precision to 5 sig digits
print 'Bmaj = ', bmaj;	# Print out the beam size
print 'Bmin = ', bmin;  print 'Bpa = ', bpa;
imgr.image(type='observed',image='1608.200mas.dirty')	# Make the dirty image
# Done with imagr for now
imgr.close()
# Use image constructor to access dirty image
im := image('1608.200mas.dirty')
im.statistics(statsout=dirty_200mas_stat,list=T) # get whole image stats
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
im.coordmeasures(pos=dirty_200mas_max,direction=direc_200mas_max)
print 'Peak at ',dm.time(direc_200mas_max.m0),dm.angle(direc_200mas_max.m1),' at ',dirty_200mas_max
# Get reference pixel in map
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 at ',dm.time(direc_200mas_ref.m0),dm.angle(direc_200mas_ref.m1),' at ',dirty_200mas_ref
# Get shift
shift_200mas_ra := dq.sub( direc_200mas_ref.m0, direc_200mas_max.m0 )
shift_200mas_dec := dq.sub( direc_200mas_ref.m1, direc_200mas_max.m1 )
shift_dec := dq.tos( dq.convert(shift_200mas_dec, 'arcsec') )
shift_ra := dq.tos( dq.convert(shift_200mas_ra, 'arcsec') )
print 'Shift x,y in 200mas map = ',shift_ra, shift_dec
im.close()
#
# Back to imager, but shift phasecenter and remake dirty image
#
imgr.open('1608.ms')
imgr.advise(takeadvice=T,phasecenter=phase_center)
imgr.setimage(cellx='0.05arcsec', celly='0.05arcsec', nx=512, ny=512, 
  stokes='I', spwid=1, fieldid=1, doshift=T, phasecenter=phase_center,
  shiftx=shift_ra, shifty=shift_dec)	
imgr.uvrange()
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')
imgr.close()
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 'Ref 50mas 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