##############################################################################
# Data reduction script for VY CMa Band 7 325-GHz 12-m data: Part 1: Calibration
# VYCMa_Band7_325GHz_Calibration.py
# Tested using CASA Version 4.2.1 (r29048)
#
# Datasets used:
#
# NOTE - THIS DATASET DOES NOT HAVE ITS OWN BANDPASS CALIBRATOR AND SHOULD
# BE USED WITH CAUTION. IN THE PRESENT SCRIPT THE BANDPASS CALIBRATOR IS
# COPIED FROM THE NEXT DATASET (X1ba1). WHILE THIS PROVIDES A SATISFACTORY RESULT
# FOR THE BANDPASS AND AN IMPROVED SNR, SOME USERS MAY WANT TO OMIT THIS DATASET.
#
# uid___A002_X6cdf83_X1488 - observed Aug 16 2013 with 16 antennas
# Flux scale calibrator Ganymede*
# Bandpass calibrator (copy from next dataset) 
# Phase reference source J0648-3044 
# Target Vy Cma
#
# uid___A002_X6d0f96_X1ba1 - observed Aug 16 2013 with 20 antennas
# Flux scale calibrator J0522-364  
# Bandpass calibrator J0522-364
# Phase reference source J0648-3044 
# Target Vy Cma
#
# NOTE - THIS DATASET DOES NOT HAVE ITS OWN BANDPASS CALIBRATOR AND SHOULD
# BE USED WITH CAUTION. IN THE PRESENT SCRIPT THE BANDPASS CALIBRATOR IS
# COPIED FROM THE PREVIOUS DATASET (X1ba1). WHILE THIS PROVIDES A SATISFACTORY RESULT
# FOR THE BANDPASS AND AN IMPROVED SNR, SOME USERS MAY WANT TO OMIT THIS DATASET.
#
# uid___A002_X6d0f96_X1de2 - observed Aug 17 2013 with 20 antennas
# Flux scale calibrator  Pallas*
# Bandpass calibrator    (copy from previous dataset)
# Phase reference source J0648-3044
# Target                 Vy Cma
#
# * This part of the spectrum is strongly affected by tropospheric
# * water so the most accurate flux scale is in fact obtained by using
# * the values and spectral indices measured at 321-GHz for J0648-3044
# * and J0522-364
#
##############################################################################
"""
See accompanying README file for details of the necessary input files
and important comments on the data.

Prior to running this script you will need to obtain the Analysis Utilities package.
See http://casaguides.nrao.edu/index.php?title=Analysis_Utilities for details and 
for how to download Analysis Utilities.

"""

#----------------------------------------------------------------------------------
#----- Optional Steps -------------------------------------------------------------
#----------------------------------------------------------------------------------

#----- Set this option to true if you want to make diagnostic plots along the way.
#----- Note that this will slow down the reduction somewhat.
#----- Default is no plots (you can make them later anyway).
calplots=T

#----------------------------------------------------------------------------------
#----- Some setup steps -----------------------------------------------------------
#----------------------------------------------------------------------------------

version = casalog.version()
print "You are using " + version
if (version < '4.2.1'):
    print "\033[91m YOUR VERSION OF CASA IS TOO OLD FOR THIS GUIDE."
    print "\033[91m PLEASE UPDATE IT BEFORE PROCEEDING."
else:
    print "Your version of CASA is appropriate for this guide."

# Repeat imports if you exit and restart CASA if using copy and paste
#import re
#import os

thesteps=[]
step_title = {0: 'Import of the ASDMs, listing the data and plotting antenna table',
              1: 'Flag autocorrelations and pointing scans; back up flags',
              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: 'Plot raw data',
              6: 'Application of the WVR, Tsys and antpos cal tables',
              7: 'Split out science SPWs and time average',
              8: 'Listobs, clear pointing table, and save original flags',
              9: 'Plot data with instrumental calibration',
              10: 'Initial flagging',
              11: 'Bandpass calibration',
              12: 'Putting a model for the flux calibrator (not used in this script)',
              13: 'Calibrate the flux scale',
              14: 'Gain calibration',
              15: 'Application of the bandpass and gain cal tables',
              16: 'Plot calibrated phase ref data',
              17: 'Split out the science data'}

# 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.

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

######################################################

