# Sample TDM ALMA Data Reduction Script

# Calibration

thesteps = []
step_title = {0: 'listobs',
              1: 'A priori flagging',
              2: 'Generation and time averaging of the WVR cal table',
              3: 'Generation of the Tsys cal table',
              4: 'Generation of the antenna position cal table',
              5: 'Application of the WVR, Tsys and antpos cal tables',
              6: 'Split out science SPWs and time average',
              7: 'Listobs, clear pointing table, and save original flags',
              8: 'Initial flagging',
              9: 'Putting a model for the flux calibrator(s)',
              10: 'Save flags before bandpass cal',
              11: 'Bandpass calibration',
              12: 'Save flags before gain cal',
              13: 'Gain calibration',
              14: 'Save flags before applycal cal',
              15: 'Application of the bandpass and gain cal tables'}

try:
  print 'List of steps to be executed ...', mysteps
  thesteps = mysteps
except:
  print 'global variable mysteps not set.'
if (thesteps==[]):
  thesteps = range(0,len(step_title))
  print 'Executing all steps: ', thesteps

# The Python variable 'mysteps' will control which steps
# are executed when you start the script using
#   execfile('scriptForCalibration.py')
# e.g. setting
#   mysteps = [2,3,4]# before starting the script will make the script execute
# only steps 2, 3, and 4
# Setting mysteps = [] will make it execute all steps.

import re

es = aU.stuffForScienceDataReduction() 


if re.search('^3.4', casadef.casa_version) == None:
 sys.exit('ERROR: PLEASE USE THE SAME VERSION OF CASA THAT YOU USED FOR GENERATING THE SCRIPT: 3.4')

# Using reference antenna = DV04

print "# A priori calibration"

# listobs
mystep = 0
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  os.system('rm -rf example.ms.listobs')
  listobs(vis = 'example.ms',
    listfile = 'example.ms.listobs')
  
  

# A priori flagging
mystep = 1
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  tflagdata(vis = 'example.ms',
    mode = 'manual',
    autocorr = T,
    flagbackup = F)
  
  tflagdata(vis = 'example.ms',
    mode = 'manual',
    intent = '*POINTING*,*SIDEBAND_RATIO*,*ATMOSPHERE*',
    flagbackup = F)
  
  flagcmd(vis = 'example.ms',
    inpmode = 'table',
    action = 'plot')
  
  flagcmd(vis = 'example.ms',
    inpmode = 'table',
    action = 'apply')
  

# Generation and time averaging of the WVR cal table
mystep = 2
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  os.system('rm -rf example.ms.wvr') 
  
  mylogfile = casalog.logfile()
  casalog.setlogfile('example.ms.wvrgcal')
  
  wvrgcal(vis = 'example.ms',
    caltable = 'example.ms.wvr',
    tie = ['phasecal1', 'source1'],
    statsource = 'source1')
  
  casalog.setlogfile(mylogfile)
  
  os.system('rm -rf example.ms.wvr.smooth') 
  
  smoothcal(vis = 'example.ms',
    tablein = 'example.ms.wvr',
    caltable = 'example.ms.wvr.smooth',
    smoothtype = 'mean',
    smoothtime = 6.048)
  
  
  aU.plotWVRSolutions(caltable='example.ms.wvr.smooth', spw='9', antenna='DV04',
    yrange=[-180,180],subplot=22, interactive=False,
    figfile='example.ms.wvr.smooth.plots/example.ms.wvr.smooth') 
  
  #Note: If you see wraps in these plots, try changing yrange or unwrap=True 
  #Note: If all plots look strange, it may be a bad WVR on the reference antenna.
  #      To check, you can set antenna='' to show all baselines.
  

