########################################################################
# Data reduction script for HD163296 Band 6:
# - Calibration script: "HD163296_Band6_Calibration.py" -
# Tested in CASA Version 3.4.0 (r19988)
########################################################################
"""
See accompanying README file for details of the necessary input files
and comments on the data

"""
#----------------------------------------------------------------------------------
#----- 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 significantly.
#----- Default is no plots (you can make them later anyway).
calplots=F

#----- Set this option to true if you want to run this script interactively to make inspection plots.
#----- You'll need to hit <enter> at various stages to continue the script.
#----- Default is true (user input is required).
#----- Set to false if you want no user interaction except clean.
interact=T 

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

version = casalog.version()
print "You are using " + version
if (int(version.split()[4][1:-1]) < 19874):
    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."

#----- Input MS names 
#-----  (the data provided have already been coverted to CASA measurement set
#-----  format (.ms) using the task importasdm)
basename=['uid___A002_X43b3b1_X47','uid___A002_X4267ef_X1b6','uid___A002_X453114_Xb3']


#----------------------------------------------------------------------------------
#----- Initial summary information ------------------------------------------------
#----------------------------------------------------------------------------------

#----- Create a summary text file for each dataset
for asdm in basename:
	os.system('rm -rf '+asdm+'.listobs.txt')
        listobs(vis=asdm+'.ms', listfile=asdm+'.listobs.txt', verbose=True)
        print "<< Printed listobs output to "+asdm+'.listobs.txt'

#----- Plot the antenna configuration for each dataset to a file
#for asdm in basename:
#        print "<< Plotting antenna configuration for dataset "+asdm+' to '+asdm+'.plotants.png'
#	os.system('rm -rf '+asdm+'.plotants.png')
#        plotants(vis=asdm+'.ms', figfile=asdm+'.plotants.png')
#        if(interact):
#		dummy_string = raw_input("Hit <Enter> to see the antenna configuration for the next data set")

#----------------------------------------------------------------------------------
#----- Inspect the data -----------------------------------------------------------
#----------------------------------------------------------------------------------

#----- Inspect amplitude versus time for each dataset
if(interact):
	for asdm in basename:
		plotms(vis=asdm+'.ms',
		       interactive=T,
		       spw='1',
		       xaxis='time', yaxis='amp',
		       correlation='XX',
		       antenna='*&*',
		       avgchannel='10000',coloraxis='field')
		dummy_string = raw_input("Once the plot has loaded, hit <Enter> to see the next data set")


#----------------------------------------------------------------------------------
#----- Generate and examine Tsys and WVR calibration tables -----------------------
#----------------------------------------------------------------------------------

#Removing existing Tsys information

for asdm in basename:
    os.system('rm -rf '+asdm+'.ms.tsys')
    print 'Creating TDM Tsys Table for '+asdm+'.ms'
    gencal(vis=asdm+'.ms',caltable=asdm+'.ms.tsys',caltype='tsys')

    if(interact):
        print "Plotting Tsys vs. time for "+asdm
        plotcal(caltable=asdm+'.ms.tsys', 
	    xaxis="time",yaxis="tsys",
	    plotsymbol=".", subplot=421,
	    antenna='0~7',
	    iteration='antenna', figfile=asdm+'.tsys_vs_time.page1.png',
	    fontsize=7.0)    
        #dummy_string = raw_input("First eight antennas for "+asdm+" . Hit <Enter> to continue.")
        plotcal(caltable=asdm+'.ms.tsys', 
	    xaxis="time",yaxis="tsys",
	    plotsymbol=".", subplot=421,
	    antenna='8~15',
	    iteration='antenna', figfile=asdm+'.tsys_vs_time.page1.png',
	    fontsize=7.0)    
        #dummy_string = raw_input("Next eight antennas for "+asdm+" . Hit <Enter> to continue.")
        plotcal(caltable=asdm+'.ms.tsys', 
	    xaxis="time",yaxis="tsys",
	    antenna='16~19',
	    spw='5:50~50',plotsymbol=".", subplot=421,
	    iteration='antenna', figfile=asdm+'.tsys_vs_time.page2.png',
	    fontsize=7.0)    
        #dummy_string = raw_input("Remaining antennas for "+asdm+" . Hit <Enter> to continue.")

#----- Generation and time averaging of the WVR cal table
#Removing existing WVR information
os.system('rm -rf *.wvrgcal')
for asdm in basename:
    wvrgcal(vis=asdm+'.ms',caltable=asdm+'.wvrgcal',toffset=-1,tie=['HD163296,J1733-130'],statsource = 'HD163296')
    smoothcal(vis=asdm+'.ms',tablein=asdm+'.wvrgcal',caltable=asdm+'.ms.wvr.smooth',
              smoothtype='mean',smoothtime=6.048)

    if(interact):
        aU.plotWVRSolutions(caltable=asdm+'.ms.wvr.smooth', spw='1',
            yrange=[-180,180],subplot=22, interactive=True,
            figfile=asdm+'.wvr.smooth.plots/'+asdm+'.ms.wvr.smooth') 