# If executing this script by copying and pasting, and you restart
# CASA, copy the definition of basename before executing steps 0-7
#
# IF YOU WISH TO ONLY USE uid___A002_X6d0f96_X1ba1 WHICH HAS ITS OWN
# BANDPASS CALIBRATOR (SEE README EXPLANATION) USE INSTEAD: 
#basename=['uid___A002_X6d0f96_X1ba1']
#
basename=['uid___A002_X6cdf83_X1488','uid___A002_X6d0f96_X1ba1','uid___A002_X6d0f96_X1de2']
# Copy the definition of refant for steps 5-15
refantenna='DV15'

# Import of the ASDMs, listing the data and plotting antenna table
mystep = 0
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  for name in basename:
      os.system('rm -rf '+name+'.ms')
      importasdm(asdm=name, vis=name+'.ms', process_flags=F, 
                 asis='Antenna Station CalWVR Receiver Source')

# Do listobs and write info for each dataset to a text file
  for name in basename:
      os.system('rm '+name+'_listobs.txt')
      listobs(vis=name+'.ms', listfile=name+'_listobs.txt', verbose=True)
      print "<< Printed listobs output to "+name+'.listobs.txt'

# Plot antenna configuration to a file
  for name in basename:
      plotants(name+'.ms', figfile='plotants_'+name+'.png')

  print 'Choose refant - with lots of short baselines but not so close other antennas to\nthat it will be shadowed. Definition is at top of script, change if desired.'

#################################
# Flag autocorrelations and pointing scans; back up flags
mystep = 1 
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

# Remove all previous flags
  for name in basename:
      flagdata(vis=name+'.ms',mode='unflag', flagbackup=F)

  for name in basename:
      flagdata(vis = name+'.ms',
               mode = 'manual',
               spw = '1~16',
               autocorr = T,
               flagbackup = F)
  
  for name in basename:
      flagdata(vis = name+'.ms',
               mode = 'manual',
               intent = '*POINTING*,*SIDEBAND_RATIO*,*ATMOSPHERE*',
               flagbackup = F)

# Store current flagging state
  for name in basename:
      flagmanager(vis = name+'.ms', mode = 'save',
                  versionname = 'Apriori')

# To restore a priori flags:
#for name in basename:
#	flagmanager(vis = name+'.ms', mode = 'restore', versionname = 'Apriori')	

###############
# 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]
  
  for name in basename:
      os.system('rm -rf ' +name+'.ms.wvr') 
      mylogfile = casalog.logfile()
      casalog.setlogfile(name+'.ms.wvrgcal')
      wvrgcal(vis = name+'.ms',
              caltable = name+'.ms.wvr',
              toffset = 0,
              tie = ['Vy Cma,J0648-3044'],
              statsource = 'Vy Cma')
      casalog.setlogfile(mylogfile)
  
  for name in basename:
      os.system('rm -rf ' +name+'.ms.wvr.smooth') 
      smoothcal(vis = name+'.ms',
                tablein = name+'.ms.wvr',
                caltable = name+'.ms.wvr.smooth',
                smoothtype = 'mean',
                smoothtime = 6.048)
  
  if(calplots):
      print 'Please ignore any messages arising from data sets with <17 antennas'
      for name in basename:
          os.system('rm -rf '+name+'.ms.wvr.smooth*.png')
          for ant in range(3):
              plotcal(caltable = name+'.ms.wvr.smooth',
                      xaxis = 'time', yaxis = 'phase',
                      plotrange = [0,0,-180,180],
                      antenna=str(ant*8)+'~'+str(ant*8+7),
                      iteration = 'antenna',
                      markersize=3.0,
                      figfile=name+'.ms.wvr.smooth_ant'+str(ant*8)+'-'+str(ant*8+7)+'.png',
                      subplot = 421,
                      fontsize=7.0)


# Flag WVR autocorrelations
  for name in basename:
      flagdata(vis = name+'.ms',
               mode = 'manual',
               spw = '0',
               autocorr = T,
               flagbackup = F)
  
#####################################	
# 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]
 
  for name in basename:
      os.system('rm -rf '+name+'.ms.tsys') 
      gencal(vis = name+'.ms',
             caltable = name+'.ms.tsys',
             caltype = 'tsys')