# Generation of the Tsys cal table
mystep = 3
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  os.system('rm -rf example.ms.tsys') 
  gencal(vis = 'example.ms',
    caltable = 'example.ms.tsys',
    caltype = 'tsys')
  
  aU.plotbandpass(caltable='example.ms.tsys', overlay='time', 
    xaxis='freq', yaxis='amp', subplot=22, buildpdf=False, interactive=False,
    showatm=True,pwv='auto',chanrange='5~123',showfdm=True, 
    field='', figfile='example.ms.tsys.plots.overlayTime/example.ms.tsys') 
  
  
  es.checkCalTable('example.ms.tsys', msName='example.ms', interactive=False) 
  

# Generation of the antenna position cal table
mystep = 4
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  # Position for antenna DV18 is derived from baseline run made on 2012-08-08 03:45:47.
  
  # Position for antenna DA44 is derived from baseline run made on 2012-08-08 03:45:47.
  
  # Position for antenna DV13 is derived from baseline run made on 2012-08-08 03:45:47.
  
  # Position for antenna DV12 is derived from baseline run made on 2012-08-08 03:45:47.
  
  # Position for antenna DA41 is derived from baseline run made on 2012-08-08 03:45:47.
  
  # Position for antenna DV14 is derived from baseline run made on 2012-08-08 03:45:47.
  
  # Position for antenna DA43 is derived from baseline run made on 2012-08-08 03:45:47.
  
  # Position for antenna DA42 is derived from baseline run made on 2012-08-08 03:45:47.
  
  # Position for antenna PM01 is derived from baseline run made on 2012-08-08 03:45:47.
  
  # Position for antenna PM03 is derived from baseline run made on 2012-08-08 03:45:47.
  
  # Position for antenna DV10 is derived from baseline run made on 2012-08-08 03:45:47.
  
  # Position for antenna DV20 is derived from baseline run made on 2012-08-08 03:45:47.
  
  # Position for antenna DV21 is derived from baseline run made on 2012-08-08 03:45:47.
  
  # Position for antenna DA46 is derived from baseline run made on 2012-08-08 03:45:47.
  
  # Position for antenna DV08 is derived from baseline run made on 2012-08-08 03:45:47.
  
  # Position for antenna DV15 is derived from baseline run made on 2012-08-08 03:45:47.
  
  # Position for antenna DV04 is derived from baseline run made on 2012-08-08 03:45:47.
  
  # Position for antenna DV02 is derived from baseline run made on 2012-08-08 03:45:47.
  
  # Position for antenna DV03 is derived from baseline run made on 2012-08-08 03:45:47.
  
  # Position for antenna DV17 is derived from baseline run made on 2012-08-08 03:45:47.
  
  # Position for antenna DV16 is derived from baseline run made on 2012-08-08 03:45:47.
  
  os.system('rm -rf example.ms.antpos') 
  gencal(vis = 'example.ms',
    caltable = 'example.ms.antpos',
    caltype = 'antpos',
    antenna = 'PM01,DV18,PM03,DA46,DV08,DV21,DA44,DV13,DV12,DA41,DV14,DA43,DA42,DV02,DV03,DV20,DV17,DV10,DV16,DV04,DV15',
    parameter = [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0])
  #  parameter = [-0.000345313535421,0.000334442860567,0.000121842255689,-4.08322492433e-05,0.000467001455385,0.000209667279493,-4.77150030442e-05,0.000298935292451,0.000214559943634,-9.7899480286e-05,0.000294809177666,3.20817786836e-05,-2.88260530472e-05,0.000316137562566,3.69743181209e-05,-7.74986637655e-05,3.09850973321e-05,9.39674457339e-05,-0.000163363823899,0.000304330190095,5.27768551355e-05,-0.000234717159163,0.000399494944287,0.000163241063104,-9.40415857865e-05,0.000356702195474,0.000182657394833,-2.99885869026e-07,5.02914190292e-08,-4.80096787214e-07,6.53000743989e-05,-0.000112107882566,0.00013462437725,-0.000398558780189,0.000467274739713,0.000117129058565,-5.57055037029e-05,0.000175815160516,0.000119949681909,-0.00052189608713,0.000618181807822,0.000249023546509,-0.000228094450489,0.000420299422235,0.000200702061088,-7.59716868921e-05,0.000175179280833,0.000167444087258,-1.53797916934e-06,-2.79378547111e-05,0.000332344867411,-0.000272655006342,0.000381348849974,0.000125304836861,-6.9154646651e-05,0.000348717685481,0.000207042533502,-3.2614916563e-05,0.000306133180857,0.000139814335853,-6.42650146273e-05,0.000378984128862,0.000286478686756])
  

