# Script to look at and reduce eMERLIN data with CASA 4.2.2 or later
# This version starts with averaged and sorted data for 3C277.1 May 2015

# See  http://www.e-merlin.ac.uk/data_red/CASA/3C277.1_cal.html
# THIS VERSION HAS ALL THE ANSWERS
# Please try filling in 3C277.1_cal_outline.py first, use this just to
# check or if you are really stuck.

# Data can be downloaded as
# wget -c http://almanas.jb.man.ac.uk/amsr/3C277p1/all_avg.ms.tar
# wget -c http://almanas.jb.man.ac.uk/amsr/3C277p1/3C286_C.clean.model.tt0.tar
# wget -c http://almanas.jb.man.ac.uk/amsr/3C277p1/all_avg_target.flags

# Steps 1-3 are not included here as they are the routine preparation
# of the full data set (many tens G) including sorting and averaging.

# Steps 4 and 5 must be run first to flag bad data.  

# The calibration starts with step 6 

#############################
# ENTER THE NUMBER OF THE STEP(s) TO RUN (e.g. [6,7]
# To run interactively, define:
# mysteps = [s1, s2, ..] with the steps you want to run

step_title = {#1: 'Import data and list',
              #2: 'Sort and fixvis data and concatenate',
              #3: 'Time and channel average data',
              4: 'Check averaged data are present, inspect, quack and flag by channel',
              5: 'Flagging by time',
              6: 'Delay calibration',
              7: 'Setting the flux of 1331+305',
              8: 'Initial bandpass calibration, pre-gaincal inspection',
              9: 'Gain calibration of calibration sources',
              10: 'Determine and set fluxes of other calibrators',
              11: 'Derive bandpass with spectral index',
              12: 'Amp gain calibration of phase ref',
              13: 'Apply calibration from cal sources to all sources',
              14: 'Split out target'}


bases=['0319+415','1252+5634','1302+5748','1331+305','1407+284']

target = '1252+5634' # 3C277.1
phref  = '1302+5748'  # J1302+5748
bpcal  = '0319+415,1407+284' # 3C84, OQ208
leakcal = '0319+415'  # 3C84 - if not observed, use BP cal
fluxcal = '1331+305'  # 3C286, also pol angle cal.

calsources=phref+','+bpcal+','+fluxcal  # list of unique sources

antref='Mk2'          # refant - usually - change if known to be bad.
eMfactor =  0.938     # Flux scale factor if VLA values for 3C286 are used

gui = True            # False makes plot files only; True displays 'live' (needs interaction)

thesteps = []
try:
    print 'List of steps to be executed ...', mysteps
    thesteps = mysteps
except:
    print 'global variable mysteps not set'

if (thesteps==[]):
    thesteps = range(0,15)
    print 'Executing all steps: ', thesteps


### 4) Check data, list and identify end channels to flag
mystep = 4

if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

# Check whether data need loading
    if not os.path.exists('all_avg.ms'):
        print 'Please check that all_avg.ms is in this directory'
        exit()

# Listobs 
# You can uncomment for Linux
    os.system('rm -rf all_avg.listobs.txt')
    listobs(vis='all_avg.ms',
            overwrite=True,
            listfile='all_avg.listobs.txt')

    plotants(vis='all_avg.ms', figfile='plotants.png')

    # Default will be to use plotms
    plotuv(vis='all_avg.ms',
           maxnpts=10000000, symb='.',
           field='1302+5748', figfile='1302+5748_uv.png')

# Back up original flag state
    flagmanager(vis='all_avg.ms',
                mode='save',
                versionname='pre_quack')

    # Flag first quackinterval of each scan, before antennas all on source
    # Ascertain required interval in plotms
    # Inspect phase ref amp v. time, LL,RR, 
    # average central chans, zoom in manually

    plotms(vis='all_avg.ms',
           field='1302+5748',
           spw='0~3:20~43',
           avgchannel='24',
           xaxis='time',yaxis='amp',
           antenna=antref+'&Pi',           # change if necessary
           correlation='LL,RR',
           coloraxis='corr',
           showgui=gui,                    # change to F if necessary
           overwrite=True,plotfile='pre-quack.png')

# Flag first 40 s of each scan for all sources
    flagdata(vis='all_avg.ms',
                 mode='quack',quackinterval=40)