# replace bad data on DV15 pol Y for uid___A002_X6d0f96_X1de2.ms
  calTable = 'uid___A002_X6d0f96_X1de2.ms.tsys'
  antenna = 11 # Antenna id
  spwsToFix = [9] # Original spw ids
  frompol = 'X'
  topol = 'Y'
  scaleFactor = 1.0
  # do not change anything below this point
  pols = ['X', 'Y']  
  frompol = pols.index(frompol)
  topol = pols.index(topol)
  tb.open(calTable,nomodify=False)
  spws = tb.getcol('SPECTRAL_WINDOW_ID')
  antennas = tb.getcol('ANTENNA1')
  replaced = 0
  for i in xrange(len(antennas)):
      if (antennas[i] == antenna and (spws[i] in spwsToFix)):
          gain = tb.getcell('FPARAM',i)
          gain[topol] = gain[frompol] * scaleFactor
          tb.putcell('FPARAM',i,gain)
          gainerr = tb.getcell('PARAMERR',i)
          gainerr[topol] = gainerr[frompol] * scaleFactor
          tb.putcell('PARAMERR',i,gainerr)
          gainflag = tb.getcell('FLAG',i)
          gainflag[topol] = gainflag[frompol]
          tb.putcell('FLAG',i,gainflag)
          gainsnr = tb.getcell('SNR',i)
          gainsnr[topol] = gainsnr[frompol]
          tb.putcell('SNR',i,gainsnr)
          replaced += 1
  tb.close()
  
  if(calplots):
      for name in basename:
          os.system('rm -rf '+name+'.ms.tsys*.png')
          for s in (9,11):
              for ant in range(3):
                  plotcal(caltable = name+'.ms.tsys',
                          xaxis = 'freq', yaxis = 'tsys',
                          spw=str(s)+':8~120',
                          markersize=3.0,
                          antenna=str(ant*8)+'~'+str(ant*8+7),
                          iteration = 'antenna',
                          figfile=name+'.ms.tsys_spw'+str(s)+'_ant'+str(ant*8)+'-'+str(ant*8+7)+'.png',
                          subplot = 421,
                          fontsize=7.0)
 

######################################
# 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]

  for name in basename:    
      os.system('rm -rf '+ name+'.ms.antpos') 

  
  os.system('rm -rf uid___A002_X6cdf83_X1488.ms.antpos') 
  gencal(vis = 'uid___A002_X6cdf83_X1488.ms',
    caltable = 'uid___A002_X6cdf83_X1488.ms.antpos',
    caltype = 'antpos',
    antenna = 'DA46,DA52,DV03,DV04,DV05,DV07,DV11,DV13,DV14,DV15,DV17,DV18,DV19,DV22,DV24,DV25',
    parameter = [0.000006, 0.000228, -0.000297, 0.000068, 0.000018, 0.000008, -0.000196, -0.000061,   0.000178, -0.000132, 0.000341, -0.000190, 0.000068, 0.000267, -0.000494, -0.000524, 0.000838, 0.000626, 0.000112, -0.000116, 0.000112, -0.000032, 0.000188, -0.000261, 0.000043, -0.000173, 0.000008, 0.000037, 0.000468, -0.000484, 0.000167, -0.000474, 0.000167, -0.000047, 0.000177, -0.000151, -0.000282,-0.000304, -0.000094, -0.000185, 0.000227, -0.000283, -0.000404, 0.003222, 0.001466, -0.000296, 0.000305, 0.000940])

  os.system('rm -rf uid___A002_X6d0f96_X1ba1.ms.antpos') 
  gencal(vis = 'uid___A002_X6d0f96_X1ba1.ms',
    caltable = 'uid___A002_X6d0f96_X1ba1.ms.antpos',
    caltype = 'antpos',
    antenna = 'DA41,DA46,DA50,DA56,DV03,DV04,DV05,DV07,DV11,DV13,DV14,DV15,DV17,DV18,DV19,DV21,DV22,DV24,DV25',
   parameter = [-0.018216, 0.067084, 0.060448, 0.000006, 0.000228, -0.000297, 0.000068, 0.000018, 0.000008, 0.000065, 0.000035, 0.000115, -0.000196, -0.000061, 0.000178, -0.000132, 0.000341, -0.000190, 0.000068, 0.000267, -0.000494, -0.000524, 0.000838, 0.000626, 0.000112, -0.000116, 0.000112, -0.000032, 0.000188, -0.000261, 0.000043, -0.000173, 0.000008, 0.000037, 0.000468, -0.000484, 0.000167, -0.000474, 0.000167, -0.000047, 0.000177, -0.000151, -0.000282, -0.000304, -0.000094, -0.000185, 0.000023, 0.000032, 0.000009, 0.000227, -0.000283, -0.000404, 0.003222, 0.001466, -0.000296, 0.000305, 0.000940])
   
  os.system('rm -rf uid___A002_X6d0f96_X1de2.ms.antpos') 
  gencal(vis = 'uid___A002_X6d0f96_X1de2.ms',
    caltable = 'uid___A002_X6d0f96_X1de2.ms.antpos',
    caltype = 'antpos',
    antenna = 'DA41,DA46,DA50,DA56,DV03,DV04,DV05,DV07,DV11,DV13,DV14,DV15,DV17,DV18,DV19,DV21,DV22,DV24,DV25',
    parameter = [-0.018216, 0.067084, 0.060448, 0.000006, 0.000228, -0.000297, 0.000068, 0.000018, 0.000008, 0.000065, 0.000035, 0.000115, -0.000196, -0.000061, 0.000178, -0.000132, 0.000341, -0.000190, 0.000068, 0.000267, -0.000494, -0.000524, 0.000838, 0.000626, 0.000112, -0.000116, 0.000112, -0.000032, 0.000188, -0.000261, 0.000043, -0.000173, 0.000008, 0.000037, 0.000468, -0.000484, 0.000167, -0.000474, 0.000167, -0.000047, 0.000177, -0.000151, -0.000282, -0.000304, -0.000094, -0.000185, 0.000023, 0.000032, 0.000009, 0.000227, -0.000283, -0.000404, 0.003222, 0.001466, -0.000296, 0.000305, 0.000940])

