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

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

#
# As in the first example, we will smooth the data before imaging,
# combining groups of 100 channels into a single channel. This
# dramatically speeds up processing.
#

avg_data_1 = "TWHYA_BAND7_AVG.ms"
os.system("rm -rf "+avg_data_1)
split(vis=data_1,
      outputvis=avg_data_1,
      datacolumn='data',
      width=100,
      spw='')

avg_data_2 = "TWHYA_BAND6_AVG.ms"
os.system("rm -rf "+avg_data_2)
split(vis=data_2,
      outputvis=avg_data_2,
      datacolumn='data',
      width=100,
      spw='')

avg_data_3 = "TWHYA_BAND3_AVG.ms"
os.system("rm -rf "+avg_data_3)
split(vis=data_3,
      outputvis=avg_data_3,
      datacolumn='data',
      width=100,
      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 avg_data_1 etc. are applied.
#

avg_data_1 = "TWHYA_BAND7_AVG.ms"
avg_data_2 = "TWHYA_BAND6_AVG.ms"
avg_data_3 = "TWHYA_BAND3_AVG.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 be 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 = [avg_data_1, avg_data_2, avg_data_3]
vis = [avg_data_2, avg_data_3]
#vis = avg_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

# 
# After you have typed "go" the viewer will pop up and you will
# proceed to clean interactively. Review the suggestions regarding
# clean boxing, proceeding, and stopping under the first script on the
# Workshop casaguide
#

# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%
# 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')