# Use this script to explore using CASA to self-calibrate ALMA Science
# Verification data. The basic procedure for self-calibration,
# outlined in the previous talk, is to build a model of your source
# and then use "gaincal" to solve for antenna dependent complex gain
# terms that best match the data to that model. This can be an
# iterative process. Here we'll step through a first iteration for
# either NGC 3256 or TW Hydra. See the casaguides for more detailed
# discussion.

# Note that this script can equally well apply to either TW Hydra or
# NGC 3256. Just juggle the comments (NGC 3256 is "Band 3", TW Hydra
# is "Band 7") as you go through the script. The hosting CASAguide
# contains a few helpful comments, review those as you got through the
# script.

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

#
# We will use two data sets here, NGC 3256 and TW Hydra. Feel free to
# put in the other TW Hydra data if you like. 
#

# 
# We point at the calibrated (but not yet *self*-calibrated) data and
# define two output images: the first made before self calibration and
# the second made after self-calibration.
#

data = "../../calib/TWHYA_BAND7_CalibratedData/TWHydra_corrected.ms"
first_image = 'TWHYA_BAND7_CONT'
selfcal_image = 'TWHYA_BAND7_SELFCAL'

#data = "../../calib/NGC3256_Band3_CalibratedData/ngc3256_line_target.ms"
#first_image = 'NGC3256_BAND3_CONT'
#selfcal_image = 'NGC3256_BAND3_SELFCAL'

# Get a basic description of the input data
listobs(vis=data)

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

# Define the output name for the averaged data:
avg_data = first_image+'_AVG.ms'

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

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

width = 100 # for TW Hydra
# width = 4 # for NGC 3256

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

# Flag the line channels after averaging. Note that we skipped this
# step before (in "basic imaging"). You could identify the line
# channels either from an integrated spectrum we used to prepare the
# continuum subtraction or the line cubes that we made (see the "line
# imaging" exercise).

# Line channels for TW Hydra.
linechans = ['0:16~16','2:21~21','3:33~37']

# Line channels for NGC 3256
linechans = ['0:15~16','1:0~16']
flagdata(vis=data, spw=linechans)

# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%
# CARRY OUT A FIRST IMAGING
# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%

# 
# This step almost exactly parallels the basic continuum imaging
# step. The most important thing here is that you make sure that
# calready is set to "True" so that the deconvolved model is read into
# the data set for future use with gaincal. Otherwise just mask and
# image as you usually would, stopping when the source is
# indistinguishable from the surrounding residuals.
#

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

# This is the *the* key to self-calibration. Setting calready=True
# tells CLEAN to place the CLEANed model in the "model" column of the
# data set. This will allow "gaincal" to operate on the source data.
calready=True

# image the frequency-smoothed data
vis=avg_data
imagename=first_image

# we want to do this interactively
interactive=True

# SPW selections to avoid bright lines
spw = ''

# use the multifrequency synthesis mode with nterms=1
mode = 'mfs'
nterms = 1

# for this threshold you will need to stop cleaning manually
niter = 10000
threshold = '0.1mJy'

# ... CELL SIZE
# Good for Bands 6 & 7
# cell = '0.3arcsec'
# cell = '0.5arcsec'
# Good for Band 3
cell = '1.0arcsec'

# ... IMAGE SIZE
# Good for Bands 6 & 7
# imsize = 100
# Good for Band 3
imsize = 300

imagermode = 'csclean'

weighting = 'briggs'
robust = 0.5

# erase previous images (skip this step if you want to continue cleaning)
os.system('rm -rf '+first_image+'.*')

# review inputs
inp

# execute CLEAN
go

# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%
# INSPECT THE OUTPUT
# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%

# 
# So far this is familiar from this morning's tutorial. Have a look at
# the image that you have made and be sure to note the noise in the
# image. You can measure this from a signal-free region by dragging
# out the region and double-clicking. Look for the output in the
# terminal.
#

viewer(first_image+'.image')

# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%
# CARRY OUT A SELF-CALIBRATION
# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%