#######
# Plot visibilities of bandpass calibrator before applying instrumental calibration
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_X6d0f96_X1ba1.ms.rawBPsource_*.png')
    plotms(vis='uid___A002_X6d0f96_X1ba1.ms', 
           field='J0522-364',
           xaxis='frequency', yaxis='amp',
           selectdata=T, spw='13,15', 
           avgtime='1e8',avgscan=T,
           coloraxis='baseline',
           antenna=refantenna+'&*',
           ydatacolumn='data',
           showgui=F,
           plotfile='uid___A002_X6d0f96_X1ba1.ms.rawBPsource_amp.png')

    plotms(vis='uid___A002_X6d0f96_X1ba1.ms', 
           field='J0522-364',
           xaxis='frequency', yaxis='phase',
           selectdata=T, spw='13,15',
           avgtime='1e8',avgscan=T,
           coloraxis='baseline',
           antenna=refantenna+'&*',
           ydatacolumn='data',
           showgui=F,
           plotfile='uid___A002_X6d0f96_X1ba1.ms.rawBPsource_phase.png')

####################################
# Apply WVR, Tsys and antenna position corrections 
mystep = 6
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 = basename[0]+'.ms', tsystable = basename[0]+'.ms.tsys')
    
    for field in ['Ganymede','J0648-3044','Vy Cma']:
        name='uid___A002_X6cdf83_X1488'
        applycal(vis='uid___A002_X6cdf83_X1488.ms',
                 spw='13,15',
                 field=field,
                 gainfield=[field,'',''],
                 interp=['linear,linear'],
                 gaintable=[name+'.ms.tsys',name+'.ms.wvr.smooth',name+'.ms.antpos'],
                 flagbackup=F, calwt=T, 
                 spwmap = [tsysmap,[],[]])

    for field in ['J0522-364','J0648-3044','Vy Cma']:
        name='uid___A002_X6d0f96_X1ba1'
        applycal(vis = 'uid___A002_X6d0f96_X1ba1.ms',
                 spw='13,15',
                 field=field,
                 gainfield=[field,'',''],
                 interp=['linear,linear'],
                 gaintable=[name+'.ms.tsys',name+'.ms.wvr.smooth',name+'.ms.antpos'],
                 flagbackup=F, calwt=T, 
                 spwmap = [tsysmap,[],[]])

    for field in ['Pallas','J0648-3044','Vy Cma']:
        name='uid___A002_X6d0f96_X1de2'
        applycal(vis = 'uid___A002_X6d0f96_X1de2.ms',
                 spw='13,15',
                 field=field,
                 gainfield=[field,'',''],
                 interp=['linear,linear'],
                 gaintable=[name+'.ms.tsys',name+'.ms.wvr.smooth',name+'.ms.antpos'],
                 flagbackup=F, calwt=T, 
                 spwmap = [tsysmap,[],[]])