# re-plot to check
    plotms(vis='all_avg.ms',
           field='',                      # check all sources
           spw='0~3:20~43',
           avgchannel='24',
           xaxis='time',yaxis='amp',
           antenna=antref+'&Pi',           # change if necessary
           correlation='LL,RR',
           coloraxis='field',            # Now separate fields
           showgui=gui,                    # change to F if necessary
           overwrite=True,plotfile='quacked.png')

# Identify end channels to flag (same for all baselines/fields)
    plotms(vis='all_avg.ms',
       field='0319+415',
       gridrows=4,gridcols=1,
       xaxis='channel',yaxis='amp',
       antenna=antref+'&Kn',           # change if necessary
       correlation='LL,RR',   
       avgtime='99999',    # if data mostly good
       avgscan=True,
       iteraxis='spw', coloraxis='corr',
       showgui=gui,                          # change to T if you want ,
       overwrite=True,plotfile='0319+415_avg_amp-chan.png')
    os.system('mv 0319+415_avg_amp-chan?*.png 0319+415_avg_amp-chan.png')

# Back up flags
    flagmanager(vis='all_avg.ms',
                mode='save',
                versionname='pre_endchans')

 # end chans
    flagdata(vis='all_avg.ms',mode='manual',
             spw='0:0~6;60~63,1:0~3;61~63,2:0~3;61~63,3:0~3;63')

#  plot again to confirm
    plotms(vis='all_avg.ms',
       field='0319+415',
       gridrows=4,gridcols=1,
       xaxis='channel',yaxis='amp',
       antenna=antref+'&Kn',           # change if necessary
       correlation='LL,RR',
       avgtime='99999',    # if data mostly good
       avgscan=True,
       iteraxis='spw', coloraxis='corr',
       showgui=gui,                          # change to T if you want ,
       overwrite=True,plotfile='0319+415_flagged_avg_amp-chan.png')
    os.system('mv 0319+415_flagged_avg_amp-chan?*.png 0319+415_flagged_avg_amp-chan.png')

# If you notice any other bad channels, check these in more detail

### 5) Flag bad data by timerange
mystep = 5

if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

    flagmanager(vis='all_avg.ms',
                mode='save',
                versionname='pre_timeflagging')

# Identify obvious bad data 
# If target is faint, don't flag until after some calibration
# Here, it is bright enough to tell good from bad!
    # Generate all plots automatically: 
    for f in bases:
        plotms(vis='all_avg.ms',
               gridrows=5,gridcols=1,
               field=f,
               avgchannel='64',
               xaxis='time',yaxis='amp',
               antenna=antref+'&*',           # change if necessary
               correlation='LL,RR',
               iteraxis='baseline', coloraxis='corr',
               showgui=False,                    # change to True to inspect closely
               overwrite=True,plotfile=f+'_avg_amp-time.png')

        os.system('mv '+f+'_avg_amp-time?*.png '+f+'_avg_amp-time.png')

    # Plot one particular source:
    plotms(vis='all_avg.ms',
           gridrows=5,gridcols=1,
           field='0319+415',
           avgchannel='64',
           xaxis='time',yaxis='amp',
           antenna=antref+'&*',           # change if necessary
           correlation='LL,RR',
           iteraxis='baseline', coloraxis='corr',
           showgui=gui)

# If a source shows bad data, plot interactively and zoom/use locate
# to find the time/antenna affected.
# Write flag command file 'all_avg_1.flags'

    flagdata(vis='all_avg.ms',
            mode='list',
            inpfile='all_avg_1.flags')

    for f in bases:
        plotms(vis='all_avg.ms',
               gridrows=5,gridcols=1,
               field=f,
               avgchannel='64',
               xaxis='time',yaxis='amp',
               antenna=antref+'&*',           # change if necessary
               correlation='LL,RR',
               iteraxis='baseline', coloraxis='corr',
               showgui=False,                 # change to True to inspect closely
               overwrite=True,plotfile=f+'_flagged_avg_amp-time.png')

        os.system('mv '+f+'_flagged_avg_amp-time?*.png '+f+'_flagged_avg_amp-time.png')


### CALIBRATION

### 6) Delay calibration
mystep = 6

if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

# If phase-ref is weak, do bp cals only first and apply to all
# Here, phase ref is OK (see all sources have similar solutions).

# Plot baselines to refant for one bright source
# Phase against freq
    plotms(vis='all_avg.ms', field='0319+415', xaxis='frequency', 
           gridrows=5,gridcols=1,iteraxis='baseline',
           yaxis='phase',ydatacolumn='data',antenna=antref+'&*',
           avgtime='60000',avgscan=True,correlation='LL,RR',
           coloraxis='corr',
           plotfile='0319+415_phase_pre-delay.png',
           showgui=gui,                    # change to T to inspect closely
           overwrite=True)
    os.system('mv 0319+415_phase_pre-delay?*.png 0319+415_phase_pre-delay.png')

