cd /home/scidata/cdebreuc/alma/TST2.0/N4826
aips++


#
# Feathering
#
img:=imagefromfits('n4826_12mchan.im','NGC4826.12motf.chan.fits');
img.srb('55arcsec','55arcsec','0deg');
img.sbu('Jy/beam');
img.done();
img:=imager('n4826_both.ms');
img.setvp(dovp=T);
img.feather(image='NGC4826.feathered',highres='n4826_mom0.im',lowres='n4826_12mmom0.im');
img.done();
# -----------------------------------------------------------------------------

#
# uv-plane combination directly converting SD image to a model
# ============================================================
#
# Start with interactively defining a mask
#
include 'interactivemask.g' 
imask:=interactivemask(refimage='n4826_mom0.im',mask='n4826.mask'); # Use BIMA moment 0 map to define mask
imask.start(); # define mask using polygon
# the mask is ready, but I had problems when using this for the first time. Mayby forgot to close imask tool?
# SEVERE ERROR: Table n4826.mask cannot be created because it is in use in another process. Don't know how to close this process -> exit aips++ and restart seems to work.
#
# Start making a model directly from the SD image
#
im:=imager(filename='n4826_both.ms'); # create imager tool
im.setdata(fieldid=[1:7],spwid=[1:3],nchan=[64,64,64],start=[1,1,1],step=[1,1,1]); # select all channels in the first 3 spwids & all fields in the mosaic
im.setimage(nx=256,ny=256,cellx='1arcsec',celly='1arcsec',stokes='I',mode='channel',nchan=30,start=47,step=4,fieldid=1,spwid=[1:3]); # set total mosaic size, start with channel 47, make 30 total channels with 4 channel averaging, fieldid=1 is the central field of the mosaic
im.setmfcontrol(scaletype="NONE",minpb=0.1); # no scaling, 10% PB correction
im.setvp(dovp=T); # Do primary beam correction
im.makemodelfromsd(sdimage='n4826_12mchan.im',modelimage='n4826_joint1',maskimage='n4826.mask');
im.clean(algorithm='mfclark',model='n4826_joint1',gain=0.2,niter=500,mask='n4826.mask'); # clean the data cube using SD as model
im.done();
img:=image('n4826_joint1.restored');
img.stats();
#
# Cleaned cube is ready, now make a moment 0 map of it
#
img_mom0:=img.moments(outfile='n4826_joint1.mom0',moments=0,axis=4,mask='indexin(4,[5:26])',includepix=[0.070,1000.0]); # create a moment 0 map
img_mom0.stats();
img_mom0.done();
img.done();
# -----------------------------------------------------------------------------


#
# uv-plane combination using SD image deconvolved with multi-scale CLEAN
# ======================================================================
#
# Regrid the SD image cube to have the same coordinates as the sysntesis data
#
ima:=image('n4826_joint1.restored'); # Use previous image as grid example
csys:=ima.coordsys(); # retrieve details of coordinate system
img:=image(infile='n4826_12mchan.im'); # SD image
imb:=img.regrid(outfile='n4826.12motf.chregrid.im',shape=[256,256,1,30],csys=csys); # Regrid SD cube to synthesis mosaic coordinate system
ima.done();
imb.done();
img.done();
#
# Deconvolve regridded SD image to correct for response of SD beam
#
include 'deconvolver.g'
decon:=deconvolver(dirtyname='n4826.12motf.chregrid.im',psfname=''); # create deconvolver tool using regridded SD image
decon.makegaussian(gaussianimage='n12m_gaussian.im',bmaj='77.1arcsec',bmin='77.1arcsec',bpa='134deg',normalize=F); # make gaussian PSF
decon.open(dirtyname='n4826.12motf.chregrid.im',psfname='n12m_gaussian.im'); # use PSF generated above
decon.setscales(scalemethod='uservector',uservector=[30,60]); # set scales for multi-scale CLEAN
decon.clean(algorithm='msclean',model='n4826_joint2',niter=100,gain=0.3); # multi-scale CLEAN
decon.done(); # close tool and start joint deconvolution
#
# Start joint deconvolution 
#
im:=imager(filename='n4826_both.ms'); # create imager tool
im.setdata(fieldid=[1:7],spwid=[1:3],nchan=[64,64,64],start=[1,1,1],step=[1,1,1]); # select all channels in the first 3 spwids & all fields in the mosaic
im.setimage(nx=256,ny=256,cellx='1arcsec',celly='1arcsec',stokes='I',mode='channel',nchan=30,start=47,step=4,fieldid=1,spwid=[1:3]); # set total mosaic size, start with channel 47, make 30 total channels with 4 channel averaging, fieldid=1 is the central field of the mosaic
im.setmfcontrol(scaletype="NONE",minpb=0.1); # no scaling, 10% PB correction
im.setvp(dovp=T); # Do primary beam correction
im.clean(algorithm='mfclark',model='n4826_joint2',gain=0.2,niter=500,mask='n4826.mask'); # clean the data cube using SD as model
im.done();
img:=image('n4826_joint2.restored');
img.stats();
img.done();
#
# Cleaned cube is ready, now make a moment 0 map of it
#
im:=image('n4826_joint2.restored');
im_mom0:=im.moments(outfile='n4826_joint2.mom0',moments=0,axis=4,mask='indexin(4,[5:26])',includepix=[0.070,1000.0]); # create a moment 0 map
im_mom0.stats();
im_mom0.done();
# -----------------------------------------------------------------------------


#
# Calculate statistics on SD image
#
# Need to first construct a moment 0 map
img:=image('n4826.12motf.chregrid.im');
img_mom0:=img.moments(outfile='n4826_12motf.chregrid.mom0',moments=0,axis=4,mask='indexin(4,[5:26])',includepix=[0.070,1000.0]);
img_mom0.srb('55arcsec','55arcsec','0deg');
img_mom0.sbu('Jy/beam'); # have to cheat here, units are Jy/beam.km/s
img_mom0.stats(); # only gives total flux when units are Jy/beam
img_mom0.done();