# SELF-CALIBRATION USING A SPECTRAL LINE

# This final script illustrates the application of self-calibration to
# a bright spectral line. We focus exclusively on the NGC 3256 data
# set and the CO line. Unlike in the previous line example (which
# built on the continum selfcal), we start from scratch.

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

# 
# 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/NGC3256_Band3_CalibratedData/ngc3256_line_target.ms"
first_image = 'NGC3256_BAND3_LINE'
selfcal_image = 'NGC3256_BAND3_SELFCAL_LINE'

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

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

# 
# Again we average the data into a single-channel data set but now we
# focus ONLY on the channels containing bright line emission.
#

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

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

# ... average together all channels across the bright part of the line
width = 5
split(vis=data,
      outputvis=avg_data,
      datacolumn='data',
      width=width,
      spw='0:60~64')

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

# start by setting the task to clean and the inputs to their defaults
default('clean')
vis=avg_data
imagename=first_image
interactive=True
spw = ''
mode = 'mfs'
nterms = 1
niter = 10000
threshold = '0.1mJy'
cell = '1.0arcsec'
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+'.*')
inp
go

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

viewer(first_image+'.image')

imview(raster={'file': first_image+'.image',
                       'range': [-0.1,0.7],
                       'colormap': 'RGB 1', 'scaling': 0.0, 'colorwedge': True},
               contour={'file': first_image+'.image',
                        'levels': [-1,1,2,3,4,6,8,10,12,14,16],
                        'unit': 0.05}, zoom=3)

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

# for NGC 3256
refant = 'DV06'
solint = '720s'
caltable = 'selfcal_ngc3256_line.gcal'

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

#
# Watch out for failed solutions noted in the terminal.
#

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

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

# It is useful to plot the X-Y solutions that give a good indication
# of the noise in the solutions

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

# The rms noise is about 4 to 8 deg, depending on antenna, but the
# phase changes are considerably larger.  This indicates that the
# application of this solution will improve the image.

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

flagmanager(vis=avg_data,
            mode='save',
            versionname='before_selfcal_apply')
applycal(vis=avg_data,
         gaintable=caltable,
         interp='linear',
         flagbackup=False)
# Use this if you have to redo something.
#flagmanager(vis=avg_data,
#            mode='restore',
#            versionname='before_selfcal_apply')

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

default('clean')
vis=avg_data
imagename=selfcal_image
interactive=True
swp = ''
mode = 'mfs'
nterms = 1
niter = 10000
threshold = '0.1mJy'
cell = '1.0arcsec'
imsize = 300
imagermode = 'csclean'
weighting = 'briggs'
robust = 0.5
os.system('rm -rf '+selfcal_image+'.*')
inp
go

# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%
# INSPECT THE NEW IMAGE AND COMPARE
# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%

imview(raster={'file': first_image+'.image',
               'range': [-0.1,0.7],
               'colormap': 'RGB 1', 
               'scaling': 0.0, 'colorwedge': True},
       contour={'file': first_image+'.image',
                'levels': [-1,1,2,3,4,6,8,10,12,14,16],
                'unit': 0.05}, zoom=3)

imview(raster={'file': selfcal_image+'.image',
               'range': [-0.1,0.8],
               'colormap': 'RGB 1', 'scaling': 0.0, 'colorwedge': True},
       contour={'file': selfcal_image+'.image',
                'levels': [-1,1,2,3,4,6,8,10,12,16,20,30,40,50],
                'unit': 0.02}, zoom=3)