# Phase against freq per 10-min on bright source to test for 
# time stability
    plotms(vis='all_avg.ms', field='0319+415', xaxis='frequency', 
           yaxis='phase',ydatacolumn='data',antenna=antref+'&Cm',
           avgtime='600',correlation='RR',
           coloraxis='time',plotfile='0319+415_phase10min_pre-delay.png',
           showgui=gui,                    # change to T to inspect closely
           overwrite=True)

# amp against time, longest baseline, average all channels
    plotms(vis='all_avg.ms', field='0319+415', xaxis='time', 
           spw='0~3',avgchannel='64',plotrange=[-1,-1,0,0.018],
           yaxis='amp',ydatacolumn='data',antenna=antref+'&Cm',
           avgtime='60',correlation='RR',
           coloraxis='spw',plotfile='0319+415_amp-allch_pre-delay.png',
           showgui=gui,                    # change to T to inspect closely
           overwrite=True)

# amp against time, longest baseline, single channel
    plotms(vis='all_avg.ms', field='0319+415', xaxis='time', 
           spw='0~3:55', plotrange=[-1,-1,0,0.018],
           yaxis='amp',ydatacolumn='data',antenna=antref+'&Cm',
           avgtime='60',correlation='RR',
           coloraxis='spw',plotfile='0319+415_amp-1ch_pre-delay.png',
           showgui=gui,                    # change to T to inspect closely
           overwrite=True)

# Derive gain solutions
    os.system('rm -rf all_avg.K')
    gaincal(vis='all_avg.ms',
            gaintype='K',
            field=calsources,
            caltable='all_avg.K',
            spw='0~3',solint='600',
            refant=antref,
            minblperant=2,
            minsnr=2)

    plotcal(caltable='all_avg.K',
            xaxis='time',
            yaxis='delay',
            subplot=321,
            iteration='antenna',
            figfile='all_avg.K.png')


### 7) Setting the flux of 1331+305
mystep = 7

if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

# This takes the e-MERLIN model (centred on 5.5 GHz) and scales it to
# the zero spacing flux density (corected for spectral index and
# curvature) at the observed frequencies, centred on 5.073 GHz.
# However, because some of the flux is resolved-out by e-MERLIN, the
# derived flux densities will later be rescaled by 0.938 to allow for
# this.

    setjy(vis='all_avg.ms',
          field='1331+305',
          standard='Perley-Butler 2013',
          model='3C286_C.clean.model.tt0')
    
    os.system('rm -rf 3C286_model.png')
    plotms(vis='all_avg.ms', field='1331+305', xaxis='uvwave', 
           yaxis='amp', coloraxis='spw',ydatacolumn='model',
           correlation='LL,RR',
           showgui=gui,                    # change to T to inspect closely
           symbolsize=4, plotfile='3C286_model.png')
 

### 8) Initial bandpass calibration
mystep = 8

if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]
    
# Pre-calibrate phase and amp as a function of time

# First, phase: plot to estimate suitable solution interval
# constrained by longest baseline
    plotms(vis='all_avg.ms', field='0319+415', xaxis='time', 
           spw='0~3:55', 
           yaxis='phase',ydatacolumn='data',antenna=antref+'&Cm',
           coloraxis='spw',plotfile='0319+415_1ch_phase-time.png',
           showgui=gui,
           overwrite=True)
  
# Zoom in to help decide what phase solution interval to use
# The similar plot made for amplitude in Step 6 (delay) shows
# that the amplitudes are stable for longer.

    os.system('rm -rf bpcals_precal.p1')
    gaincal(vis='all_avg.ms',
        calmode='p',
        field=bpcal,
        caltable='bpcals_precal.p1',
        solint='30s',
        refant=antref,
        minblperant=2,
        gaintable='all_avg.K',
        minsnr=2)

# Check solutions are not just noise
    plotcal(caltable='bpcals_precal.p1',
            xaxis='time',
            yaxis='phase',
            subplot=321,
            plotrange=[-1,-1,-180,180],
            iteration='antenna',
            figfile='bpcals_precal.p1.png')

