# 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'

# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%
# (OPTIONAL) INSPECT THE INPUT DATA
# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%

#
# List the basic parameters of the input data set.
#

listobs(vis=data)

# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%
# SMOOTH THE DATA
# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%

#
# Smooth the data in frequency.
#

# Define the output name for the smoothed data:
smooth_data = 'Antennae_North.cal.ms'

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

# ... first removing any previous version
os.system("rm -rf "+smooth_data)

width = 10
split(vis=data,
      outputvis=smooth_data,
      datacolumn='data',
      width=width,
      spw='')
#clearcal(smooth_data)

#os.system("cp -r "+data+" .")

# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%
# (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=smooth_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=smooth_data
#vis="Antennae_North.cal.ms"
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'
#spw='0:1~30;120~164'

# for continuum we will use the multifrequency synthesis mode, here
# pick nterms=1 but we will move on to nterms > 1 in the next example.
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.5mJy'
# thresh = '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'
imagermode = 'csclean'
phasecenter = '12'

# set the u-v weighting scheme. This will determine how the u-v data
# are weighted, affecting the resultant PSF.

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

# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%
# 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.
#