# 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.
# Edit the commented lines to set the dataset you wish to work on. 

#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. You can check the spectral resolution using
# listobs.  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 with these data 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. This will only be necessary
# if you have run the smoothing step before and want to recreate the
# smooth data from the original.
os.system("rm -rf "+smooth_data)

# Edit the smoothing width for the dataset you are working on:
# width = 10 # for TW Hydra ("FDM" mode)
width = 4 # for NGC 3256 ("TDM" mode)

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

# COMMENTARY: We gloss over a couple potentially important steps for
# expediency. We could omit noisy edge channels or flag line data that
# may contaminate the continuum map. If you wanted a science-quality
# image you should do this. See the CASA guide on TW Hydra for a good
# example of how to perform these steps:
# http://casaguides.nrao.edu/index.php?title=TWHydraBand7_Imaging

# 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 the datasets 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=''

# To image the continuum we will use the multifrequency synthesis
# mode.  Here we set nterms=1, but will move on to nterms > 1 in the
# next example.
mode = 'mfs'
nterms = 1

# Next, we set the conditions that will cause CLEAN to stop.  CLEAN
# will stop when either the maximum iterations are run, or the value
# of the peak residual is less than the threshold value.  For the
# values chosen in this initial run-through, 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' # a condition that is likely to be met

# ... CELL SIZE (WANT ~3-5 CELLS PER SYNTHESIZED BEAM)
# Good for Bands 6 & 7
# cell = '0.3arcsec'
# cell = '0.5arcsec'
# Good for Band 3
cell = '1.0arcsec'

# ... IMAGE SIZE (TRY TO IMAGE OUT TO ~ PRIMARY BEAM BY DEFAULT)
# Good for Bands 6 & 7
# imsize = 100
# Good for Band 3
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'
# robust = 0.0
weighting = 'briggs'
robust = 0.5

# Erase previous images (skip this step if you want to continue
# cleaning from a previous run). Note that this will delete the .mask
# file as well. If you want to start over but keep your mask, you can
# rename this file and feed it in to the mask parameter for clean.
os.system('rm -rf '+image_out+'.*')

# review inputs
inp

# execute CLEAN
go

# If you are using the TW Hydra Band 7 dataset, for example, the
# interactive viewer will look like this:
# http://casaguides.nrao.edu/index.php?title=File:Viewer_clean1.png
# Try setting a CLEAN mask and interactively running CLEAN until the
# residuals are consistent with noise as described on the wiki. You
# can also edit the weighting scheme, etc. and re-run. Remember to
# remove previous images if you want to CLEAN from the beginning.

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