###############################################
# Split out science spw with time averaging
mystep = 7
if(mystep in thesteps):
    casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
    print 'Step ', mystep, step_title[mystep]

    for name in basename:
        os.system('rm -rf '+ name+'.ms.split') 
        split(vis = name+'.ms',
              outputvis = name+'.ms.split',
              datacolumn = 'corrected',
              spw = '13,15',
              timebin = '6.048s',
              keepflags = T)
  
#########################################
# Listobs, clear pointing table, and save original flags
mystep=8
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  for name in basename:
      os.system('rm -rf '+name+'.ms.split.listobs.txt')
      listobs(vis = name+'.ms.split',
              listfile = name+'.ms.split.listobs.txt')

# Remove POINTING tables
  for name in basename:
      tb.open(name+'.ms.split/POINTING', nomodify = False)
      a = tb.rownumbers()
      tb.removerows(a)
      tb.close()

# Store pre-flagging state
  for name in basename:
      flagmanager(vis = name+'.ms.split',
                  mode = 'save',
                  versionname = 'Pre-flagging')

#######
# Plot visibilities of bandpass calibration source after applying instrumental calibration
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_X6d0f96_X1ba1.ms.split.corrBPsource_*.png')
    plotms(vis='uid___A002_X6d0f96_X1ba1.ms.split', 
           field='J0522-364',
           xaxis='frequency', yaxis='amp',
           selectdata=T, spw='0,1', 
           avgtime='1e8',avgscan=T,
           coloraxis='baseline',
           antenna=refantenna+'&*',
           showgui=F,
           plotfile='uid___A002_X6d0f96_X1ba1.ms.split.corrBPsource_amp.png')
    plotms(vis='uid___A002_X6d0f96_X1ba1.ms.split', 
           field='J0522-364',
           xaxis='frequency', yaxis='phase',
           selectdata=T, spw='0,1',
           avgtime='1e8',avgscan=T,
           coloraxis='baseline',
           antenna=refantenna+'&*',
           showgui=F,
           plotfile='uid___A002_X6d0f96_X1ba1.ms.split.corrBPsource_phase.png')

    print 'You should see that the scatter in the data is slightly less.'
#, but a few bad\badpoints appear. These are included in the flagging commands below, bu#t you can\ninspect manually in plotms and modify flagging commands if wanted.'

#################
# Flagging shadowed data and other bad data and back up flags
mystep = 10
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  for name in basename:
      flagdata(vis = name+'.ms.split',
               mode = 'shadow',
               flagbackup = F)

# First channel is bad
  for name in basename:
      flagdata(vis=name+'.ms.split',
               mode='manual',
               spw='0:0,1:0',
               flagbackup = F)

  for name in basename:
      flagmanager(vis = name+'.ms.split',
                  mode = 'save',
                  versionname = 'BeforeBandpassCalibration')

 