#----------------------------------------------------------------------------------
#----- Inspect the data again -----------------------------------------------------
#----------------------------------------------------------------------------------

#----- Plot spectral plot in amplitude and phase
if(interact):
    for asdm in basename:
	plotms(vis=asdm+'.ms', 
	       field='3c279',
	       xaxis='frequency', yaxis='amp',
	       selectdata=T, spw='1,3,5,7', 
	       avgtime='1e8',avgscan=T,
	       coloraxis='corr',
	       iteraxis='baseline',
	       antenna='*&*',
	       ydatacolumn='data')
	dummy_string = raw_input("Next plot for "+asdm+" . Hit <Enter> to continue.")
	plotms(vis=asdm+'.ms', 
	       field='3c279',
	       xaxis='frequency', yaxis='phase',
	       selectdata=T, spw='1,3,5,7', 
	       avgtime='1e8',avgscan=T,
	       coloraxis='corr',
	       iteraxis='baseline',
	       antenna='*&*',
	       ydatacolumn='data')
 
#----- Plot continuum plot in amplitude and phase
if(interact):
    for asdm in basename:
	plotms(vis=asdm+'.ms', 
	       field='',
	       xaxis='time', yaxis='amp',
	       selectdata=T, spw='1,3,5,7', 
	       avgchannel='3840',avgscan=F,
	       coloraxis='corr',
	       antenna='*&*',
	       ydatacolumn='data')
	dummy_string = raw_input("Next plot for "+asdm+" . Hit <Enter> to continue.")
	plotms(vis=asdm+'.ms', 
	       field='',
	       xaxis='time', yaxis='phase',
	       selectdata=T, spw='1', 
	       avgchannel='3840',avgscan=F,
	       coloraxis='corr',
	       iteraxis='baseline',
	       antenna='*&*',
	       ydatacolumn='data')


#----------------------------------------------------------------------------------
#----- Calibration of uid___A002_X43b3b1_X47.ms -----------------------------------
#----------------------------------------------------------------------------------

# Using reference antenna = DV04

asdm = 'uid___A002_X43b3b1_X47'

print "# A priori calibration"
   
fixplanets(vis = 'uid___A002_X43b3b1_X47.ms', field = '1', fixuvw = T)

# A priori flagging
tflagdata(vis = 'uid___A002_X43b3b1_X47.ms',
          mode = 'manual',autocorr = T,flagbackup = F)

# Bad Tsys - flagging DV19 XX data
tflagdata(vis = 'uid___A002_X43b3b1_X47.ms', 
          mode = 'manual',antenna = 'DV19',
          correlation = 'XX',flagbackup = F)

# Bad Tsys -  flagging DA42 spw13 data
tflagdata(vis = 'uid___A002_X43b3b1_X47.ms',
          mode = 'manual',antenna = 'DA42',
          spw = '13', flagbackup = F)
  
tflagdata(vis = 'uid___A002_X43b3b1_X47.ms',
          mode = 'manual',
          intent = '*POINTING*,*SIDEBAND_RATIO*,*ATMOSPHERE*',flagbackup = F)
  
flagcmd(vis = 'uid___A002_X43b3b1_X47.ms',
        inpmode = 'table',action = 'plot')
  
flagcmd(vis = 'uid___A002_X43b3b1_X47.ms',
        inpmode = 'table',action = 'apply')
  
 
# Generation of the antenna position cal table
  
