###########################
# Initial Data Inspection #
###########################
listobs(vis='TwoQuasarSplit.ms',verbose=True)
plotants(vis='TwoQuasarSplit.ms')

####################
# Initial Flagging #
####################
    
# shadowed data
flagdata(vis='TwoQuasarSplit.ms', flagbackup=False, mode='shadow')

# autocorrelations
flagdata(vis='TwoQuasarSplit.ms', flagbackup=False, autocorr=True)

# save flags
flagmanager(vis='TwoQuasarSplit.ms', mode='save', versionname='Apriori')


#################
# More Flagging #
#################

# Discover Bad Edge Channels
plotms(vis='TwoQuasarSplit.ms', xaxis='freq', yaxis='amp', 
       averagedata=True, avgscan=True, avgtime='1e8',
       antenna='0&1',field='0') 

# Flag Edge Channels - 
flagdata(vis='TwoQuasarSplit.ms', flagbackup=False, mode='manualflag', 
         spw='*:0~5;61~63')

# Save Flag State
flagmanager(vis='TwoQuasarSplit.ms', mode='save', versionname='BeforeCal')

########################
# Bandpass Calibration #
########################

# Do Short-Timescale Phase Solution
gaincal(vis='TwoQuasarSplit.ms', caltable='bpphase.gcal',
                field='0', selectdata=True, gaintype='G',
                calmode='p', solint='int', refant='DV11',
                spw='0:15~25')

# Do frequency-dependent amplitude and phase solution
bandpass(vis='TwoQuasarSplit.ms', caltable='bandpass.bcal',
                 field='0', solint='inf', 
                 combine='scan', refant='DV11', bandtype='B', 
                 fillgaps=1, solnorm=True,
                 gaintable='bpphase.gcal', gainfield='0', interp='nearest')

####################
# Gain Calibration #
####################

setjy(vis='TwoQuasarSplit.ms', field='0', fluxdensity=[8.1, 0, 0, 0])

gaincal(vis='TwoQuasarSplit.ms', caltable='amp_phase_1.gcal',
       field='0', calmode='ap', solint='inf', refant='DV11', 
       gaintable='bandpass.bcal')

##################################
# Apply Calibrations And Inspect #
##################################

applycal(vis='TwoQuasarSplit.ms',field='0',
         gaintable=['bandpass.bcal','amp_phase_1.gcal'],
         gainfield=['0','0'],
         interp=['nearest','nearest'],
         flagbackup=False,calwt=True)

applycal(vis='TwoQuasarSplit.ms',field='1',
         gaintable=['bandpass.bcal','amp_phase_1.gcal'],
         gainfield=['0','0'],
         interp=['nearest','nearest'],
         flagbackup=False,calwt=True)

##################
# Flag More Data #
##################

# DV05
flagdata(vis='TwoQuasarSplit.ms', flagbackup=False, mode='manualflag', 
         antenna='DV05')

flagdata(vis='TwoQuasarSplit.ms', flagbackup=False, mode='manualflag', 
         antenna='DV11&DV12',spw='0:15~16')

# Save Flag State
flagmanager(vis='TwoQuasarSplit.ms', mode='save', versionname='AfterCal1')

####################
# Redo Calibration #
####################

bandpass(vis='TwoQuasarSplit.ms', caltable='bandpass2.bcal',
                 field='0', solint='inf', 
                 combine='scan', refant='DV10', bandtype='B', 
                 fillgaps=1, solnorm=True,
                 gaintable='bpphase.gcal', gainfield='0', interp='nearest')

gaincal(vis='TwoQuasarSplit.ms', caltable='amp_phase_2.gcal',
       field='0', calmode='ap', solint='inf', refant='DV10', 
       gaintable='bandpass2.bcal')

applycal(vis='TwoQuasarSplit.ms',field='0',
         gaintable=['bandpass2.bcal','amp_phase_2.gcal'],
         gainfield=['0','0'],
         interp=['nearest','nearest'],
         flagbackup=False,calwt=True)

applycal(vis='TwoQuasarSplit.ms',field='1',
         gaintable=['bandpass2.bcal','amp_phase_2.gcal'],
         gainfield=['0','0'],
         interp=['nearest','nearest'],
         flagbackup=False,calwt=True)

###########################
# Split Out The Good Data #
###########################

default(split)
split(vis='TwoQuasarSplit.ms', outputvis='1924-292.cal.ms',field='0')
split(vis='TwoQuasarSplit.ms', outputvis='1923-210.cal.ms',field='1')

#######################
# Selfcal on 1923-210 #
#######################

plotms(vis='1923-210.cal.ms')
gaincal(vis='1923-210.cal.ms', caltable='1923_selfcal.gcal',
       field='0', calmode='ap', solint='inf', refant='DV10')

plotcal(caltable='1923_selfcal.gcal',yaxis='phase')
plotcal(caltable='1923_selfcal.gcal',yaxis='amp')