#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 low-res dirty image
#
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
imgr.weight('robust')
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 'Bmin = ', bmin;  print 'Bpa = ', bpa;
imgr.image(type='observed',image='1608.200mas.dirty')	# Make the dirty image
#
# 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) # 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
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