os.system('rm -rf uid___A002_X43b3b1_X47.ms.antpos') 
gencal(vis = 'uid___A002_X43b3b1_X47.ms',
    caltable = 'uid___A002_X43b3b1_X47.ms.antpos',
    caltype = 'antpos',
    antenna = 'DV05,DV08,DV09,DV11,DA44,DV13,DA46,DA41,DV03,DA43,DA42,DV17,DV10,DV16,DV14,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 = [0.000407105675771,-6.81411070583e-05,-0.000444496072293,0.000428024806118,0.000192386926053,-7.13865255571e-05,0.000397703996595,0.000303499271338,-0.000252154165917,-0.000407067127526,0.000235328450799,-1.69882550836e-05,0.000465348935348,0.000344158009965,0.000228367319351,0.000584870703864,0.000952647723338,-0.000227043112602,0.000247000250965,-0.000508999451995,-0.000311999581754,0.00060776527971,-0.000525096431375,-0.000386487226933,0.000227385910925,-0.000120083500757,-7.74828167994e-07,0.000478540400055,0.000105861353228,3.76092994824e-05,0.000157838716115,0.000211078409778,-0.00037259206194,-7.6349849694e-05,0.000854968838704,1.63702564808e-05,0.000486001936163,7.96910315645e-05,0.000160790103968,6.20906828402e-05,0.000561427652214,0.000133509319652,0.000293836351762,0.000269524465858,-8.05776686625e-05,0.00023328837383,0.000377916222577,-0.000436685084871])
  

# Application of the WVR, Tsys and antpos cal tables
  
from recipes.almahelpers import tsysspwmap
tsysmap=tsysspwmap(vis = 'uid___A002_X43b3b1_X47.ms', tsystable = 'uid___A002_X43b3b1_X47.ms.tsys')
  
applycal(vis = 'uid___A002_X43b3b1_X47.ms',
    field = '0',
    spw = '17,19,21,23',
    gaintable = ['uid___A002_X43b3b1_X47.ms.tsys', 'uid___A002_X43b3b1_X47.ms.wvr.smooth', 'uid___A002_X43b3b1_X47.ms.antpos'],
    gainfield = ['0', '', ''],
    interp = 'linear,linear',
    spwmap = [tsysmap,'',''],
    calwt = T,
    flagbackup = F)

applycal(vis = 'uid___A002_X43b3b1_X47.ms',
    field = '1',
    spw = '17,19,21,23',
    gaintable = ['uid___A002_X43b3b1_X47.ms.tsys', 'uid___A002_X43b3b1_X47.ms.wvr.smooth', 'uid___A002_X43b3b1_X47.ms.antpos'],
    gainfield = ['1', '', ''],
    interp = 'linear,linear',
    spwmap = [tsysmap,'',''],
    calwt = T,
    flagbackup = F)

applycal(vis = 'uid___A002_X43b3b1_X47.ms',
    field = '2',
    spw = '17,19,21,23',
    gaintable = ['uid___A002_X43b3b1_X47.ms.tsys', 'uid___A002_X43b3b1_X47.ms.wvr.smooth', 'uid___A002_X43b3b1_X47.ms.antpos'],
    gainfield = ['2', '', ''],
    interp = 'linear,linear',
    spwmap = [tsysmap,'',''],
    calwt = T,
    flagbackup = F)

applycal(vis = 'uid___A002_X43b3b1_X47.ms',
    field = '3',
    spw = '17,19,21,23',
    gaintable = ['uid___A002_X43b3b1_X47.ms.tsys', 'uid___A002_X43b3b1_X47.ms.wvr.smooth', 'uid___A002_X43b3b1_X47.ms.antpos'],
    gainfield = ['3', '', ''],
    interp = 'linear,linear',
    spwmap = [tsysmap,'',''],
    calwt = T,
    flagbackup = F)


# Split out science SPWs and time average

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


print "# Calibration"

# Listobs, clear pointing table, and save original flags

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

tb.open('uid___A002_X43b3b1_X47.ms.split/POINTING', nomodify = False)
a = tb.rownumbers()
tb.removerows(a)
tb.close()

flagmanager(vis = 'uid___A002_X43b3b1_X47.ms.split',
    mode = 'save',
    versionname = 'Original')


# Initial flagging

# Flagging shadowed data
tflagdata(vis = 'uid___A002_X43b3b1_X47.ms.split',
    mode = 'shadow',
    flagbackup = F)

tflagdata(vis = 'uid___A002_X43b3b1_X47.ms.split',
    mode = 'manual', antenna='PM03',
    flagbackup = F)

tflagdata(vis = 'uid___A002_X43b3b1_X47.ms.split',
    mode = 'manual', antenna='DV10',correlation='YY', 
    flagbackup = F)


# Putting a model for the flux calibrator(s)

setjy(vis = 'uid___A002_X43b3b1_X47.ms.split',
    field = '1', # Neptune
    spw = '0,1,2,3',
    standard = 'Butler-JPL-Horizons 2010')

os.system('rm -rf uid___A002_X43b3b1_X47.ms.split.setjy.field*.png')

if(interact):
  for i in ['1']:
    plotms(vis = 'uid___A002_X43b3b1_X47.ms.split',
      xaxis = 'uvdist',
      yaxis = 'amp',
      ydatacolumn = 'model',
      field = i,
      spw = '0,1,2,3',
      avgchannel = '9999',
      coloraxis = 'spw',
      plotfile = 'uid___A002_X43b3b1_X47.ms.split.setjy.field'+i+'.png')


# Save flags before bandpass cal
flagmanager(vis = 'uid___A002_X43b3b1_X47.ms.split',
    mode = 'save',
    versionname = 'BeforeBandpassCalibration')


# Bandpass calibration

os.system('rm -rf uid___A002_X43b3b1_X47.ms.split.ap_pre_bandpass') 
gaincal(vis = 'uid___A002_X43b3b1_X47.ms.split',
    caltable = 'uid___A002_X43b3b1_X47.ms.split.ap_pre_bandpass',
    field = '0', # J1924-292
    spw = '0:1536~2304,1:1536~2304,2:1536~2304,3:1536~2304',
    solint = 'int',
    refant = 'DV04',
    calmode = 'ap')

os.system('rm -rf uid___A002_X43b3b1_X47.ms.split.bandpass') 
bandpass(vis = 'uid___A002_X43b3b1_X47.ms.split',
    caltable = 'uid___A002_X43b3b1_X47.ms.split.bandpass',
    field = '0', # J1924-292
    solint = 'inf',
    combine = 'scan',
    refant = 'DV04',
    solnorm = T,
    bandtype = 'B',
    gaintable = 'uid___A002_X43b3b1_X47.ms.split.ap_pre_bandpass')

  #aU.smoothbandpass('uid___A002_X43b3b1_X47.ms.split.bandpass') 


# Save flags before gain cal

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


# Gain calibration

clearcal('uid___A002_X43b3b1_X47.ms.split',field='J1733-130')

os.system('rm -rf uid___A002_X43b3b1_X47.ms.split.phase_int') 
gaincal(vis = 'uid___A002_X43b3b1_X47.ms.split',
    caltable = 'uid___A002_X43b3b1_X47.ms.split.phase_int',
    field = '0~2', # J1924-292,Neptune,J1733-130
    solint = 'int',
    refant = 'DV04',
    gaintype = 'G',
    calmode = 'p',
    gaintable = 'uid___A002_X43b3b1_X47.ms.split.bandpass')

os.system('rm -rf uid___A002_X43b3b1_X47.ms.split.ampli_inf') 
gaincal(vis = 'uid___A002_X43b3b1_X47.ms.split',
    caltable = 'uid___A002_X43b3b1_X47.ms.split.ampli_inf',
    field = '0~2', # J1924-292,Neptune,J1733-130
    solint = 'inf',
    refant = 'DV04',
    gaintype = 'T',
    calmode = 'ap',
    gaintable = ['uid___A002_X43b3b1_X47.ms.split.bandpass', 'uid___A002_X43b3b1_X47.ms.split.phase_int'])

os.system('rm -rf uid___A002_X43b3b1_X47.ms.split.flux_inf') 
os.system('rm -rf uid___A002_X43b3b1_X47.ms.split.fluxscale') 
casalog.setlogfile('uid___A002_X43b3b1_X47.ms.split.fluxscale')

fluxscale(vis = 'uid___A002_X43b3b1_X47.ms.split',
    caltable = 'uid___A002_X43b3b1_X47.ms.split.ampli_inf',
    fluxtable = 'uid___A002_X43b3b1_X47.ms.split.flux_inf',
    reference = '1',
    refspwmap=[0,1,3,3]) # Neptune

casalog.setlogfile('')

os.system('rm -rf uid___A002_X43b3b1_X47.ms.split.phase_inf') 
gaincal(vis = 'uid___A002_X43b3b1_X47.ms.split',
    caltable = 'uid___A002_X43b3b1_X47.ms.split.phase_inf',
    field = '0~2', # J1924-292,Neptune,J1733-130
    solint = 'inf',
    refant = 'DV04',
    gaintype = 'G',
    calmode = 'p',
    gaintable = 'uid___A002_X43b3b1_X47.ms.split.bandpass')


# Save flags before applycal cal
flagmanager(vis = 'uid___A002_X43b3b1_X47.ms.split',
    mode = 'save',
    versionname = 'BeforeApplycal')


# Application of the bandpass and gain cal tables

for i in ['0', '1']: # J1924-292,Neptune
    applycal(vis = 'uid___A002_X43b3b1_X47.ms.split',
      field = i,
      gaintable = ['uid___A002_X43b3b1_X47.ms.split.bandpass', 'uid___A002_X43b3b1_X47.ms.split.phase_int', 'uid___A002_X43b3b1_X47.ms.split.flux_inf'],
      gainfield = ['', i, i],
      interp = 'linear,linear',
      calwt = F,
      flagbackup = F)

applycal(vis = 'uid___A002_X43b3b1_X47.ms.split',
    field = '2,3', # HD163296
    gaintable = ['uid___A002_X43b3b1_X47.ms.split.bandpass', 'uid___A002_X43b3b1_X47.ms.split.phase_inf', 'uid___A002_X43b3b1_X47.ms.split.flux_inf'],
    gainfield = ['', '2', '2'], # J1733-130
    interp = 'linear,linear',
    calwt = F,
    flagbackup = F)


#----------------------------------------------------------------------------------
#----- Calibration of uid___A002_X4267ef_X1b6.ms ----------------------------------
#----------------------------------------------------------------------------------

# Using reference antenna = DV04

asdm = ' uid___A002_X4267ef_X1b6'

print "# A priori calibration"


# A priori flagging
tflagdata(vis = 'uid___A002_X4267ef_X1b6.ms',
    mode = 'manual',
    autocorr = T,
    flagbackup = F)


tflagdata(vis = 'uid___A002_X4267ef_X1b6.ms',
    mode = 'manual',
    intent = '*POINTING*,*SIDEBAND_RATIO*,*ATMOSPHERE*',
    flagbackup = F)


flagcmd(vis = 'uid___A002_X4267ef_X1b6.ms',
    inpmode = 'table',
    action = 'plot')


flagcmd(vis = 'uid___A002_X4267ef_X1b6.ms',
    inpmode = 'table',
    action = 'apply')


# Generation of the antenna position cal table

os.system('rm -rf uid___A002_X4267ef_X1b6.ms.antpos') 
gencal(vis = 'uid___A002_X4267ef_X1b6.ms',
    caltable = 'uid___A002_X4267ef_X1b6.ms.antpos',
    caltype = 'antpos',
    antenna = 'DV18,DV05,CM05,DV08,DV09,DA44,DV12,DV02,DV03,DA43,DA42,DV17,DV10,DV16,DV14,DV15,DA41',
    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])
  #  parameter = [-0.000234402626516,0.000356489476784,0.000113836117399,-0.000153000877315,-0.000399853715125,-8.07882295961e-05,-1.22678466141e-05,8.66865739226e-05,-4.28566709161e-05,-0.000153605189238,9.29090249286e-05,6.73747647555e-05,-9.62091062928e-05,0.000211448146537,0.000118640166858,-3.49322730402e-05,0.000333722228151,7.52144317149e-05,-0.000133162034778,0.00023574780856,0.000212197651384,-1.03950910152e-05,-0.000460584326197,-0.000282536759635,0.000227385910925,-0.000120083500757,-7.74828167994e-07,-8.6869038668e-05,0.000257294852406,-2.23411961739e-05,-0.000309263198507,0.000119513280554,-0.000341309023937,-0.000114322037339,0.000245166722321,-0.000154021914172,-0.000136382804963,8.25043158853e-05,6.70654036784e-05,-0.000110827047308,0.000231352682506,-0.00011547857759,-8.68336887457e-05,0.000241361754262,8.97519586356e-05,-0.000113920733063,-1.20425004929e-05,-0.000206850295528,0.00060776527971,-0.000525096431375,-0.000386487226933])


