# The high signal-to-noise and wide bandwidth in ALMA data can lead to
# the case where the approximation of a single flux at all frequencies
# is an indequate model for continuum imaging. This script explores
# CASA's multi-frequency synthesis capabilities using the three TW
# Hydra data sets.

# A caveat before we begin: the data need to be pretty carefully
# calibrated for this to work well (just like for any ratio) and this
# script deals with the data before any self-calibration has been
# applied. The approach here should be correct but the results should
# not be overemphasized.

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

#
# These are the three TW Hydra Science verification sets. We have data
# at Band 3 (~100 GHz), 6 (~230 GHz), and 7 (~350 GHz). Use 'listobs'
# to get the details of each. We can image these with a spectral term
# either separately or together.
#

data_1 = "../../calib/TWHYA_BAND7_CalibratedData/TWHydra_corrected.ms"
data_2 = "../../calib/TWHYA_BAND6_CalibratedData/TWHydra_corrected.ms"
data_3 = "../../calib/TWHYA_BAND3_CalibratedData/TWHydra_corrected.ms"

# listobs(vis=data_1)
# listobs(vis=data_2)
# listobs(vis=data_3)

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

#
# As in the first example, we will smooth the data before imaging,
# combining groups of 10 channels into a single channel. This
# dramatically speeds up processing (note that we already smoothed a
# bit before providing you the data, too). 
#
# (Don't let the warnings about incomplete edge channels worry you -
#  CASA is just complaining that it doesn't find a full 10 channels in
#  the last batch.)
#

smooth_data_1 = "TWHYA_BAND7_SMOOTH.ms"
os.system("rm -rf "+smooth_data_1)
split(vis=data_1,
      outputvis=smooth_data_1,
      datacolumn='data',
      width=10,
      spw='')

smooth_data_2 = "TWHYA_BAND6_SMOOTH.ms"
os.system("rm -rf "+smooth_data_2)
split(vis=data_2,
      outputvis=smooth_data_2,
      datacolumn='data',
      width=10,
      spw='')

smooth_data_3 = "TWHYA_BAND3_SMOOTH.ms"
os.system("rm -rf "+smooth_data_3)
split(vis=data_3,
      outputvis=smooth_data_3,
      datacolumn='data',
      width=10,
      spw='')

#
# If you have already smoothed the data and are starting a new session
# just paste these three lines in to the session to ensure that the
# names smooth_data_1 etc. are applied.
#

smooth_data_1 = "TWHYA_BAND7_SMOOTH.ms"
smooth_data_2 = "TWHYA_BAND6_SMOOTH.ms"
smooth_data_3 = "TWHYA_BAND3_SMOOTH.ms"

# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%
# SET PARAMETERS FOR THE CLEAN
# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%

#
# Now set the parameters for CLEAN. The main difference from before is
# that we will use nterms = 2 (or you can try higher orders) and that
# we will combining multiple data sets taken at different bands. We're
# working under the assumption that the emission at these bands is the
# same related only by a simple spectral behavior. For the continuum
# emission at Bands 6 and 7 (and probably Band 3) this is almost
# certainly a reasonable assumption.
#

default('clean')

# set the set of data to use in the CLEAN. Note that CLEAN will accept
# a list of data sets in and image them all together.

#vis = [smooth_data_1, smooth_data_2, smooth_data_3]
vis = [smooth_data_2, smooth_data_3]
#vis = smooth_data_3

# set the output name
#imagename = 'TWHYA_MFS_B367'
imagename = 'TWHYA_MFS_B67'
#imagename = 'TWHYA_MFS_B7'

# set interactive to true
interactive = True

# set the mode to multifrequency synthesis but now set nterms to 2
# instead of 1. The result will be CLEAN components that have both a
# frequency behavior (spectral index) and a flux.
mode = 'mfs'
nterms = 2

# set the stopping conditions
threshold = '1.0mJy'
# niter = 10000
niter = 100

# set the cell size, aiming for 3-5 cells per synthesized beam
cell = '0.3arcsec'
#cell = '0.5arcsec'

# ... set the image size
imsize = 300

imagermode = 'csclean'

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

# erase previous output images
os.system('rm -rf '+imagename+'*')

# review inputs
inp

# start clean
go

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

# Note that the output of the MFS clean with nterms 2 or more is a bit
# different than the original CLEAN output. The intensity image at a
# fiducial frequency near the center of the your range has extension
# .tt0. There is also a spectral index image with extension .alpha.

viewer(imagename+'.image.tt0')
viewer(imagename+'.image.alpha')