# Application of the WVR, Tsys and antpos cal tables
mystep = 5
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  
  
  from recipes.almahelpers import tsysspwmap
  tsysmap = tsysspwmap(vis = 'example.ms', tsystable = 'example.ms.tsys')
  
  
  
  applycal(vis = 'example.ms',
    field = '0',
    spw = '9,11,13,15',
    gaintable = ['example.ms.tsys', 'example.ms.wvr.smooth', 'example.ms.antpos'],
    gainfield = ['0', '', ''],
    interp = 'linear,linear',
    spwmap = [tsysmap,[],[]],
    calwt = T,
    flagbackup = F)
  
  
  
  applycal(vis = 'example.ms',
    field = '1',
    spw = '9,11,13,15',
    gaintable = ['example.ms.tsys', 'example.ms.wvr.smooth', 'example.ms.antpos'],
    gainfield = ['1', '', ''],
    interp = 'linear,linear',
    spwmap = [tsysmap,[],[]],
    calwt = T,
    flagbackup = F)
  
  
  
  applycal(vis = 'example.ms',
    field = '2',
    spw = '9,11,13,15',
    gaintable = ['example.ms.tsys', 'example.ms.wvr.smooth', 'example.ms.antpos'],
    gainfield = ['2', '', ''],
    interp = 'linear,linear',
    spwmap = [tsysmap,[],[]],
    calwt = T,
    flagbackup = F)
  

  
  applycal(vis = 'example.ms',
    field = '3',
    spw = '9,11,13,15',
    gaintable = ['example.ms.tsys', 'example.ms.wvr.smooth', 'example.ms.antpos'],
    gainfield = ['2', '', ''],
    interp = 'linear,linear',
    spwmap = [tsysmap,[],[]],
    calwt = T,
    flagbackup = F)
  
  
  
  applycal(vis = 'example.ms',
    field = '4',
    spw = '9,11,13,15',
    gaintable = ['example.ms.tsys', 'example.ms.wvr.smooth', 'example.ms.antpos'],
    gainfield = ['2', '', ''],
    interp = 'linear,linear',
    spwmap = [tsysmap,[],[]],
    calwt = T,
    flagbackup = F)
  
  
  
  
  
  applycal(vis = 'example.ms',
    field = '5',
    spw = '9,11,13,15',
    gaintable = ['example.ms.tsys', 'example.ms.wvr.smooth', 'example.ms.antpos'],
    gainfield = ['2', '', ''],
    interp = 'linear,linear',
    spwmap = [tsysmap,[],[]],
    calwt = T,
    flagbackup = F)
  
  
  
 
  
  applycal(vis = 'example.ms',
    field = '6',
    spw = '9,11,13,15',
    gaintable = ['example.ms.tsys', 'example.ms.wvr.smooth', 'example.ms.antpos'],
    gainfield = ['2', '', ''],
    interp = 'linear,linear',
    spwmap = [tsysmap,[],[]],
    calwt = T,
    flagbackup = F)
  
  
  
 
  
  applycal(vis = 'example.ms',
    field = '7',
    spw = '9,11,13,15',
    gaintable = ['example.ms.tsys', 'example.ms.wvr.smooth', 'example.ms.antpos'],
    gainfield = ['2', '', ''],
    interp = 'linear,linear',
    spwmap = [tsysmap,[],[]],
    calwt = T,
    flagbackup = F)
  
  
  
 
  
  applycal(vis = 'example.ms',
    field = '8',
    spw = '9,11,13,15',
    gaintable = ['example.ms.tsys', 'example.ms.wvr.smooth', 'example.ms.antpos'],
    gainfield = ['2', '', ''],
    interp = 'linear,linear',
    spwmap = [tsysmap,[],[]],
    calwt = T,
    flagbackup = F)
  
  
  
  
  
  applycal(vis = 'example.ms',
    field = '9',
    spw = '9,11,13,15',
    gaintable = ['example.ms.tsys', 'example.ms.wvr.smooth', 'example.ms.antpos'],
    gainfield = ['2', '', ''],
    interp = 'linear,linear',
    spwmap = [tsysmap,[],[]],
    calwt = T,
    flagbackup = F)
  
  
  
 
  
  applycal(vis = 'example.ms',
    field = '10',
    spw = '9,11,13,15',
    gaintable = ['example.ms.tsys', 'example.ms.wvr.smooth', 'example.ms.antpos'],
    gainfield = ['2', '', ''],
    interp = 'linear,linear',
    spwmap = [tsysmap,[],[]],
    calwt = T,
    flagbackup = F)
  
  
  
 
  
  applycal(vis = 'example.ms',
    field = '11',
    spw = '9,11,13,15',
    gaintable = ['example.ms.tsys', 'example.ms.wvr.smooth', 'example.ms.antpos'],
    gainfield = ['2', '', ''],
    interp = 'linear,linear',
    spwmap = [tsysmap,[],[]],
    calwt = T,
    flagbackup = F)
  
  
  
 
  
  applycal(vis = 'example.ms',
    field = '12',
    spw = '9,11,13,15',
    gaintable = ['example.ms.tsys', 'example.ms.wvr.smooth', 'example.ms.antpos'],
    gainfield = ['2', '', ''],
    interp = 'linear,linear',
    spwmap = [tsysmap,[],[]],
    calwt = T,
    flagbackup = F)
  
  es.getCalWeightStats('example.ms') 
  