# Application of the WVR, Tsys and antpos cal tables

from recipes.almahelpers import tsysspwmap
tsysmap = tsysspwmap(vis = 'uid___A002_X4267ef_X1b6.ms', tsystable = 'uid___A002_X4267ef_X1b6.ms.tsys')


applycal(vis = 'uid___A002_X4267ef_X1b6.ms',
    field = '0',
    spw = '17,19,21,23',
    gaintable = ['uid___A002_X4267ef_X1b6.ms.tsys', 'uid___A002_X4267ef_X1b6.ms.wvr.smooth', 'uid___A002_X4267ef_X1b6.ms.antpos'],
    gainfield = ['0', '', ''],
    interp = 'linear,linear',
    spwmap = [tsysmap,'',''],
    calwt = T,
    flagbackup = F)


applycal(vis = 'uid___A002_X4267ef_X1b6.ms',
    field = '1',
    spw = '17,19,21,23',
    gaintable = ['uid___A002_X4267ef_X1b6.ms.tsys', 'uid___A002_X4267ef_X1b6.ms.wvr.smooth', 'uid___A002_X4267ef_X1b6.ms.antpos'],
    gainfield = ['1', '', ''],
    interp = 'linear,linear',
    spwmap = [tsysmap,'',''],
    calwt = T,
    flagbackup = F)