#################################
# Bandpass calibration
# Only one dataset has the required source. See caveats in Readme file.
# First perform time-dependent calibration, averaging the central channels
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 cal-uid___A002_X6d0f96_X1ba1-BPint.Gap')
  gaincal(vis='uid___A002_X6d0f96_X1ba1.ms.split',
          caltable='cal-'+'uid___A002_X6d0f96_X1ba1-BPint.Gap',
          spw='*:1536~2304',field='J0522-364',
          selectdata=F, solint='int',
          refant=refantenna, calmode='ap')
  
  if(calplots): # make some plots if calplots=True     
      os.system('rm -rf cal-uid___A002_X6d0f96_X1ba1-phase-time-?-BPint-Gap.png')
      os.system('rm -rf cal-uid___A002_X6d0f96_X1ba1-amp-time-?-BPint-Gap.png')
      plotcal(caltable='cal-uid___A002_X6d0f96_X1ba1-BPint.Gap',
              xaxis = 'time', yaxis = 'phase',poln='X',
              plotsymbol='o', plotrange = [0,0,-180,180],
              iteration = 'spw',
              figfile='cal-uid___A002_X6d0f96_X1ba1-phase-time-X-BPint-Gap.png',
              subplot = 211)
      plotcal(caltable='cal-uid___A002_X6d0f96_X1ba1-BPint.Gap',
              xaxis = 'time', yaxis = 'phase',poln='Y',
              plotsymbol='o', plotrange = [0,0,-180,180],
              iteration = 'spw',
              figfile='cal-uid___A002_X6d0f96_X1ba1-phase-time-Y-BPint-Gap.png',
              subplot = 211)
      plotcal(caltable='cal-'+'uid___A002_X6d0f96_X1ba1-BPint.Gap',
              xaxis = 'time', yaxis = 'amp',poln='X',
              plotsymbol='o',
              iteration = 'spw',
              figfile='cal-uid___A002_X6d0f96_X1ba1-amp-time-X-BPint-Gap.png',
              subplot = 211)
      plotcal(caltable='cal-'+'uid___A002_X6d0f96_X1ba1-BPint.Gap',
              xaxis = 'time', yaxis = 'amp',poln='Y',
              plotsymbol='o', 
              iteration = 'spw',
              figfile='cal-uid___A002_X6d0f96_X1ba1-amp-time-Y-BPint-Gap.png',
              subplot = 211)
      
      os.system('rm -rf cal-uid___A002_X6d0f96_X1ba1.B1')
      bandpass(vis='uid___A002_X6d0f96_X1ba1.ms.split',
               caltable='cal-uid___A002_X6d0f96_X1ba1.B1',
               field='J0522-364', 
               bandtype='B', 
               fillgaps=10, solnorm = T, 
               solint='inf,20ch', 
               refant=refantenna,
               gaintable='cal-uid___A002_X6d0f96_X1ba1-BPint.Gap')


  if(calplots): # make some plots if calplots=True
      os.system('rm -rf cal-uid___A002_X6d0f96_X1ba1*-B1.png')
      for spw in ['0','1']:
          for ant in range(3):
              plotcal(caltable = 'cal-uid___A002_X6d0f96_X1ba1.B1',
                      xaxis='freq', yaxis='phase', spw=spw, 
                      antenna=str(ant*8)+'~'+str(ant*8+7),
                      iteration='antenna',subplot=421, 
                      fontsize=7.0,
                      plotrange = [0,0,-180,180],
                      figfile='cal-uid___A002_X6d0f96_X1ba1_-phase_spw'+spw+'_ant'+str(ant*8)+'-'+str(ant*8+7)+'-B1.png')
              plotcal(caltable = 'cal-uid___A002_X6d0f96_X1ba1.B1',
                      xaxis='freq', yaxis='amp', spw=spw, 
                      antenna=str(ant*8)+'~'+str(ant*8+7),
                      iteration='antenna',subplot=421, 
                      fontsize=7.0,
                      figfile='cal-uid___A002_X6d0f96_X1ba1_-amp_spw'+spw+'_ant'+str(ant*8)+'-'+str(ant*8+7)+'-B1.png')


#################
# Putting in model for flux calibrators Pallas and Ganymede
# This step can be omitted 
mystep = 12
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  setjy(vis = 'uid___A002_X6cdf83_X1488.ms.split',
    field = 'Ganymede', 
    spw = '0,1',
    standard = 'Butler-JPL-Horizons 2012')
  os.system('rm -rf uid___A002_X6cdf83_X1488.ms.split.setjy.Pallas.png') 
  plotms(vis = 'uid___A002_X6cdf83_X1488.ms.split',
         xaxis = 'uvdist',
         yaxis = 'amp',
         ydatacolumn = 'model',
         field = 'Pallas',
         spw = '0,1',
         avgchannel = '9999',
         coloraxis = 'spw',
         plotfile = 'uid___A002_X6cdf83_X1488.ms.split.setjy.Ganymede.png')
 
  setjy(vis = 'uid___A002_X6d0f96_X1de2.ms.split',
    field = 'Pallas', 
    spw = '0,1',
    standard = 'Butler-JPL-Horizons 2012')
  os.system('rm -rf uid___A002_X6d0f96_X1de2.ms.split.setjy.Pallas.png') 
  plotms(vis = 'uid___A002_X6d0f96_X1de2.ms.split',
         xaxis = 'uvdist',
         yaxis = 'amp',
         ydatacolumn = 'model',
         field = 'Pallas',
         spw = '0,1',
         avgchannel = '9999',
         coloraxis = 'spw',
         plotfile = 'uid___A002_X6d0f96_X1de2.ms.split.setjy.Pallas.png')