# Split out science SPWs and time average
mystep = 6
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  # Important note: the correlator mode for this dataset was TDM, you may want not to do any time averaging.
  
  os.system('rm -rf example.ms.split') 
  split(vis = 'example.ms',
    outputvis = 'example.ms.split',
    datacolumn = 'corrected',
    spw = '9,11,13,15',
    timebin = '6.048s',
    keepflags = T)
  
  

print "# Calibration"

# Listobs, clear pointing table, and save original flags
mystep = 7
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  os.system('rm -rf example.ms.split.listobs')
  listobs(vis = 'example.ms.split',
    listfile = 'example.ms.split.listobs')
  
  tb.open('example.ms.split/POINTING', nomodify = False)
  a = tb.rownumbers()
  tb.removerows(a)
  tb.close()
  flagmanager(vis = 'example.ms.split',
    mode = 'save',
    versionname = 'Original')
  
  

# Initial flagging
mystep = 8
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  # Flagging shadowed data
  
  tflagdata(vis = 'example.ms.split',
    mode = 'shadow',
    flagbackup = F)
  
  # Flagging edge channels
  
  tflagdata(vis = 'example.ms.split',
    mode = 'manual',
    spw = '0:0~7;120~127,1:0~7;120~127,2:0~7;120~127,3:0~7;120~127',
    flagbackup = F)
  
  # Flagging atmospheric line(s)
  
  tflagdata(vis = 'example.ms.split',
    mode = 'manual',
    spw = '0:0~109',
    field = '1',
    flagbackup = F)
  
  # Flagging atmospheric line(s)
  
  tflagdata(vis = 'example.ms.split',
    mode = 'manual',
    spw = '1:0~127',
    field = '1',
    flagbackup = F)
  
  

