# Use this script to explore CASA imaging using ALMA Science
# Verification data.

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

#
# First define variables to hold the name of the input data set and
# the root name of the output data. These are read into variables that
# are then fed to parameters in the call to CLEAN below. Only the most
# recent definition of data and image_out will be used.
#

#data = "../../calib/TWHYA_BAND7_CalibratedData/TWHydra_corrected.ms"
#image_out = 'TWHYA_BAND7_CONT'

#data = "../../calib/TWHYA_BAND6_CalibratedData/TWHydra_corrected.ms"
#image_out = 'TWHYA_BAND6_CONT'

#data = "../../calib/TWHYA_BAND3_CalibratedData/TWHydra_corrected.ms"
#image_out = 'TWHYA_BAND3_CONT'

data = "../../calib/NGC3256_Band3_CalibratedData/ngc3256_line_target.ms"
image_out = 'NGC3256_BAND3_CONT'

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

#
# List the basic parameters of the input data set. Pay attention to
# the spectral setup here, as it will inform how (and if) we smooth
# the data in just a moment.
#

listobs(vis=data)

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

#
# Smooth the data in frequency. This will make the imaging and
# cleaning operation MUCH faster. You do this using SPLIT with the
# width parameter. For the TW Hydra data sets, smoothing together 100
# channels makes sense. The native channel width in the NGC 3256 data
# set is coarser, so we smooth only by a factor of 4.
#

# Define the output name for the smoothed data:
smooth_data = image_out+'_SMOOTH.ms'

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

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

# width = 100
width = 4
split(vis=data,
      outputvis=smooth_data,
      datacolumn='data',
      width=width,
      spw='')

#
# COMMENTARY: We gloss over a couple potentially important steps. We
# could omit noisy edge channels or flag line data that may
# contaminate the continuum map. See the CASA guide on TW Hydra for a
# good example of this. We're stepping past this step for
# expediency. If you wanted a science-quality image you should do
# this.
#
# In principle, you also need to be careful to avoid averaging
# together data that have distinct u-v positions due to their
# different frequencies. In practice the smoothing width that we use
# for data sets that we consider here will be fine.

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

#
# Plot the u-v coverage of your data. Note that you also get the
# complement of these.
#

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. After smoothing, CLEAN should run
# fairly quickly so you can experiment with these and get an idea of
# their impact on the final image.
#

#
# A NOTE ON SYNTAX: As you have already seen, CASA has two modes for
# calling tasks. As a general rule the "in-line" call, e.g., task(x=x,
# y=y), is more robust. However for exploration, the "tget" / "inp" /
# "go" syntax can be very useful and we suggest using that here.
#

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

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

# we want to do this interactively
interactive=True

# combine all frequencies
spw=''

# 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 = 10000
threshold = '0.1mJy'
# thresh = '1mJy'

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

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

imagermode = 'csclean'

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

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