#################################
# Setting the flux scale
# We use the flux scale derived at 321 GHz for J0522-364 and J0648-3044
# These are the values I got
# Fitted spectrum for J0522-364 with fitorder=1: 
# Flux density = 4.77004 +/- 0.00314719 (freq=316.093 GHz) 
# spidx=-0.53027 +/- 0.0397935
# Fitted spectrum for J0648-3044 with fitorder=1: 
# Flux density = 0.440988 +/- 0.000683724 (freq=316.093 GHz) 
# spidx=-0.776429 +/- 0.0934349

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

# Remove any previous models for QSOs
  for name in basename:
      delmod(vis=name+'.ms.split', field='J*')

  J0522_flux=4.750479
  J0522_rfreq=316.093
  J0522_spidx=-0.53027

  J0648_flux=0.440988
  J0648_rfreq=316.093
  J0648_spidx=-0.776429

  setjy(vis='uid___A002_X6d0f96_X1ba1.ms.split',
        field='J0522-364',
        standard='manual',
        fluxdensity=[J0522_flux,0,0,0],
        spix=J0522_spidx,
        reffreq=str(J0522_rfreq)+'GHz')      

  for name in basename:
      setjy(vis=name+'.ms.split',
            field='J0648-3044',
            standard='manual',
            fluxdensity=[J0648_flux,0,0,0],
            spix=J0648_spidx,
            reffreq=str(J0648_rfreq)+'GHz')

#################
# Time-dependent calibration 
mystep = 14
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

# phase calibration per integration to apply to cal sources
  for name in basename:
      os.system('rm -rf cal-'+name+'.Gp_int')
      gaincal(vis=name+'.ms.split',
              caltable = 'cal-'+name+'.Gp_int',
              field='J*',
              solint='int',
              refant=refantenna,
              calmode='p',
              gaintable = 'cal-uid___A002_X6d0f96_X1ba1.B1')

  if(calplots): # make some plots if calplots=True
      for name in basename:
          os.system('rm -rf cal-'+name+'-phase-time-?-Gp_int.png')
          plotcal(caltable='cal-'+name+'.Gp_int',
                  xaxis = 'time', yaxis = 'phase',poln='X',
                  plotsymbol='o', plotrange = [0,0,-180,180],
                  iteration = 'spw',
                  figfile='cal-'+name+'-phase-time-X-Gp_int.png',
                  subplot = 211)
          plotcal(caltable='cal-'+name+'.Gp_int',
                  xaxis = 'time', yaxis = 'phase',poln='Y',
                  plotsymbol='o', plotrange = [0,0,-180,180],
                  iteration = 'spw',
                  figfile='cal-'+name+'-phase-time-Y-Gp_int.png',
                  subplot = 211)
 
# phase calibration per scan for phase-ref to apply to target
  for name in basename:
      os.system('rm -rf cal-'+name+'.Gp_inf')
      gaincal(vis=name+'.ms.split',
              caltable = 'cal-'+name+'.Gp_inf',
              field='J0648-3044',
              solint='inf',
              refant=refantenna,
              calmode='p',
              gaintable = 'cal-uid___A002_X6d0f96_X1ba1.B1')

# amp cal per scan
  for name in basename:
      os.system('rm -rf cal-'+name+'.Gap_inf')
      gaincal(vis=name+'.ms.split',
              caltable = 'cal-'+name+'.Gap_inf',
              field='J*',
              solint='inf',
              refant=refantenna,
              gaintype='T',
              calmode='ap',
              gaintable = ['cal-uid___A002_X6d0f96_X1ba1.B1','cal-'+name+'.Gp_int'])

  if(calplots): # make some plots if calplots=True
      for name in basename:
          os.system('rm -rf cal-'+name+'-amp-time-Gap_inf.png')
          plotcal(caltable='cal-'+name+'.Gap_inf',
                  xaxis = 'time', yaxis = 'amp',
                  plotsymbol='o',
                  iteration = 'spw',
                  figfile='cal-'+name+'-amp-time-Gap_inf.png',
                  subplot = 211)
 