# Putting a model for the flux calibrator(s)
mystep = 9
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  setjy(vis = 'example.ms.split',
    field = '1', # Neptune
    spw = '0,1,2,3',
    standard = 'Butler-JPL-Horizons 2010')
  
  os.system('rm -rf example.ms.split.setjy.field*.png') 
  for i in ['1']:
    plotms(vis = 'example.ms.split',
      xaxis = 'uvdist',
      yaxis = 'amp',
      ydatacolumn = 'model',
      field = i,
      spw = '0,1,2,3',
      avgchannel = '9999',
      coloraxis = 'spw',
      plotfile = 'example.ms.split.setjy.field'+i+'.png')
  

# Save flags before bandpass cal
mystep = 10
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  flagmanager(vis = 'example.ms.split',
    mode = 'save',
    versionname = 'BeforeBandpassCalibration')
  
  

# Bandpass calibration
mystep = 11
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  os.system('rm -rf example.ms.split.ap_pre_bandpass') 
  
  gaincal(vis = 'example.ms.split',
    caltable = 'example.ms.split.ap_pre_bandpass',
    field = '0', # J1924-292
    spw = '0:51~76,1:51~76,2:51~76,3:51~76',
    solint = 'int',
    refant = 'DV04',
    calmode = 'ap')
  
  es.checkCalTable('example.ms.split.ap_pre_bandpass', msName='example.ms.split', interactive=False) 
  
  os.system('rm -rf example.ms.split.bandpass') 
  bandpass(vis = 'example.ms.split',
    caltable = 'example.ms.split.bandpass',
    field = '0', # J1924-292
    solint = 'inf',
    combine = 'scan',
    refant = 'DV04',
    solnorm = T,
    bandtype = 'B',
    gaintable = 'example.ms.split.ap_pre_bandpass')
  
  es.checkCalTable('example.ms.split.bandpass', msName='example.ms.split', interactive=False) 
  

# Save flags before gain cal
mystep = 12
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  flagmanager(vis = 'example.ms.split',
    mode = 'save',
    versionname = 'BeforeGainCalibration')
  
  