# Then amp and phase.
# Have to do sources separately because of normalisation

    os.system('rm -rf bpcals_precal.ap1')
    for bp in bpcal.split(','):
        gaincal(vis='all_avg.ms',
                calmode='ap',
                field=bp,
                caltable='bpcals_precal.ap1',
                solint='120s',
                solnorm=True,
                refant=antref,
                minblperant=2,
                append=True,
                gaintable=['all_avg.K','bpcals_precal.p1'],
                minsnr=2)

# Check there are no wild solutions
    plotcal(caltable='bpcals_precal.ap1',
            xaxis='time',
            yaxis='amp',
            subplot=321,
            iteration='antenna',
            figfile='bpcals_precal.ap1.png')

    os.system('rm -rf bpcal.B1')
    bandpass(vis='all_avg.ms',
        caltable='bpcal.B1',
        field=bpcal,
        fillgaps=16,
        solint='inf',combine='scan',
        solnorm=True,
        refant=antref,
        minblperant=2,
        gaintable=['all_avg.K','bpcals_precal.p1','bpcals_precal.ap1'],
        minsnr=3)

    # Plot bandpass solutions
    plotcal(caltable='bpcal.B1',xaxis='freq',
        plotrange=[-1,-1,-180,180],
        iteration='antenna',yaxis='phase',subplot=321,
        figfile='bpcal.B1_phase.png')
    plotcal(caltable='bpcal.B1',xaxis='freq',
        iteration='antenna',yaxis='amp',subplot=321,
        figfile='bpcal.B1_amp.png')    
 
# Apply delay and bandpass solutions to calibration sources to check
# bandpass calibration effect.

    applycal(vis='all_avg.ms',
             field=calsources,
             calwt=False,
             applymode='calonly',
             gaintable=['all_avg.K','bpcal.B1'],
             interp=['linear','nearest,linear'])
             
# Plot baselines to refant
    plotms(vis='all_avg.ms', field='1302+5748', xaxis='frequency', 
           gridrows=5,gridcols=1,iteraxis='baseline',
           yaxis='phase',ydatacolumn='data',antenna=antref+'&*',
           avgtime='600',correlation='LL,RR',overwrite=True,
           showgui=gui,                    # change to T to inspect closely
           coloraxis='scan',plotfile='1302+5748_beforeBP_phase.png')
    os.system('mv 1302+5748_beforeBP_phase?*.png 1302+5748_beforeBP_phase.png')

    plotms(vis='all_avg.ms', field='1302+5748', xaxis='frequency', 
           gridrows=5,gridcols=1,iteraxis='baseline',
           yaxis='amp',ydatacolumn='data',antenna=antref+'&*',
           avgtime='600',correlation='LL,RR',overwrite=True,
           showgui=gui,                    # change to T to inspect closely
           coloraxis='scan',plotfile='1302+5748_beforeBP_amp.png')
    os.system('mv 1302+5748_beforeBP_amp?*.png 1302+5748_beforeBP_amp.png')
 
    plotms(vis='all_avg.ms', field='1302+5748', xaxis='frequency', 
           gridrows=5,gridcols=1,iteraxis='baseline',
           yaxis='phase',ydatacolumn='corrected',antenna=antref+'&*',
           avgtime='600',correlation='LL,RR',overwrite=True,
           showgui=gui,                    # change to T to inspect closely
           coloraxis='scan',plotfile='1302+5748_afterBP_phase.png')
    os.system('mv 1302+5748_afterBP_phase?*.png 1302+5748_afterBP_phase.png')

    plotms(vis='all_avg.ms', field='1302+5748', xaxis='frequency', 
           gridrows=5,gridcols=1,iteraxis='baseline',
           yaxis='amp',ydatacolumn='corrected',antenna=antref+'&*',
           avgtime='600',correlation='LL,RR',overwrite=True,
           showgui=gui,                    # change to T to inspect closely
           coloraxis='scan',plotfile='1302+5748_afterBP_amp.png')
    os.system('mv 1302+5748_afterBP_amp?*.png 1302+5748_afterBP_amp.png')

# Use BP-corrected data to examine time-dependent fluctuations
# Use showgui=True to check for reasonable solution intervals
    plotms(vis='all_avg.ms', field='1302+5748', xaxis='time', 
           yaxis='phase',ydatacolumn='corrected',antenna=antref+'&*',
           iteraxis='baseline',avgchannel='64',
           correlation='LL',overwrite=True,
           showgui=gui,                    # change to T to inspect closely           
           coloraxis='spw',plotfile='1302+5748_afterBP_phase-time.png' )
    os.system('mv 1302+5748_afterBP_phase-time?*.png 1302+5748_afterBP_phase-time.png')