#################
# Apply calibration
mystep = 15
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  applycal(vis='uid___A002_X6d0f96_X1ba1.ms.split',
           field='J0522-364',
           gaintable=['cal-uid___A002_X6d0f96_X1ba1.B1','cal-uid___A002_X6d0f96_X1ba1.Gp_int','cal-uid___A002_X6d0f96_X1ba1.Gap_inf'],
           gainfield=['','J0522-364','J0522-364'],
           interp='linear,linear',
           calwt=F,flagbackup=F)


  for name in basename:
      applycal(vis=name+'.ms.split',
               field='J0648-3044',
               gaintable=['cal-uid___A002_X6d0f96_X1ba1.B1','cal-'+name+'.Gp_int','cal-'+name+'.Gap_inf'],
               gainfield=['','J0648-3044','J0648-3044'],
               interp='linear,linear',
               calwt=F,flagbackup=F)

  for name in basename:
      applycal(vis=name+'.ms.split',
               field='V*',
               gaintable=['cal-uid___A002_X6d0f96_X1ba1.B1','cal-'+name+'.Gp_inf','cal-'+name+'.Gap_inf'],
               gainfield=['','J0648-3044','J0648-3044'],
               interp='linear,linear',
               calwt=F,flagbackup=F)

#######
# Plot phase-ref visibilities after applying instrumental calibration
mystep = 16
if(mystep in thesteps):
    casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
    print 'Step ', mystep, step_title[mystep]

    for name in basename:
        os.system('rm -rf '+name+'.ms.split.J0648-3044_*.png')
	plotms(vis=name+'.ms.split', 
	       field='J0648-3044',
	       xaxis='time', yaxis='amp',
	       selectdata=T, spw='0,1', 
	       avgtime='30s',avgchannel='3840',
	       coloraxis='spw',
               ydatacolumn='corrected',
               antenna=refantenna+'&*',
               showgui=F,
               plotfile=name+'.ms.split.J0648-3044_amp.png')
	plotms(vis=name+'.ms.split', 
	       field='J0648-3044',
	       xaxis='time', yaxis='phase',
	       selectdata=T, spw='0,1',
	       avgtime='30s',avgchannel='3840',
	       coloraxis='spw',
               ydatacolumn='corrected',
               antenna=refantenna+'&*',
               showgui=F,
               plotfile=name+'.ms.split.J0648-3044_phase.png')


###########################
# Split out the science data, concatenate and transform to lsrk to allow self-calibration
mystep = 17
if(mystep in thesteps):
    casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
    print 'Step ', mystep, step_title[mystep]

    vislist=[]
    for name in basename:
            os.system('rm -rf '+name+'_VYCMa_325.ms')
            mstransform(vis=name+'.ms.split',
            outputvis=name+'_VYCMa_325.ms',datacolumn='corrected',
            field='V*') 
            vislist.append(name+'_VYCMa_325.ms')

    os.system('rm -rf VYCMa_325_origspw.ms')
    concat(vis=vislist,
           concatvis='VYCMa_325_origspw.ms')
       
# Combine and regrid each group of spw with similar frequencies 
    vislist=[]
    sps=['0,2,4','1,3,5']
    for s in range(2):
        os.system('rm -rf VYCMa_325_spw'+str(s)+'.ms')
        mstransform(vis='VYCMa_325_origspw.ms',
            outputvis='VYCMa_325_spw'+str(s)+'.ms',datacolumn='data',
            spw=sps[s],
            field='V*',combinespws=T,regridms=T,
            mode='frequency',
            outframe='lsrk')
        vislist.append('VYCMa_325_spw'+str(s)+'.ms')

# Combine all regridded data
    os.system('rm -rf VYCMa_Band7_325.ms')
    concat(vis=vislist,
           concatvis='VYCMa_Band7_325.ms')
    
# Do listobs and write info to a text file
    os.system('rm VYCMa_Band7_325_listobs.txt')
    listobs(vis='VYCMa_Band7_325.ms', listfile='VYCMa_Band7_325_listobs.txt', verbose=True)
    print '<< Printed listobs output to VYCMa_Band7_325_listobs.txt'

#----------------------------------------------------------------------------------
#----- End of calibration script.
#----- To continue, see imaging script "VYCMa_Band7_325GHz_Imaging.py"
#----------------------------------------------------------------------------------