# ALMA Data Reduction Script

# Calibration

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

if 'applyonly' not in globals(): applyonly = False
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

import os

if applyonly != True: es = aU.stuffForScienceDataReduction() 


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


# CALIBRATE_AMPLI: J1517-243
# CALIBRATE_ATMOSPHERE: IRAS16293-2422,J1517-243,J1700-2610
# CALIBRATE_BANDPASS: J1700-2610
# CALIBRATE_FLUX: J1517-243
# CALIBRATE_FOCUS: 
# CALIBRATE_PHASE: J1625-2527
# CALIBRATE_POINTING: J1700-2610
# OBSERVE_TARGET: IRAS16293-2422

# Using reference antenna = DA59

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

  if os.path.exists('uid___A002_X867766_Xa7.ms') == False:
    importasdm('uid___A002_X867766_Xa7', asis='Antenna Station Receiver Source CalAtmosphere CalWVR')
  if applyonly != True: es.fixForCSV2555('uid___A002_X867766_Xa7.ms')

# Fix of SYSCAL table times
mystep = 1
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  from recipes.almahelpers import fixsyscaltimes
  fixsyscaltimes(vis = 'uid___A002_X867766_Xa7.ms')

print "# A priori calibration"

# listobs
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 uid___A002_X867766_Xa7.ms.listobs')
  listobs(vis = 'uid___A002_X867766_Xa7.ms',
    listfile = 'uid___A002_X867766_Xa7.ms.listobs')
  
  

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

  flagdata(vis = 'uid___A002_X867766_Xa7.ms',
    mode = 'manual',
    spw = '1~24',
    autocorr = T,
    flagbackup = F)
  
  flagdata(vis = 'uid___A002_X867766_Xa7.ms',
    mode = 'manual',
    intent = '*POINTING*,*SIDEBAND_RATIO*,*ATMOSPHERE*',
    flagbackup = F)
  
  flagcmd(vis = 'uid___A002_X867766_Xa7.ms',
    inpmode = 'table',
    useapplied = True,
    action = 'plot',
    plotfile = 'uid___A002_X867766_Xa7.ms.flagcmd.png')
  
  flagcmd(vis = 'uid___A002_X867766_Xa7.ms',
    inpmode = 'table',
    useapplied = True,
    action = 'apply')
  

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

  os.system('rm -rf uid___A002_X867766_Xa7.ms.wvr') 
  
  os.system('rm -rf uid___A002_X867766_Xa7.ms.wvrgcal') 
  
  mylogfile = casalog.logfile()
  casalog.setlogfile('uid___A002_X867766_Xa7.ms.wvrgcal')
  
  wvrgcal(vis = 'uid___A002_X867766_Xa7.ms',
    caltable = 'uid___A002_X867766_Xa7.ms.wvr',
    toffset = 0,
    tie = ['IRAS16293-2422,J1625-2527'],
    statsource = 'IRAS16293-2422')
  
  casalog.setlogfile(mylogfile)
  
  # This is a temporary workaround, which will be included in a future version of CASA
  
  tb.open('uid___A002_X867766_Xa7.ms.wvr', nomodify=False)
  count = 0
  numrows = tb.nrows()
  mycparamcol = tb.getcol('CPARAM')
  for i in range(0, numrows):
      if mycparamcol[0][0][i] == (1.+0.j):
          tb.putcell('FLAG', i, [[True]])
          count += 1
  tb.close()
  del mycparamcol
  if(numrows>0):
      print 'Flagged', count, 'of', numrows, 'solutions =', 100.*count/float(numrows),'%'
  
  
  os.system('rm -rf uid___A002_X867766_Xa7.ms.wvr.smooth') 
  
  smoothcal(vis = 'uid___A002_X867766_Xa7.ms',
    tablein = 'uid___A002_X867766_Xa7.ms.wvr',
    caltable = 'uid___A002_X867766_Xa7.ms.wvr.smooth',
    smoothtype = 'mean',
    smoothtime = 6.048)
  
  
  if applyonly != True: aU.plotWVRSolutions(caltable='uid___A002_X867766_Xa7.ms.wvr.smooth', spw='17', antenna='DA59',
    yrange=[-399,399],subplot=22, interactive=False,
    figfile='uid___A002_X867766_Xa7.ms.wvr.smooth.plots/uid___A002_X867766_Xa7.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 = 5
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  os.system('rm -rf uid___A002_X867766_Xa7.ms.tsys') 
  gencal(vis = 'uid___A002_X867766_Xa7.ms',
    caltable = 'uid___A002_X867766_Xa7.ms.tsys',
    caltype = 'tsys')
  
  if applyonly != True: aU.plotbandpass(caltable='uid___A002_X867766_Xa7.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='uid___A002_X867766_Xa7.ms.tsys.plots.overlayTime/uid___A002_X867766_Xa7.ms.tsys') 
  
  
  if applyonly != True: es.checkCalTable('uid___A002_X867766_Xa7.ms.tsys', msName='uid___A002_X867766_Xa7.ms', interactive=False) 
  

# Generation of the antenna position cal table
mystep = 6
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 2014-06-11 04:48:50.
  
  # Note: no baseline run found for antenna DA65.
  
  # Position for antenna DA64 is derived from baseline run made on 2014-06-02 04:15:36.
  
  # Position for antenna DA62 is derived from baseline run made on 2014-06-11 04:48:50.
  
  # Position for antenna DA61 is derived from baseline run made on 2014-06-11 04:48:50.
  
  # Position for antenna DA60 is derived from baseline run made on 2014-06-11 04:48:50.
  
  # Position for antenna DV10 is derived from baseline run made on 2014-06-11 04:48:50.
  
  # Position for antenna DV13 is derived from baseline run made on 2014-06-11 04:48:50.
  
  # Position for antenna DA41 is derived from baseline run made on 2014-06-11 04:48:50.
  
  # Position for antenna DV14 is derived from baseline run made on 2014-06-11 04:48:50.
  
  # Position for antenna DA43 is derived from baseline run made on 2014-06-11 04:48:50.
  
  # Position for antenna DV16 is derived from baseline run made on 2014-06-11 04:48:50.
  
  # Position for antenna DV22 is derived from baseline run made on 2014-06-11 04:48:50.
  
  # Note: the correction for antenna DA56 is larger than 2mm.
  
  # Position for antenna DA56 is derived from baseline run made on 2013-08-17 08:36:51.
  
  # Position for antenna DA55 is derived from baseline run made on 2014-06-11 04:48:50.
  
  # Position for antenna DV15 is derived from baseline run made on 2014-06-11 04:48:50.
  
  # Position for antenna DV01 is derived from baseline run made on 2014-06-11 04:48:50.
  
  os.system('rm -rf uid___A002_X867766_Xa7.ms.antpos') 
  gencal(vis = 'uid___A002_X867766_Xa7.ms',
    caltable = 'uid___A002_X867766_Xa7.ms.antpos',
    caltype = 'antpos',
    antenna = 'DV18,DV22,DA64,DA56,DA62,DA61,DA60,DV10,DV13,DA55,DA41,DV14,DA43,DV16,DV01,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])
  #  parameter = [3.38209701852e-05,0.000162975135086,0.000104963075981,5.11011722337e-05,-6.68214898192e-05,8.94626105963e-05,2.32582675032e-05,0.000577052575402,-0.00020807301174,0.0007190015167,0.00154299940914,0.00110400049016,-0.000118697035043,8.94472765878e-05,-0.000153398821974,-2.19487514641e-06,0.000121593701559,-1.4005508361e-05,8.7917182889e-05,-0.000539778145669,-0.00061467074309,-6.12014207699e-05,3.34072436067e-05,0.000144410325705,-0.000207869911354,-0.000128378804152,-8.74984719765e-05,-3.06732186321e-05,-0.00039518403727,-0.000294974118095,-2.22678337143e-05,4.38745189799e-05,-0.000185537516869,3.20842733064e-05,0.000270253174919,0.000279583298594,-0.000189503312959,0.000590201703297,0.000258763296864,2.70871584607e-05,0.00025869211119,9.93365262099e-06,8.56135263886e-06,-0.000118660870292,-0.000148688604482,-2.69875650599e-05,0.00042259905275,0.000214001663422])
  

# Application of the WVR, Tsys and antpos cal tables
mystep = 7
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 = 'uid___A002_X867766_Xa7.ms', tsystable = 'uid___A002_X867766_Xa7.ms.tsys', tsysChanTol = 1)
  
  
  
  applycal(vis = 'uid___A002_X867766_Xa7.ms',
    field = '0',
    spw = '17,19,21,23',
    gaintable = ['uid___A002_X867766_Xa7.ms.tsys', 'uid___A002_X867766_Xa7.ms.wvr.smooth', 'uid___A002_X867766_Xa7.ms.antpos'],
    gainfield = ['0', '', ''],
    interp = 'linear,linear',
    spwmap = [tsysmap,[],[]],
    calwt = T,
    flagbackup = F)
  
  
  
  applycal(vis = 'uid___A002_X867766_Xa7.ms',
    field = '1',
    spw = '17,19,21,23',
    gaintable = ['uid___A002_X867766_Xa7.ms.tsys', 'uid___A002_X867766_Xa7.ms.wvr.smooth', 'uid___A002_X867766_Xa7.ms.antpos'],
    gainfield = ['1', '', ''],
    interp = 'linear,linear',
    spwmap = [tsysmap,[],[]],
    calwt = T,
    flagbackup = F)
  
  
  
  # Note: J1625-2527 didn't have any Tsys measurement, so I used the one made on IRAS16293-2422. This is probably Ok.
  
  applycal(vis = 'uid___A002_X867766_Xa7.ms',
    field = '2',
    spw = '17,19,21,23',
    gaintable = ['uid___A002_X867766_Xa7.ms.tsys', 'uid___A002_X867766_Xa7.ms.wvr.smooth', 'uid___A002_X867766_Xa7.ms.antpos'],
    gainfield = ['3', '', ''],
    interp = 'linear,linear',
    spwmap = [tsysmap,[],[]],
    calwt = T,
    flagbackup = F)
  
  
  
  applycal(vis = 'uid___A002_X867766_Xa7.ms',
    field = '3',
    spw = '17,19,21,23',
    gaintable = ['uid___A002_X867766_Xa7.ms.tsys', 'uid___A002_X867766_Xa7.ms.wvr.smooth', 'uid___A002_X867766_Xa7.ms.antpos'],
    gainfield = ['3', '', ''],
    interp = 'linear,linear',
    spwmap = [tsysmap,[],[]],
    calwt = T,
    flagbackup = F)
  
  if applyonly != True: es.getCalWeightStats('uid___A002_X867766_Xa7.ms') 
  

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

  os.system('rm -rf uid___A002_X867766_Xa7.ms.split') 
  split(vis = 'uid___A002_X867766_Xa7.ms',
    outputvis = 'uid___A002_X867766_Xa7.ms.split',
    datacolumn = 'corrected',
    spw = '17,19,21,23',
    keepflags = T)
  
  

print "# Calibration"

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

  os.system('rm -rf uid___A002_X867766_Xa7.ms.split.listobs')
  listobs(vis = 'uid___A002_X867766_Xa7.ms.split',
    listfile = 'uid___A002_X867766_Xa7.ms.split.listobs')
  
  tb.open('uid___A002_X867766_Xa7.ms.split/POINTING', nomodify = False)
  a = tb.rownumbers()
  tb.removerows(a)
  tb.close()
  
  if not os.path.exists('uid___A002_X867766_Xa7.ms.split.flagversions/Original.flags'):
    flagmanager(vis = 'uid___A002_X867766_Xa7.ms.split',
      mode = 'save',
      versionname = 'Original')
  
  

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

  # Flagging shadowed data
  
  flagdata(vis = 'uid___A002_X867766_Xa7.ms.split',
    mode = 'shadow',
    flagbackup = F)
  
  flagdata(vis = 'uid___A002_X867766_Xa7.ms.split',
    mode = 'manual',
    field = '2',
    scan = '30',
    antenna = 'DV10',
    flagbackup = F)

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

  setjy(vis = 'uid___A002_X867766_Xa7.ms.split',
    field = '1', # source name = J1517-243
    spw = '0', # center frequency of spw = 157.067484291GHz
    standard = 'manual',
    fluxdensity = [1.08134701397, 0, 0, 0]) # frequency of measurement = 148.760010655GHz
  
  setjy(vis = 'uid___A002_X867766_Xa7.ms.split',
    field = '1', # source name = J1517-243
    spw = '1', # center frequency of spw = 158.224154724GHz
    standard = 'manual',
    fluxdensity = [1.08134701397, 0, 0, 0]) # frequency of measurement = 148.760010655GHz
  
  setjy(vis = 'uid___A002_X867766_Xa7.ms.split',
    field = '1', # source name = J1517-243
    spw = '2', # center frequency of spw = 147.11975309GHz
    standard = 'manual',
    fluxdensity = [1.08134701397, 0, 0, 0]) # frequency of measurement = 148.760010655GHz
  
  setjy(vis = 'uid___A002_X867766_Xa7.ms.split',
    field = '1', # source name = J1517-243
    spw = '3', # center frequency of spw = 146.044604673GHz
    standard = 'manual',
    fluxdensity = [1.08134701397, 0, 0, 0]) # frequency of measurement = 148.760010655GHz
  
  

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

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

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

  os.system('rm -rf uid___A002_X867766_Xa7.ms.split.ap_pre_bandpass') 
  
  gaincal(vis = 'uid___A002_X867766_Xa7.ms.split',
    caltable = 'uid___A002_X867766_Xa7.ms.split.ap_pre_bandpass',
    field = '0', # J1700-2610
    spw = '0:1536~2304,1:1536~2304,2:1536~2304,3:1536~2304',
    solint = 'int',
    refant = 'DA59',
    calmode = 'p')
  
  if applyonly != True: es.checkCalTable('uid___A002_X867766_Xa7.ms.split.ap_pre_bandpass', msName='uid___A002_X867766_Xa7.ms.split', interactive=False) 
  
  os.system('rm -rf uid___A002_X867766_Xa7.ms.split.bandpass') 
  bandpass(vis = 'uid___A002_X867766_Xa7.ms.split',
    caltable = 'uid___A002_X867766_Xa7.ms.split.bandpass',
    field = '0', # J1700-2610
    solint = 'inf',
    combine = 'scan',
    refant = 'DA59',
    solnorm = True,
    bandtype = 'B',
    gaintable = 'uid___A002_X867766_Xa7.ms.split.ap_pre_bandpass')
  
  os.system('rm -rf uid___A002_X867766_Xa7.ms.split.bandpass_smooth20ch') 
  
  bandpass(vis = 'uid___A002_X867766_Xa7.ms.split',
    caltable = 'uid___A002_X867766_Xa7.ms.split.bandpass_smooth20ch',
    field = '0', # J1700-2610
    solint = 'inf,20ch',
    combine = 'scan',
    refant = 'DA59',
    solnorm = True,
    bandtype = 'B',
    gaintable = 'uid___A002_X867766_Xa7.ms.split.ap_pre_bandpass')
  
  if applyonly != True: es.checkCalTable('uid___A002_X867766_Xa7.ms.split.bandpass_smooth20ch', msName='uid___A002_X867766_Xa7.ms.split', interactive=False) 
  
  if applyonly != True: es.checkCalTable('uid___A002_X867766_Xa7.ms.split.bandpass', msName='uid___A002_X867766_Xa7.ms.split', interactive=False) 
  

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

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

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

  os.system('rm -rf uid___A002_X867766_Xa7.ms.split.phase_int') 
  gaincal(vis = 'uid___A002_X867766_Xa7.ms.split',
    caltable = 'uid___A002_X867766_Xa7.ms.split.phase_int',
    field = '0~2', # J1700-2610,J1517-243,J1625-2527
    solint = 'int',
    refant = 'DA59',
    gaintype = 'G',
    calmode = 'p',
    gaintable = 'uid___A002_X867766_Xa7.ms.split.bandpass_smooth20ch')
  
  if applyonly != True: es.checkCalTable('uid___A002_X867766_Xa7.ms.split.phase_int', msName='uid___A002_X867766_Xa7.ms.split', interactive=False) 
  
  os.system('rm -rf uid___A002_X867766_Xa7.ms.split.ampli_inf') 
  gaincal(vis = 'uid___A002_X867766_Xa7.ms.split',
    caltable = 'uid___A002_X867766_Xa7.ms.split.ampli_inf',
    field = '0~2', # J1700-2610,J1517-243,J1625-2527
    solint = 'inf',
    refant = 'DA59',
    gaintype = 'T',
    calmode = 'a',
    gaintable = ['uid___A002_X867766_Xa7.ms.split.bandpass_smooth20ch', 'uid___A002_X867766_Xa7.ms.split.phase_int'])
  
  if applyonly != True: es.checkCalTable('uid___A002_X867766_Xa7.ms.split.ampli_inf', msName='uid___A002_X867766_Xa7.ms.split', interactive=False) 
  
  os.system('rm -rf uid___A002_X867766_Xa7.ms.split.flux_inf') 
  os.system('rm -rf uid___A002_X867766_Xa7.ms.split.fluxscale') 
  mylogfile = casalog.logfile()
  casalog.setlogfile('uid___A002_X867766_Xa7.ms.split.fluxscale')
  
  fluxscaleDict = fluxscale(vis = 'uid___A002_X867766_Xa7.ms.split',
    caltable = 'uid___A002_X867766_Xa7.ms.split.ampli_inf',
    fluxtable = 'uid___A002_X867766_Xa7.ms.split.flux_inf',
    reference = '1') # J1517-243
  
  casalog.setlogfile(mylogfile)
  
  if applyonly != True: es.fluxscale2(caltable = 'uid___A002_X867766_Xa7.ms.split.ampli_inf', removeOutliers=True, msName='uid___A002_X867766_Xa7.ms', writeToFile=True, preavg=10000)
  
  os.system('rm -rf uid___A002_X867766_Xa7.ms.split.phase_inf') 
  gaincal(vis = 'uid___A002_X867766_Xa7.ms.split',
    caltable = 'uid___A002_X867766_Xa7.ms.split.phase_inf',
    field = '0~2', # J1700-2610,J1517-243,J1625-2527
    solint = 'inf',
    refant = 'DA59',
    gaintype = 'G',
    calmode = 'p',
    gaintable = 'uid___A002_X867766_Xa7.ms.split.bandpass_smooth20ch')
  
  if applyonly != True: es.checkCalTable('uid___A002_X867766_Xa7.ms.split.phase_inf', msName='uid___A002_X867766_Xa7.ms.split', interactive=False) 
  

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

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

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

  for i in ['0', '1']: # J1700-2610,J1517-243
    applycal(vis = 'uid___A002_X867766_Xa7.ms.split',
      field = i,
      gaintable = ['uid___A002_X867766_Xa7.ms.split.bandpass_smooth20ch', 'uid___A002_X867766_Xa7.ms.split.phase_int', 'uid___A002_X867766_Xa7.ms.split.flux_inf'],
      gainfield = ['', i, i],
      interp = 'linear,linear',
      calwt = F,
      flagbackup = F)
  
  applycal(vis = 'uid___A002_X867766_Xa7.ms.split',
    field = '2,3', # IRAS16293-2422
    gaintable = ['uid___A002_X867766_Xa7.ms.split.bandpass_smooth20ch', 'uid___A002_X867766_Xa7.ms.split.phase_inf', 'uid___A002_X867766_Xa7.ms.split.flux_inf'],
    gainfield = ['', '2', '2'], # J1625-2527
    interp = 'linear,linear',
    calwt = F,
    flagbackup = F)
  

# Split out corrected column
mystep = 18
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  os.system('rm -rf uid___A002_X867766_Xa7.ms.split.cal') 
  split(vis = 'uid___A002_X867766_Xa7.ms.split',
    outputvis = 'uid___A002_X867766_Xa7.ms.split.cal',
    datacolumn = 'corrected',
    keepflags = T)