# Zoom in on a few scans.  Inspect several baselines including Mk2-Cm.
# Also inspect amplitude

### 9) Gain calibration of calibration sources
mystep = 9

if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

# Delete any previous models for sources other than 1331+305
    delmod(vis='all_avg.ms',
          field='0319+415,1407+284,1302+5748,1252+5634')

# Derive phase solutions
    
    os.system('rm -rf calsources.p1')
    gaincal(vis='all_avg.ms',
            calmode='p',
            caltable='calsources.p1',
            field=calsources,
            solint='16s',
            refant=antref,
            minblperant=2,minsnr=2,
            gaintable=['all_avg.K','bpcal.B1'],
            interp=['linear','nearest,linear'])
 
    # Plot solutions
    plotcal(caltable='calsources.p1',xaxis='time',
        plotrange=[-1,-1,-180,180],
        iteration='antenna',yaxis='phase',subplot=321,
        poln='L',figfile='calsources.p1_L.png')

    plotcal(caltable='calsources.p1',xaxis='time',
        plotrange=[-1,-1,-180,180],
        iteration='antenna',yaxis='phase',subplot=321,
        poln='R',figfile='calsources.p1_R.png')

# Derive amp solutions

    os.system('rm -rf calsources.ap1')
    gaincal(vis='all_avg.ms',
            calmode='ap',
            caltable='calsources.ap1',
            field=calsources,
            solint='180s',
            refant=antref,
            minblperant=2,minsnr=2,
            gaintable=['all_avg.K','bpcal.B1','calsources.p1'],
            interp=['linear','nearest,linear','nearest'])
 
    # Plot solutions
    plotcal(caltable='calsources.ap1',xaxis='time',
        iteration='antenna',yaxis='amp',subplot=321,
        poln='L',figfile='calsources.ap1_L.png')

    plotcal(caltable='calsources.ap1',xaxis='time',
        iteration='antenna',yaxis='amp',subplot=321,
        poln='R',figfile='calsources.ap1_R.png')

### 10) Determine and set the fluxes of the calibration sources
mystep = 10

if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

# Gains for 3C286 contain scale from correlator counts to known
# flux density in Jy
# Compare these with gains for other calibrators to estimate
# their flux densities
# Choose gainthreshold to get good SNR (trial and error)
# Too high a threshold - divergent gains add noise
# Too low - little data left, also worsens noise

    os.system('rm -rf calsources.ap1_flux calsources_flux.txt')
    calfluxes=fluxscale(vis='all_avg.ms',
                        caltable='calsources.ap1',
                        fluxtable='calsources.ap1_flux',
                        listfile='calsources_flux.txt',
                        gainthreshold=0.3,
                        antenna='!De',      # Least sensitive
                        reference='1331+305')

    eMcalfluxes={} 
    for k in calfluxes.keys():
        if len(calfluxes[k]) > 4:
            a=[]
            a.append(calfluxes[k]['fitFluxd']*eMfactor)
            a.append(calfluxes[k]['spidx'][0])
            a.append(calfluxes[k]['fitRefFreq'])
            eMcalfluxes[calfluxes[k]['fieldName']]=a
# eMcalfluxes now contains the scaled flux density for each source

# Use these to set flux densities of calibrators
# See messages

    for f in eMcalfluxes.keys():
        setjy(vis='all_avg.ms',
              field=f,
              standard='manual',
              fluxdensity=eMcalfluxes[f][0],
              spix=eMcalfluxes[f][1],
              reffreq=str(eMcalfluxes[f][2])+'Hz')

# Plot the  models we have just set

    plotms(vis='all_avg.ms', field=bpcal, xaxis='frequency', 
           yaxis='amp',ydatacolumn='model',
           coloraxis='spw',correlation='LL',
           customsymbol=True,symbolshape='circle',
           showgui=gui,                    # change to T to inspect closely
           symbolsize=5,plotfile='bpcal_model.png',overwrite=True)
 
 
### 11) Derive bandpass with spectral index
mystep = 11

if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

    os.system('rm -rf bpcal.B2')
    bandpass(vis='all_avg.ms',
        caltable='bpcal.B2',
        field=bpcal,
        fillgaps=16,
        solint='inf',combine='scan',
        refant=antref,
        minblperant=2,
        solnorm=False,
        gaintable=['all_avg.K','calsources.p1','calsources.ap1_flux'],
        gainfield=bpcal,
        minsnr=3)

    # Plot bandpass solutions
    plotcal(caltable='bpcal.B2',xaxis='freq',
        plotrange=[-1,-1,-180,180],
        iteration='antenna',yaxis='phase',subplot=321,
        figfile='bpcal.B2_phase.png')
    plotcal(caltable='bpcal.B2',xaxis='freq',
        iteration='antenna',yaxis='amp',subplot=321,
        figfile='bpcal.B2_amp.png')    
 