applycal(vis = 'uid___A002_X4267ef_X1b6.ms',
    field = '2',
    spw = '17,19,21,23',
    gaintable = ['uid___A002_X4267ef_X1b6.ms.tsys', 'uid___A002_X4267ef_X1b6.ms.wvr.smooth', 'uid___A002_X4267ef_X1b6.ms.antpos'],
    gainfield = ['2', '', ''],
    interp = 'linear,linear',
    spwmap = [tsysmap,'',''],
    calwt = T,
    flagbackup = F)


applycal(vis = 'uid___A002_X4267ef_X1b6.ms',
    field = '3',
    spw = '17,19,21,23',
    gaintable = ['uid___A002_X4267ef_X1b6.ms.tsys', 'uid___A002_X4267ef_X1b6.ms.wvr.smooth', 'uid___A002_X4267ef_X1b6.ms.antpos'],
    gainfield = ['3', '', ''],
    interp = 'linear,linear',
    spwmap = [tsysmap,'',''],
    calwt = T,
    flagbackup = F)


# Split out science SPWs and time average

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


# Listobs, clear pointing table, and save original flags

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

tb.open('uid___A002_X4267ef_X1b6.ms.split/POINTING', nomodify = False)
a = tb.rownumbers()
tb.removerows(a)
tb.close()

flagmanager(vis = 'uid___A002_X4267ef_X1b6.ms.split',
    mode = 'save',
    versionname = 'Original')


# Initial flagging

# Flagging shadowed data
tflagdata(vis = 'uid___A002_X4267ef_X1b6.ms.split',
    mode = 'shadow',
    flagbackup = F)


# Putting a model for the flux calibrator(s)

setjy(vis = 'uid___A002_X4267ef_X1b6.ms.split',
    field = '1', # Juno
    spw = '0,1,2,3',
    standard = 'Butler-JPL-Horizons 2010')


os.system('rm -rf uid___A002_X4267ef_X1b6.ms.split.setjy.field*.png')

for i in ['1']:
    plotms(vis = 'uid___A002_X4267ef_X1b6.ms.split',
      xaxis = 'uvdist',
      yaxis = 'amp',
      ydatacolumn = 'model',
      field = i,
      spw = '0,1,2,3',
      avgchannel = '9999',
      coloraxis = 'spw',
      plotfile = 'uid___A002_X4267ef_X1b6.ms.split.setjy.field'+i+'.png')


# Save flags before bandpass cal

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


# Bandpass calibration

os.system('rm -rf uid___A002_X4267ef_X1b6.ms.split.ap_pre_bandpass') 

