# CASA can combine multiple pointing centers into a single image in a
# procedure, creating a "mosaic" image. Here we see the basic syntax
# to do this. We create a continuum image of the northern part of the
# Antennae galaxies.

# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%
# SET INPUTS AND OUTPUTS
# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%

#
# First define variables to hold the name of the input data set and
# the root name of the output data. 
#

data = "../../calib/Antennae_Band7_CalibratedData/Antennae_North.cal.ms"
image_out = 'ANTENNAE_NORTH_BAND7_CONT'

# List the basic parameters of the input data set. Have a particular
# look at the fields.
listobs(vis=data)

# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%
# AVERAGE THE DATA IN FREQUENCY
# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%

# Define the output name for the frequency-averaged data. For reasons
# that we don't want to talk about it's very important that you not
# change this.
avg_data = 'Antennae_North.cal.ms'

# Use the SPLIT task to average the data in velocity.

# ... first removing any previous version
os.system("rm -rf "+avg_data)
width = 10
split(vis=data,
      outputvis=avg_data,
      datacolumn='data',
      width=width,
      spw='')

# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%
# (OPTIONAL) PLOT U-V COVERAGE
# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%

#
# Plot the u-v coverage of your data. To see the many distinct
# pointings, go to "display" and colorize by "field".
#

plotms(vis=avg_data, 
       xaxis='u', 
       yaxis='v',
       avgchannel='10000',
       avgspw=False,
       avgtime='1e9',
       avgscan=False)       

# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%
# SET TUNING PARAMETERS FOR CLEAN
# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%

#
# Set the parameters that CLEAN will use to grid and deconvolve the
# u-v data into your final image. 
#

# start by setting the task to clean and the inputs to their defaults
default('clean')

# image the frequency-smoothed data
vis=avg_data
imagename=image_out

# we want to do this interactively
interactive=True

# combine all frequencies where the CO line isn't obviously bright
spw='0:0~6;10~16'

# use the multifrequency synthesis (but not spectral index) mode
mode = 'mfs'
nterms = 1

# set the conditions that will cause CLEAN to stop. Either the maximum
# iterations or the threshold value of the peak residual. Neither of
# these conditions is likely to be met, so this effectively tells
# CLEAN that we will stop it interactively. After you have an idea of
# the RMS in your image, you can come back and set 'thresh' to a
# reasonable value (say ~ 2-3 sigma).
niter = 1000
threshold = '0.1mJy'

# ... CELL SIZE (WANT ~3-5 CELLS PER SYNTHESIZED BEAM)
cell = '0.15arcsec'

# ... IMAGE SIZE (TRY TO IMAGE OUT TO ~ PRIMARY BEAM BY DEFAULT)
imsize = 500

# ***Tell CLEAN to image the data as a mosaic***
imagermode = 'mosaic'
phasecenter = '12'

# set the u-v weighting scheme.
# weighting = 'natural'
weighting = 'briggs'
robust = 0.5
# robust = 0.0

# erase previous images (skip this step if you want to continue cleaning)
os.system('rm -rf '+image_out+'.*')

# review inputs
inp

# execute CLEAN
go

#
# This will send you to interactive cleaning. Review the notes in the
# CASAguide for the first script for details of interactive cleaning.

# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%
# VIEW
# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%

#
# Display the resulting image. Note that if you are using nterms > 1
# this image is instead called .image.tt0
#

viewer(image_out+'.image')

viewer(image_out+'.flux')

#
# You may also want to inspect the residual, model, mask, flux
# (primary beam response), and (if using nterms > 1) alpha images.
#