### 12) Derive new amp solutions for all cal sources  with improved bandpass table
mystep = 12

if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

    os.system('rm -rf calsources.ap2')
    gaincal(vis='all_avg.ms',
            calmode='ap',
            caltable='calsources.ap2',
            field=calsources,
            solint='inf',
            refant=antref,
            minblperant=2,minsnr=2,
            gaintable=['all_avg.K','bpcal.B2','calsources.p1'],
            interp=['linear','nearest,linear','nearest'])
 
    # Plot solutions
    plotcal(caltable='calsources.ap2',xaxis='time',
        iteration='antenna',yaxis='amp',subplot=321,
        plotrange=[-1,-1,0.01,0.06],
        poln='L',figfile='calsources.ap2_L.png')

    plotcal(caltable='calsources.ap2',xaxis='time',
        iteration='antenna',yaxis='amp',subplot=321,
        plotrange=[-1,-1,0.01,0.06],
        poln='R',figfile='calsources.ap2_R.png')


    # Generate long solution interval phase solutions to be applied to the
    # targets only
    os.system('rm -rf calsources.p2')
    gaincal(vis='all_avg.ms',
            calmode='p',
            caltable='calsources.p2',
            field=calsources,
            solint='inf',
            refant=antref,
            minblperant=2,minsnr=2,
            gaintable=['all_avg.K','bpcal.B2'],
            interp=['linear','nearest,linear'])

    # Plot solutions
    plotcal(caltable='calsources.p2',xaxis='time',
        plotrange=[-1,-1,-180,180],
        iteration='antenna',yaxis='phase',subplot=321,
        poln='L',figfile='calsources.p2_L.png')

    plotcal(caltable='calsources.p2',xaxis='time',
        plotrange=[-1,-1,-180,180],
        iteration='antenna',yaxis='phase',subplot=321,
        poln='R',figfile='calsources.p2_R.png')

### 13) Apply solutions to target sources and calibrators
mystep = 13

if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]

    cals=bases
    cals.remove(target)
    for c in cals:
        applycal(vis='all_avg.ms',
                 field=c,
                 gainfield=[c, '',c,c],
                 calwt=False,
                 applymode='calflag',
                 gaintable=['all_avg.K','bpcal.B2','calsources.p1','calsources.ap2'],
                 interp=['linear','nearest,linear','linear','linear'],
                 flagbackup=False)

# Now apply corrections to the target
    applycal(vis='all_avg.ms',
             field='1252+5634',
             gainfield=['1302+5748','','1302+5748','1302+5748'],
             calwt=False,
             applymode='calflag',  # Not too many failed solutions - OK to flag target data
             gaintable=['all_avg.K','bpcal.B2','calsources.p2','calsources.ap2'],
             interp=['linear','nearest,linear','linear','linear'],
             flagbackup=True)

# Plotms - first check phase-ref - is it point-like?
    plotms(vis='all_avg.ms', field='1302+5748', xaxis='uvdist', 
           yaxis='amp',ydatacolumn='corrected',
           avgchannel='64',correlation='LL,RR',
           overwrite=True,
           showgui=gui,                    # change to T to inspect closely
           coloraxis='spw',plotfile='1302+5748_amp-uvdist.png')

# If there is just a little bad data, no need to worry if the target
# is bright enough to self-calibrate

# Check target
    plotms(vis='all_avg.ms', field='1252+5634', xaxis='uvdist', 
           yaxis='amp',ydatacolumn='corrected',
           avgchannel='64',correlation='LL,RR',
           overwrite=True,
           showgui=gui,                    # change to T to inspect closely
           coloraxis='spw',plotfile='1252+5634_amp-uvdist.png')

### 14) Split out corrected target
mystep = 14

if(mystep in thesteps):
    print 'Step ', mystep, step_title[mystep]
    
    os.system('rm -rf 1252+5634.ms*')
    split(vis='all_avg.ms', 
          field='1252+5634', 
          outputvis='1252+5634.ms',
          keepflags=True)


    listobs(vis='1252+5634.ms',
            overwrite=True,
            listfile='1252+5634.listobs.txt')

### Now make sure that you save 1252+5634.ms for the imaging stages