gaincal(vis = 'uid___A002_X4267ef_X1b6.ms.split',
    caltable = 'uid___A002_X4267ef_X1b6.ms.split.ap_pre_bandpass',
    field = '0', # J1924-292
    spw = '0:1536~2304,1:1536~2304,2:1536~2304,3:1536~2304',
    solint = 'int',
    refant = 'DV04',
    calmode = 'ap')


os.system('rm -rf uid___A002_X4267ef_X1b6.ms.split.bandpass') 
bandpass(vis = 'uid___A002_X4267ef_X1b6.ms.split',
    caltable = 'uid___A002_X4267ef_X1b6.ms.split.bandpass',
    field = '0', # J1924-292
    solint = 'inf',
    combine = 'scan',
    refant = 'DV04',
    solnorm = T,
    bandtype = 'B',
    gaintable = 'uid___A002_X4267ef_X1b6.ms.split.ap_pre_bandpass')

#aU.smoothbandpass('uid___A002_X4267ef_X1b6.ms.split.bandpass') 


# Save flags before gain cal
flagmanager(vis = 'uid___A002_X4267ef_X1b6.ms.split',
    mode = 'save',
    versionname = 'BeforeGainCalibration')


# Gain calibration

os.system('rm -rf uid___A002_X4267ef_X1b6.ms.split.phase_int') 
gaincal(vis = 'uid___A002_X4267ef_X1b6.ms.split',
    caltable = 'uid___A002_X4267ef_X1b6.ms.split.phase_int',
    field = '0~2', # J1924-292,Juno,J1733-130
    solint = 'int',
    refant = 'DV04',
    gaintype = 'G',
    calmode = 'p',
    gaintable = 'uid___A002_X4267ef_X1b6.ms.split.bandpass')


os.system('rm -rf uid___A002_X4267ef_X1b6.ms.split.ampli_inf') 
gaincal(vis = 'uid___A002_X4267ef_X1b6.ms.split',
    caltable = 'uid___A002_X4267ef_X1b6.ms.split.ampli_inf',
    field = '0~2', # J1924-292,Juno,J1733-130
    solint = 'inf',
    refant = 'DV04',
    gaintype = 'T',
    calmode = 'ap',
    gaintable = ['uid___A002_X4267ef_X1b6.ms.split.bandpass', 'uid___A002_X4267ef_X1b6.ms.split.phase_int'])


os.system('rm -rf uid___A002_X4267ef_X1b6.ms.split.flux_inf') 
os.system('rm -rf uid___A002_X4267ef_X1b6.ms.split.fluxscale') 
casalog.setlogfile('uid___A002_X4267ef_X1b6.ms.split.fluxscale')

fluxscale(vis = 'uid___A002_X4267ef_X1b6.ms.split',
    caltable = 'uid___A002_X4267ef_X1b6.ms.split.ampli_inf',
    fluxtable = 'uid___A002_X4267ef_X1b6.ms.split.flux_inf',
    reference = '1') # Juno

casalog.setlogfile('')

os.system('rm -rf uid___A002_X4267ef_X1b6.ms.split.phase_inf') 
gaincal(vis = 'uid___A002_X4267ef_X1b6.ms.split',
    caltable = 'uid___A002_X4267ef_X1b6.ms.split.phase_inf',
    field = '0~2', # J1924-292,Juno,J1733-130
    solint = 'inf',
    refant = 'DV04',
    gaintype = 'G',
    calmode = 'p',
    gaintable = 'uid___A002_X4267ef_X1b6.ms.split.bandpass')


# Save flags before applycal cal
flagmanager(vis = 'uid___A002_X4267ef_X1b6.ms.split',
    mode = 'save',
    versionname = 'BeforeApplycal')


# Application of the bandpass and gain cal tables

for i in ['0', '1']: # J1924-292,Juno
    applycal(vis = 'uid___A002_X4267ef_X1b6.ms.split',
      field = i,
      gaintable = ['uid___A002_X4267ef_X1b6.ms.split.bandpass', 'uid___A002_X4267ef_X1b6.ms.split.phase_int', 'uid___A002_X4267ef_X1b6.ms.split.flux_inf'],
      gainfield = ['', i, i],
      interp = 'linear,linear',
      calwt = F,
      flagbackup = F)


applycal(vis = 'uid___A002_X4267ef_X1b6.ms.split',
    field = '2,3', # HD163296
    gaintable = ['uid___A002_X4267ef_X1b6.ms.split.bandpass', 'uid___A002_X4267ef_X1b6.ms.split.phase_inf', 'uid___A002_X4267ef_X1b6.ms.split.flux_inf'],
    gainfield = ['', '2', '2'], # J1733-130
    interp = 'linear,linear',
    calwt = F,
    flagbackup = F)


#----------------------------------------------------------------------------------
#----- Calibration of uid___A002_X453114_Xb3.ms ----------------------------------
#----------------------------------------------------------------------------------

# Using reference antenna = DV04

asdm = 'uid___A002_X453114_Xb3'

print "# A priori calibration"


# A priori flagging
tflagdata(vis = 'uid___A002_X453114_Xb3.ms',
    mode = 'manual',
    autocorr = T,
    flagbackup = F)