# Gain calibration
mystep = 13
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  # Note: the Solar system object used for flux calibration is highly resolved on some baselines.
  # Note: we will first determine the flux of the phase calibrator(s) on a subset of antennas.
  
  clearcal('example.ms.split',field='J2056-472')
  
  os.system('rm -rf example.ms.split.phase_short_int') 
  gaincal(vis = 'example.ms.split',
    caltable = 'example.ms.split.phase_short_int',
    field = '1', # Neptune
    selectdata = T,
    antenna = 'DA41,DV04,DV08,DV19,DV21&',
    solint = 'int',
    refant = 'DV04',
    gaintype = 'G',
    calmode = 'p',
    gaintable = 'example.ms.split.bandpass')
  
  gaincal(vis = 'example.ms.split',
    caltable = 'example.ms.split.phase_short_int',
    field = '0,2', # J1924-292,J2056-472
    selectdata = T,
    solint = 'int',
    refant = 'DV04',
    gaintype = 'G',
    calmode = 'p',
    append = T,
    gaintable = 'example.ms.split.bandpass')
  
  es.checkCalTable('example.ms.split.phase_short_int', msName='example.ms.split', interactive=False) 
  
  os.system('rm -rf example.ms.split.ampli_short_inf') 
  gaincal(vis = 'example.ms.split',
    caltable = 'example.ms.split.ampli_short_inf',
    field = '0,1,2', # J1924-292,Neptune,J2056-472
    selectdata = T,
    solint = 'inf',
    refant = 'DV04',
    gaintype = 'T',
    calmode = 'ap',
    gaintable = ['example.ms.split.bandpass', 'example.ms.split.phase_short_int'])
  
  es.checkCalTable('example.ms.split.ampli_short_inf', msName='example.ms.split', interactive=False) 
  
  os.system('rm -rf example.ms.split.flux_short_inf') 
  os.system('rm -rf example.ms.split.fluxscale') 
  mylogfile = casalog.logfile()
  casalog.setlogfile('example.ms.split.fluxscale')
  
  fluxscale(vis = 'example.ms.split',
    caltable = 'example.ms.split.ampli_short_inf',
    fluxtable = 'example.ms.split.flux_short_inf',
    reference = '1') # Neptune
  
  casalog.setlogfile(mylogfile)
  
  f = open('example.ms.split.fluxscale')
  fc = f.readlines()
  f.close()
  
  for phaseCalName in ['J2056-472']:
    for i in range(len(fc)):
      if fc[i].find('Flux density for '+phaseCalName) != -1 and re.search('in SpW=[0-9]+(?: \(ref SpW=[0-9]+\))? is: [0-9]+\.[0-9]+', fc[i]) != None:
        line = (re.search('in SpW=[0-9]+(?: \(ref SpW=[0-9]+\))? is: [0-9]+\.[0-9]+', fc[i])).group(0)
        spwId = (line.split('='))[1].split()[0]
        flux = float((line.split(':'))[1].split()[0])
        setjy(vis = 'example.ms.split',
          field = phaseCalName.replace(';','*;').split(';')[0],
          spw = spwId,
          fluxdensity = [flux,0,0,0])
  
  os.system('rm -rf example.ms.split.phase_int') 
  gaincal(vis = 'example.ms.split',
    caltable = 'example.ms.split.phase_int',
    field = '0,1,2', # J1924-292,Neptune,J2056-472
    solint = 'int',
    refant = 'DV04',
    gaintype = 'G',
    calmode = 'p',
    gaintable = 'example.ms.split.bandpass')
  
  es.checkCalTable('example.ms.split.phase_int', msName='example.ms.split', interactive=False) 
  
  os.system('rm -rf example.ms.split.flux_inf') 
  gaincal(vis = 'example.ms.split',
    caltable = 'example.ms.split.flux_inf',
    field = '0,1,2', # J1924-292,Neptune,J2056-472
    solint = 'inf',
    refant = 'DV04',
    gaintype = 'T',
    calmode = 'ap',
    gaintable = ['example.ms.split.bandpass', 'example.ms.split.phase_int'])
  
  es.checkCalTable('example.ms.split.flux_inf', msName='example.ms.split', interactive=False) 
  
  os.system('rm -rf example.ms.split.phase_inf') 
  gaincal(vis = 'example.ms.split',
    caltable = 'example.ms.split.phase_inf',
    field = '0,1,2', # J1924-292,Neptune,J2056-472
    solint = 'inf',
    refant = 'DV04',
    gaintype = 'G',
    calmode = 'p',
    gaintable = 'example.ms.split.bandpass')
  
  es.checkCalTable('example.ms.split.phase_inf', msName='example.ms.split', interactive=False) 
  

# Save flags before applycal cal
mystep = 14
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  flagmanager(vis = 'example.ms.split',
    mode = 'save',
    versionname = 'BeforeApplycal')
  
  

# Application of the bandpass and gain cal tables
mystep = 15
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  for i in ['0', '1']: # J1924-292,Neptune
    applycal(vis = 'example.ms.split',
      field = i,
      gaintable = ['example.ms.split.bandpass', 'example.ms.split.phase_int', 'example.ms.split.flux_inf'],
      gainfield = ['', i, i],
      interp = 'linear,linear',
      calwt = F,
      flagbackup = F)
  
  applycal(vis = 'example.ms.split',
    field = '2,3~12', # Science targets
    gaintable = ['example.ms.split.bandpass', 'example.ms.split.phase_inf', 'example.ms.split.flux_inf'],
    gainfield = ['', '2', '2'], # J2056-472
    interp = 'linear,linear',
    calwt = F,
    flagbackup = F)