# Use this script to see the workflow for applying self-calibration
# derived from continuum emission to improve line imaging. As in the
# other scripts, we use ALMA Science Verification data. This script
# builds on the imaging, line imaging, and self-calibration
# scripts. See those for more details.

# 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 not-self-calibrated continuum and the
# self-calibrated line image.
#

#data = "../../calib/TWHYA_BAND7_CalibratedData/TWHydra_corrected.ms"
#cont_image = 'TWHYA_BAND7_CONT'
#line_image = 'TWHYA_BAND7_SELFCAL_LINE'

data = "../../calib/NGC3256_Band3_CalibratedData/ngc3256_line_target.ms"
cont_image = 'NGC3256_BAND3_CONT'
line_image = 'NGC3256_BAND3_SELFCAL_LINE'

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

# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%
# FOCUS ONLY ON THE CO DATA
# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%

# 
# To keep thing simple, we will split out only the spectral window
# containing CO emission. The continuum from the other spectral
# windows may be helpful to improve the signal-to-noise (via
# "combine=spw" in gaincal), we only focus on one window for
# simplicity. Feel free to experiment with that by changing the
# parameters here.

# NGC 3256 CO
local = 'ngc3256_co.ms'
split(vis=data,
      outputvis=local,
      spw='0',
      datacolumn='data')

# TW HYDRA BAND 7 CO AND HCO+
#local = 'tw_hydra_band7_co.ms'
#split(vis=data,
#      outputvis=local,
#      spw='2',
#      datacolumn='data')

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

#
# We will still image the local data for the line imaging, but there's
# no reason not to average the data in frequency here to speed up
# building the continuum image. So here we use SPLIT to build an
# data set averaged in velocity.
#

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

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

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

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

# Flag the line channels after averaging.

# Line channels for TW Hydra.
#linechans = ['0:21~21']
#flagdata(vis=data, spw=linechans)

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

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

# 
# As in the basic selfcal example, image the continuum and read the
# results back into the model column.
#

# start by setting the task to clean and the inputs to their defaults
default('clean')
calready=True
vis=avg_data
imagename=cont_image
interactive=True
spw = ''
mode = 'mfs'
nterms = 1
niter = 10000
threshold = '0.1mJy'
# ... CELL SIZE
# Good for Bands 6 & 7
#cell = '0.3arcsec'
# Good for Band 3
cell = '1.0arcsec'
# ... IMAGE SIZE
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 '+cont_image+'.*')
# review inputs
inp
# execute CLEAN
go

# (OPTIONALLY) INSPECT THE OUTPUT
viewer(first_image+'.image')

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

# Self-calibrate using the model of the continuum data.

# 
# CLEAN has placed the model of the source into the data
# 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)

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

# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%
# APPLY THE CALIBRATION TO THE LINE DATA
# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%

# 
# Here is the trick to the line application. Apply the calibration
# table derived from the continuum data to the LINE data set (copied
# to the "local" file.
#

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

# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%
# CONTINUUM SUBTRACTION
# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%

#
# Now carry out a continuum subtraction as we did in the line imaging
# example this morning.
#

# ... for NGC 3256 CO
uvcontsub(vis = local,
           fitspw='0:20~53;71~120',
           solint ='int',
           fitorder = 1,
           combine='')

# ... for TW Hydra Band 7 CO
#uvcontsub(vis=local,
#          fitspw='0:20~2000,0:2400~3800',
#          solint='int',
#          fitorder=1,
#          combine='')

# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%
# IMAGE THE LINE DATA
# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%

# 
# Now image the line data. Use the parameters from this morning so
# that you can make a ready comparison.
#

default('clean')
# image the continuum-subtracted data
vis=local+'.contsub'
imagename=line_image
# stop manually
interactive=True
threshold = '0.1mJy'
niter = 10000
# ... CELL SIZE
#cell = '0.3arcsec' # for Bands 6 and 7
cell = '1.0arcsec' # for Band 3
# ... IMAGE SIZE
imsize = 300
imagermode = 'csclean'
psfmode='hogbom'
weighting = 'briggs'
robust = 0.5
# CHANNEL MODE ... for NGC 3256 CO
restreq = '115.271201800GHz'
mode='channel'
start=''
spw='0:38~87'
width=''
nchan=50
# VELOCITY MODE ... for TW Hydra Band 7 (either HCO+ or CO 3-2)
#restfreq='345.79599GHz'
#mode='velocity'
#start='-4km/s'
#width='0.12km/s'
#nchan=118
#outframe='LSRK'
# erase previous images (skip this step if you want to continue cleaning)
os.system('rm -rf '+imagename+'.image')
os.system('rm -rf '+imagename+'.model')
os.system('rm -rf '+imagename+'.psf')
os.system('rm -rf '+imagename+'.mask')
os.system('rm -rf '+imagename+'.residual')
os.system('rm -rf '+imagename+'.flux')
# Inspect input and go.
inp
go

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

viewer(line_image+'.image')