tflagdata(vis = 'uid___A002_X453114_Xb3.ms',
    mode = 'manual',
    intent = '*POINTING*,*SIDEBAND_RATIO*,*ATMOSPHERE*',
    flagbackup = F)

flagcmd(vis = 'uid___A002_X453114_Xb3.ms',
    inpmode = 'table',
    action = 'plot')

flagcmd(vis = 'uid___A002_X453114_Xb3.ms',
    inpmode = 'table',
    action = 'apply')


# Generation of the antenna position cal table

os.system('rm -rf uid___A002_X453114_Xb3.ms.antpos') 
gencal(vis = 'uid___A002_X453114_Xb3.ms',
    caltable = 'uid___A002_X453114_Xb3.ms.antpos',
    caltype = 'antpos',
    antenna = 'DV05,DV08,DV09,DA44,DV13,DA46,DA41,DV03,DA43,DA42,DV17,DV10,DV16,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])
  #  parameter = [0.000407105675771,-6.81411070583e-05,-0.000444496072293,0.000428024806118,0.000192386926053,-7.13865255571e-05,0.000397703996595,0.000303499271338,-0.000252154165917,0.000465348935348,0.000344158009965,0.000228367319351,0.000584870703864,0.000952647723338,-0.000227043112602,0.000247000250965,-0.000508999451995,-0.000311999581754,0.00060776527971,-0.000525096431375,-0.000386487226933,0.000227385910925,-0.000120083500757,-7.74828167994e-07,0.000478540400055,0.000105861353228,3.76092994824e-05,0.000157838716115,0.000211078409778,-0.00037259206194,-7.6349849694e-05,0.000854968838704,1.63702564808e-05,0.000486001936163,7.96910315645e-05,0.000160790103968,6.20906828402e-05,0.000561427652214,0.000133509319652,0.00023328837383,0.000377916222577,-0.000436685084871])


# Application of the WVR, Tsys and antpos cal tables

from recipes.almahelpers import tsysspwmap
tsysmap = tsysspwmap(vis = 'uid___A002_X453114_Xb3.ms', tsystable = 'uid___A002_X453114_Xb3.ms.tsys')

applycal(vis = 'uid___A002_X453114_Xb3.ms',
    field = '0',
    spw = '17,19,21,23',
    gaintable = ['uid___A002_X453114_Xb3.ms.tsys', 'uid___A002_X453114_Xb3.ms.wvr.smooth', 'uid___A002_X453114_Xb3.ms.antpos'],
    gainfield = ['0', '', ''],
    interp = 'linear,linear',
    spwmap = [tsysmap,[],[]],
    calwt = T,
    flagbackup = F)


applycal(vis = 'uid___A002_X453114_Xb3.ms',
    field = '1',
    spw = '17,19,21,23',
    gaintable = ['uid___A002_X453114_Xb3.ms.tsys', 'uid___A002_X453114_Xb3.ms.wvr.smooth', 'uid___A002_X453114_Xb3.ms.antpos'],
    gainfield = ['1', '', ''],
    interp = 'linear,linear',
    spwmap = [tsysmap,[],[]],
    calwt = T,
    flagbackup = F)


applycal(vis = 'uid___A002_X453114_Xb3.ms',
    field = '2',
    spw = '17,19,21,23',
    gaintable = ['uid___A002_X453114_Xb3.ms.tsys', 'uid___A002_X453114_Xb3.ms.wvr.smooth', 'uid___A002_X453114_Xb3.ms.antpos'],
    gainfield = ['2', '', ''],
    interp = 'linear,linear',
    spwmap = [tsysmap,[],[]],
    calwt = T,
    flagbackup = F)


applycal(vis = 'uid___A002_X453114_Xb3.ms',
    field = '3',
    spw = '17,19,21,23',
    gaintable = ['uid___A002_X453114_Xb3.ms.tsys', 'uid___A002_X453114_Xb3.ms.wvr.smooth', 'uid___A002_X453114_Xb3.ms.antpos'],
    gainfield = ['3', '', ''],
    interp = 'linear,linear',
    spwmap = [tsysmap,[],[]],
    calwt = T,
    flagbackup = F)


# Split out science SPWs and time average

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

print "# Calibration"

# Listobs, clear pointing table, and save original flags
os.system('rm -rf uid___A002_X453114_Xb3.ms.split.listobs')
listobs(vis = 'uid___A002_X453114_Xb3.ms.split',
    listfile = 'uid___A002_X453114_Xb3.ms.split.listobs')

tb.open('uid___A002_X453114_Xb3.ms.split/POINTING', nomodify = False)
a = tb.rownumbers()
tb.removerows(a)
tb.close()

flagmanager(vis = 'uid___A002_X453114_Xb3.ms.split',
    mode = 'save',
    versionname = 'Original')


# Initial flagging

# Flagging shadowed data
tflagdata(vis = 'uid___A002_X453114_Xb3.ms.split',
    mode = 'shadow',
    flagbackup = F)

tflagdata(vis = 'uid___A002_X453114_Xb3.ms.split',
    mode = 'manual',
    scan='10',
    flagbackup = F)

# Putting a model for the flux calibrator(s)

setjy(vis = 'uid___A002_X453114_Xb3.ms.split',
    field = '1', # Mars
    spw = '0,1,2,3',
    standard = 'Butler-JPL-Horizons 2010')

os.system('rm -rf uid___A002_X453114_Xb3.ms.split.setjy.field*.png')

if(interact):
  for i in ['1']:
    plotms(vis = 'uid___A002_X453114_Xb3.ms.split',
      xaxis = 'uvdist',
      yaxis = 'amp',
      ydatacolumn = 'model',
      field = i,
      spw = '0,1,2,3',
      avgchannel = '9999',
      coloraxis = 'spw',
      plotfile = 'uid___A002_X453114_Xb3.ms.split.setjy.field'+i+'.png')


# Save flags before bandpass cal
flagmanager(vis = 'uid___A002_X453114_Xb3.ms.split',
    mode = 'save',
    versionname = 'BeforeBandpassCalibration')


# Bandpass calibration

os.system('rm -rf uid___A002_X453114_Xb3.ms.split.ap_pre_bandpass') 

gaincal(vis = 'uid___A002_X453114_Xb3.ms.split',
    caltable = 'uid___A002_X453114_Xb3.ms.split.ap_pre_bandpass',
    field = '0', # J1924-292
    spw = '0:1536~2304,1:1536~2304,2:1536~2304,3:1536~2304',
    solint = 'int',
    refant = 'DV04',
    calmode = 'ap')


os.system('rm -rf uid___A002_X453114_Xb3.ms.split.bandpass') 
bandpass(vis = 'uid___A002_X453114_Xb3.ms.split',
    caltable = 'uid___A002_X453114_Xb3.ms.split.bandpass',
    field = '0', # J1924-292
    solint = 'inf',
    combine = 'scan',
    refant = 'DV04',
    solnorm = T,
    bandtype = 'B',
    gaintable = 'uid___A002_X453114_Xb3.ms.split.ap_pre_bandpass')

os.system('rm -rf uid___A002_X453114_Xb3.ms.split.bandpass_smooth20flat_ri') 

#aU.smoothbandpass('uid___A002_X453114_Xb3.ms.split.bandpass') 


# Save flags before gain cal
flagmanager(vis = 'uid___A002_X453114_Xb3.ms.split',
    mode = 'save',
    versionname = 'BeforeGainCalibration')


# Gain calibration

  # 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('uid___A002_X453114_Xb3.ms.split',field='')

taskname = 'setjy'
default(taskname)
vis = 'uid___A002_X453114_Xb3.ms.split'
field = '2';spw='0';fluxdensity=[1.55,0,0,0];setjy()
field = '2';spw='1';fluxdensity=[1.55,0,0,0];setjy()
field = '2';spw='2';fluxdensity=[1.55,0,0,0];setjy()
field = '2';spw='3';fluxdensity=[1.55,0,0,0];setjy()


os.system('rm -rf uid___A002_X453114_Xb3.ms.split.phase_inf') 
gaincal(vis = 'uid___A002_X453114_Xb3.ms.split',
    caltable = 'uid___A002_X453114_Xb3.ms.split.phase_inf',
    field = '2', # Just J1733-130
    solint = 'inf',
    refant = 'DV04',
    gaintype = 'G',
    calmode = 'ap',
    gaintable = 'uid___A002_X453114_Xb3.ms.split.bandpass')


# Save flags before applycal cal
flagmanager(vis = 'uid___A002_X453114_Xb3.ms.split',
    mode = 'save',
    versionname = 'BeforeApplycal')


# Application of the bandpass and gain cal tables

applycal(vis = 'uid___A002_X453114_Xb3.ms.split',
    field = '0,1,2,3', # HD163296
    #gaintable = ['uid___A002_X453114_Xb3.ms.split.bandpass', 'uid___A002_X453114_Xb3.ms.split.phase_inf', 'uid___A002_X453114_Xb3.ms.split.flux_inf'],
    gaintable = ['uid___A002_X453114_Xb3.ms.split.bandpass', 'uid___A002_X453114_Xb3.ms.split.phase_inf'],
    gainfield = ['', '2'], # J1733-130
    interp = 'linear,linear',
    calwt = F,
    flagbackup = F)


#----------------------------------------------------------------------------------
#----- Concatenate the three datasets and split out science target ----------------
#----------------------------------------------------------------------------------

os.system('rm -rf Concat.ms') 
concat(vis=['uid___A002_X453114_Xb3.ms.split',
            'uid___A002_X4267ef_X1b6.ms.split',
            'uid___A002_X43b3b1_X47.ms.split'],
  dirtol='10arcsec',
  concatvis='Concat.ms')


#----- Flag CM antennas 
tflagdata(vis = 'Concat.ms',
  mode = 'manual',
  antenna = 'CM08,CM05',
  flagbackup = F)


# split out science target

os.system('rm -rf HD163296_Band6.Calib.ms')
split(vis='Concat.ms',
outputvis='HD163296_Band6.Calib.ms',
datacolumn='corrected',
field='HD163296',
#width=128
timebin='',
spw = '')

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