# 
# CLEAN has placed the model of the source into the data (if you want
# to see this for yourself use either 'plotms' and select to visualize
# the "model" column or use 'browsetable' and simply look at that
# column in your data in grid form). 
#
# Now that you have a model of the source, you can use the 'gaincal'
# command to solve for the antenna-based phase or amplitude terms
# needed to correct the data to match the model.
#
# To do this, we'll specify a reference antenna (here just picked from
# the CASA guides). We'll also need to pick a time interval over which
# we average to derive solutions. The shorter your time interval, the
# more you are able to capture rapid variations in phase and amplitude
# response of the telescopes or atmosphere. The time cannot be
# arbitrarily short because you need to achieve enough signal-to-noise
# within each solution interval to constrain the complex gain
# solutions. Therefore the solution interval, specified via the
# "solint" keyword, is a key knob in making selfcal work. You want it
# to be as short as possible while still obtaining good solutions.
#
# The CASA guides go over the considerations for this in some
# detail. In practice you can also experiment with different solints
# to see what works.
#
# Note that we will only run a phase based selfcal in this first
# call. After you are satisfied that the phase selfcal works, you can
# change the calmode to 'a' or 'ap' to also solve for amplitude
# corrections.

# for NGC 3256
#refant = 'DV07'
#solint = '1800s'
#caltable = 'selfcal_ngc3256.gcal'

# for TW Hydra
refant='DV06'
solint='30s'
caltable = 'selfcal_twhy_band7.gcal'

gaincal(vis=avg_data,
        field='',
        caltable=caltable,
        spw='',
        solint=solint,
        refant=refant,
        calmode='p',
        minblperant=4)

#
# Watch out for failed solutions noted in the terminal during this
# solution. If you see a large fraction (really more than 1 or 2) of
# your antennas failing to converge in many time intervals then you
# may need to lengthen the solution interval.
#

# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%
# INSPECT THE CALIBRATION
# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%

#
# After you have run the gaincal, you want to inspect the
# solution. Use PLOTCAL to look at the solution (here broken into
# panels by SPW with individual antennas mapped to colors). Look at
# the overall magnitude of the correction to get an idea of how
# important the selfcal is and at how quickly it changes with time to
# get an idea of how stable the instrument and atmosphere were.
#

plotcal(caltable=caltable,
        xaxis='time', yaxis='phase',
        plotrange=[0,0,-180,180],
        iteration='spw', subplot = 221)

# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%
# APPLY THE CALIBRATION
# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%

#
# If you are satisfied with your solution, you can now apply it to the
# data to generate a new corrected data column, which you can then
# image. Be sure to save the previous flags before you do so because
# applycal will flag data without good solutions. The commented
# command after the applycal will roll back to the saved solution in
# case you get in trouble.
#

flagmanager(vis=avg_data,
            mode='save',
            versionname='before_selfcal_apply')
applycal(vis=avg_data,
         gaintable=caltable,
         interp='linear',
         flagbackup=False)
#flagmanager(vis=avg_data,
#            mode='restore',
#            versionname='before_selfcal_apply')

# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%
# RE-IMAGE THE SELF-CALIBRATED DATA
# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%

# 
# Now that you have applied your self-calibration to the data you are
# ready to build an improved image. Repeat your CLEAN call but now
# directing the output to a new image so that you can compare the two
# cases. Note that CLEAN automatically images the corrected data
# column if that column is present, so that it will focus on the
# results of applycal automatically.
#


# start by setting the task to clean and the inputs to their defaults
default('clean')
calready=True
vis=avg_data
imagename=selfcal_image
interactive=True
swp = ''
mode = 'mfs'
nterms = 1
niter = 10000
threshold = '0.1mJy'
# cell = '1.0arcsec' # for NGC 3256
cell = '0.3arcsec' # for TW Hydra
imsize = 300
imagermode = 'csclean'
weighting = 'briggs'
robust = 0.5
os.system('rm -rf '+selfcal_image+'.*')
inp
go

# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%
# INSPECT THE NEW IMAGE
# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%

#
# Inspect the output image, noting the noise. Compare it to the
# previous, not self-calibrated image.
#

viewer(selfcal_image+'.image')