########################################################################
# Data reduction script for HD163296 Band 7:
# - Calibration script: "HD163296_Band7_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_X42aee8_X6a","uid___A002_X4518c5_X3c3","uid___A002_X439aba_X385","uid___A002_X4267ef_X358","uid___A002_X4518c5_X221"]


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

# Using reference antenna = DV04

asdm = 'uid___A002_X4267ef_X358'

print "# A priori calibration"

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

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

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

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


#----- Antenna positions

os.system('rm -rf uid___A002_X4267ef_X358.ms.antpos') 
gencal(vis = 'uid___A002_X4267ef_X358.ms',
    caltable = 'uid___A002_X4267ef_X358.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_X358.ms', tsystable = 'uid___A002_X4267ef_X358.ms.tsys')

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

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

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

applycal(vis = asdm+'.ms',
    field = '3',
    spw = '17,19,21,23',
    gaintable = [asdm+'.ms.tsys', asdm+'.ms.wvr.smooth', asdm+'.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_X358.ms.split') 
split(vis = 'uid___A002_X4267ef_X358.ms',
    outputvis = 'uid___A002_X4267ef_X358.ms.split',
    datacolumn = 'corrected',
    spw = '17,19,21,23',
    timebin = '6.048s',
    keepflags = T)


#----- Calibration Steps
print "# Calibration"

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

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

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


#----- Initial flagging

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


#----- Putting the flux for the phase calibrator to be used as flux calibrator(s). The BP cal., although more intense is more variable in flux. 

setjy(vis = 'uid___A002_X4267ef_X358.ms.split',
    field = '2', # J1733-130
    spw = '0',
    fluxdensity = [1.16555,0,0,0])

setjy(vis = 'uid___A002_X4267ef_X358.ms.split',
    field = '2', # J1733-130
    spw = '1',
    fluxdensity = [1.16729,0,0,0])

setjy(vis = 'uid___A002_X4267ef_X358.ms.split',
    field = '2', # J1733-130
    spw = '2',
    fluxdensity = [1.08586,0,0,0])  

setjy(vis = 'uid___A002_X4267ef_X358.ms.split',
    field = '2', # J1733-130
    spw = '3',
    fluxdensity = [1.05983,0,0,0])  


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


#----- Bandpass calibration
os.system('rm -rf uid___A002_X4267ef_X358.ms.split.ap_pre_bandpass') 
gaincal(vis = 'uid___A002_X4267ef_X358.ms.split',
    caltable = 'uid___A002_X4267ef_X358.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_X358.ms.split.bandpass') 
bandpass(vis = 'uid___A002_X4267ef_X358.ms.split',
    caltable = 'uid___A002_X4267ef_X358.ms.split.bandpass',
    field = '0', # J1924-292
    solint = 'inf',
    combine = 'scan',
    refant = 'DV04',
    solnorm = T,
    bandtype = 'B',
    gaintable = 'uid___A002_X4267ef_X358.ms.split.ap_pre_bandpass')


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

os.system('cp -r uid___A002_X4267ef_X358.ms.split.bandpass uid___A002_X4267ef_X358.ms.split.bandpass_smooth20flat_ri')

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


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

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

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


fluxscale(vis = 'uid___A002_X4267ef_X358.ms.split',
    caltable = 'uid___A002_X4267ef_X358.ms.split.ampli_inf',
    fluxtable = 'uid___A002_X4267ef_X358.ms.split.flux_inf',
    reference = '2') # J1733-130


casalog.setlogfile('')


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


#----- Save flags before applycal cal
flagmanager(vis = 'uid___A002_X4267ef_X358.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_X358.ms.split',
      field = i,
      gaintable = ['uid___A002_X4267ef_X358.ms.split.bandpass_smooth20flat_ri', 'uid___A002_X4267ef_X358.ms.split.phase_int', 'uid___A002_X4267ef_X358.ms.split.flux_inf'],
      gainfield = ['', i, i],
      interp = 'linear,linear',
      calwt = F,
      flagbackup = F)

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


#No good solutions solutions for swp 0. Skipt it

#----- Inspecting calibrated data
if(interact):
	for asdm in basename:
		plotms(vis = asdm+'.ms.split', xaxis='uvdist', yaxis='amp',
		        ydatacolumn='corrected', field='3c279',
		        averagedata=True, avgchannel='3840', avgtime='',
		        avgscan=F, avgbaseline=F, coloraxis='corr')
		dummy_string = raw_input("Next plot for "+asdm+" . Hit <Enter> to continue.")
		plotms(vis = asdm+'.ms.split', xaxis='uvdist', yaxis='phase',
		        ydatacolumn='corrected', field='3c279',
			avgchannel='3840', avgscan=F, avgbaseline=F, coloraxis='corr')
		dummy_string = raw_input("Next plot for "+asdm+" . Hit <Enter> to continue.")
		plotms(vis = asdm+'.ms.split', xaxis='time', yaxis='amp',
		        ydatacolumn='corrected', field='',
			averagedata=True, avgchannel='3840', avgtime='',
			avgscan=F, avgbaseline=F, coloraxis='corr')
		dummy_string = raw_input("Next plot for "+asdm+" . Hit <Enter> to continue.")
		plotms(vis = asdm+'.ms.split', xaxis='time', yaxis='phase',
		        ydatacolumn='corrected', field='',
			avgchannel='3840', avgscan=F, avgbaseline=F, coloraxis='corr')


#----- Split out calibrated data for X358
split(vis = 'uid___A002_X4267ef_X358.ms.split',
    outputvis = 'uid___A002_X4267ef_X358.corrected_spw123.ms',
    spw='1,2,3',
    datacolumn = 'corrected',
    keepflags = T)

split(vis = 'uid___A002_X4267ef_X358.ms.split',
    outputvis = 'uid___A002_X4267ef_X358.HD163296_spw123.calibrated',
    spw='1,2,3',
    field = '3',
    datacolumn = 'corrected',
    keepflags = T)


#----------------------------------------------------------------------------------
#----- Calibration of uid___A002_X42aee8_X6a.ms -----------------------------------
#----------------------------------------------------------------------------------

# Using reference antenna = DV04

print "# A priori calibration"

#----- Running fixplanets on fields with 0,0 coordinates
fixplanets(vis = 'uid___A002_X42aee8_X6a.ms',
    field = '1', # Neptune
    fixuvw = T)


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

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

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

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


#----- Generation of the antenna position cal table

os.system('rm -rf uid___A002_X42aee8_X6a.ms.antpos') 
gencal(vis = 'uid___A002_X42aee8_X6a.ms',
    caltable = 'uid___A002_X42aee8_X6a.ms.antpos',
    caltype = 'antpos',
    antenna = 'DV12,DV08,DV09,DV11,DA44,DV13,DV05,DA41,DV14,DA43,DA42,DV02,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,0,0,0,0,0,0])
    parameter = [-0.000133162034778,0.00023574780856,0.000212197651384,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,3.7086635797e-05,-0.000230659921057,-0.000312076244311,0.000407105675771,-6.81411070583e-05,-0.000444496072293,0.00060776527971,-0.000525096431375,-0.000386487226933,0.000293836351762,0.000269524465858,-8.05776686625e-05,0.000478540400055,0.000105861353228,3.76092994824e-05,0.000157838716115,0.000211078409778,-0.00037259206194,-1.03950910152e-05,-0.000460584326197,-0.000282536759635,-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_X42aee8_X6a.ms', tsystable = 'uid___A002_X42aee8_X6a.ms.tsys')

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

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

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

applycal(vis = 'uid___A002_X42aee8_X6a.ms',
    field = '3',
    spw = '17,19,21,23',
    gaintable = ['uid___A002_X42aee8_X6a.ms.tsys', 'uid___A002_X42aee8_X6a.ms.wvr.smooth', 'uid___A002_X42aee8_X6a.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_X42aee8_X6a.ms.split') 
split(vis = 'uid___A002_X42aee8_X6a.ms',
    outputvis = 'uid___A002_X42aee8_X6a.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_X42aee8_X6a.ms.split.listobs')
listobs(vis = 'uid___A002_X42aee8_X6a.ms.split',
    listfile = 'uid___A002_X42aee8_X6a.ms.split.listobs')


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


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


#----- Initial flagging

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


# Flagging atmospheric line(s)
tflagdata(vis = 'uid___A002_X42aee8_X6a.ms.split',
    mode = 'manual',
    spw = '2:0~3839',
    field = '1',
    flagbackup = F)

# Flagging atmospheric line(s)
tflagdata(vis = 'uid___A002_X42aee8_X6a.ms.split',
    mode = 'manual',
    spw = '3:0~3839',
    field = '1',
    flagbackup = F)

# Flagging atmospheric line(s antenna with very high Tsys.
tflagdata(vis = 'uid___A002_X42aee8_X6a.ms.split',
    mode = 'manual',
    spw = '3',
    field = '',
    antenna = 'DA43',
    flagbackup = F)


#----- Putting a model for the flux calibrator(s)

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


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

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


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


#----- Bandpass calibration

os.system('rm -rf uid___A002_X42aee8_X6a.ms.split.ap_pre_bandpass') 
gaincal(vis = 'uid___A002_X42aee8_X6a.ms.split',
    caltable = 'uid___A002_X42aee8_X6a.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_X42aee8_X6a.ms.split.bandpass') 
bandpass(vis = 'uid___A002_X42aee8_X6a.ms.split',
    caltable = 'uid___A002_X42aee8_X6a.ms.split.bandpass',
    field = '0', # J1924-292
    solint = 'inf',
    combine = 'scan',
    refant = 'DV04',
    solnorm = T,
    bandtype = 'B',
    gaintable = 'uid___A002_X42aee8_X6a.ms.split.ap_pre_bandpass')
os.system('rm -rf uid___A002_X42aee8_X6a.ms.split.bandpass_smooth20flat_ri') 

os.system('cp -r  uid___A002_X42aee8_X6a.ms.split.bandpass uid___A002_X42aee8_X6a.ms.split.bandpass_smooth20flat_ri')


#----- Save flags before gain cal
flagmanager(vis = 'uid___A002_X42aee8_X6a.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_X42aee8_X6a.ms.split',field='J1733-130')

os.system('rm -rf uid___A002_X42aee8_X6a.ms.split.phase_short_int') 
gaincal(vis = 'uid___A002_X42aee8_X6a.ms.split',
    caltable = 'uid___A002_X42aee8_X6a.ms.split.phase_short_int',
    field = '0~2', # J1924-292,Neptune,J1733-130
    selectdata = T,
    antenna = 'DA41,DA42,DV04,DV09&',
    solint = 'int',
    refant = 'DV04',
    minblperant = 3,
    minsnr = 2.0,
    gaintype = 'G',
    calmode = 'p',
    gaintable = 'uid___A002_X42aee8_X6a.ms.split.bandpass_smooth20flat_ri')

os.system('rm -rf uid___A002_X42aee8_X6a.ms.split.ampli_short_inf') 
gaincal(vis = 'uid___A002_X42aee8_X6a.ms.split',
    caltable = 'uid___A002_X42aee8_X6a.ms.split.ampli_short_inf',
    field = '0~2', # J1924-292,Neptune,J1733-130
    selectdata = T,
    antenna = 'DA41,DA42,DV04,DV09&',
    solint = 'inf',
    refant = 'DV04',
    minblperant = 3,
    minsnr = 2.0,
    gaintype = 'T',
    calmode = 'ap',
    gaintable = ['uid___A002_X42aee8_X6a.ms.split.bandpass_smooth20flat_ri', 'uid___A002_X42aee8_X6a.ms.split.phase_short_int'])

os.system('rm -rf uid___A002_X42aee8_X6a.ms.split.flux_short_inf') 
os.system('rm -rf uid___A002_X42aee8_X6a.ms.split.fluxscale') 
casalog.setlogfile('uid___A002_X42aee8_X6a.ms.split.fluxscale')


fluxscale(vis = 'uid___A002_X42aee8_X6a.ms.split',
    caltable = 'uid___A002_X42aee8_X6a.ms.split.ampli_short_inf',
    fluxtable = 'uid___A002_X42aee8_X6a.ms.split.flux_short_inf', 
    reference = '1',
    refspwmap=[0,1,1,1]) # Neptune

casalog.setlogfile('')

f = open('uid___A002_X42aee8_X6a.ms.split.fluxscale')
fc = f.readlines()
f.close()

for phaseCalName in ['J1733-130']:
    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 = 'uid___A002_X42aee8_X6a.ms.split',
          field = phaseCalName.replace(';','*;').split(';')[0],
          spw = spwId,
          fluxdensity = [flux,0,0,0])


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

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

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


#----- Save flags before applycal cal
flagmanager(vis = 'uid___A002_X42aee8_X6a.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_X42aee8_X6a.ms.split',
      field = i,
      gaintable = ['uid___A002_X42aee8_X6a.ms.split.bandpass_smooth20flat_ri', 'uid___A002_X42aee8_X6a.ms.split.phase_int', 'uid___A002_X42aee8_X6a.ms.split.flux_inf'],
      gainfield = ['', i, i],
      interp = 'linear,linear',
      calwt = F,
      flagbackup = F)

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


tflagdata(vis = 'uid___A002_X42aee8_X6a.ms.split',
    mode = 'manual',
    spw = '',
    field = '',
    antenna = 'CM04,CM08',
    flagbackup = F)

split(vis = 'uid___A002_X42aee8_X6a.ms.split',
    outputvis = 'uid___A002_X42aee8_X6a.corrected.ms',
    datacolumn = 'corrected',
    keepflags = T)

split(vis = 'uid___A002_X42aee8_X6a.ms.split',
    outputvis = 'uid___A002_X42aee8_X6a.HD163296.calibrated',
    datacolumn = 'corrected',
    field='3',   
    keepflags = T)


#----------------------------------------------------------------------------------
#----- Calibration of uid___A002_X439aba_X385.ms ----------------------------------
#----------------------------------------------------------------------------------

# Using reference antenna = DV04

asdm='uid___A002_X439aba_X385'

print "# A priori calibration"

#----- Running fixplanets on fields with 0,0 coordinates
fixplanets(vis = 'uid___A002_X439aba_X385.ms',
    field = '1', # Neptune
    fixuvw = T)

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

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

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

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


#----- Generation of the antenna position cal table

os.system('rm -rf uid___A002_X439aba_X385.ms.antpos') 
gencal(vis = 'uid___A002_X439aba_X385.ms',
    caltable = 'uid___A002_X439aba_X385.ms.antpos',
    caltype = 'antpos',
    antenna = 'DV05,DV08,DV09,DV11,DA44,DV13,DA46,DA41,DV03,DA43,DA42,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])
    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,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_X439aba_X385.ms', tsystable = 'uid___A002_X439aba_X385.ms.tsys')

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

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

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

applycal(vis = 'uid___A002_X439aba_X385.ms',
    field = '3',
    spw = '17,19,21,23',
    gaintable = ['uid___A002_X439aba_X385.ms.tsys', 'uid___A002_X439aba_X385.ms.wvr.smooth', 'uid___A002_X439aba_X385.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_X439aba_X385.ms.split') 
split(vis = 'uid___A002_X439aba_X385.ms',
    outputvis = 'uid___A002_X439aba_X385.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_X439aba_X385.ms.split.listobs')
listobs(vis = 'uid___A002_X439aba_X385.ms.split',
    listfile = 'uid___A002_X439aba_X385.ms.split.listobs')

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

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


#----- Initial flagging

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

# Flagging atmospheric line(s)
tflagdata(vis = 'uid___A002_X439aba_X385.ms.split',
    mode = 'manual',
    spw = '2:0~3839',
    field = '1',
    flagbackup = F)

# Flagging atmospheric line(s)
tflagdata(vis = 'uid___A002_X439aba_X385.ms.split',
    mode = 'manual',
    spw = '3:0~3839',
    field = '1',
    flagbackup = F)


#----- Putting a model for the flux calibrator(s)
setjy(vis = 'uid___A002_X439aba_X385.ms.split',
    field = '1', # Neptune
    spw = '0,1,2,3',
    standard = 'Butler-JPL-Horizons 2010')

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

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


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


#----- Bandpass calibration

os.system('rm -rf uid___A002_X439aba_X385.ms.split.ap_pre_bandpass') 
gaincal(vis = 'uid___A002_X439aba_X385.ms.split',
    caltable = 'uid___A002_X439aba_X385.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_X439aba_X385.ms.split.bandpass') 
bandpass(vis = 'uid___A002_X439aba_X385.ms.split',
    caltable = 'uid___A002_X439aba_X385.ms.split.bandpass',
    field = '0', # J1924-292
    solint = 'inf',
    combine = 'scan',
    refant = 'DV04',
    solnorm = T,
    bandtype = 'B',
    gaintable = 'uid___A002_X439aba_X385.ms.split.ap_pre_bandpass')
os.system('rm -rf uid___A002_X439aba_X385.ms.split.bandpass_smooth20flat_ri') 

os.system('cp -r uid___A002_X439aba_X385.ms.split.bandpass uid___A002_X439aba_X385.ms.split.bandpass_smooth20flat_ri')


#----- Save flags before gain cal
flagmanager(vis = 'uid___A002_X439aba_X385.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_X439aba_X385.ms.split',field='J1733-130')

os.system('rm -rf uid___A002_X439aba_X385.ms.split.phase_short_int') 
gaincal(vis = 'uid___A002_X439aba_X385.ms.split',
    caltable = 'uid___A002_X439aba_X385.ms.split.phase_short_int',
    field = '0~2', # J1924-292,Neptune,J1733-130
    selectdata = T,
    antenna = 'DA41,DA42,DV04,DV09&',
    solint = 'int',
    refant = 'DV04',
    minblperant = 3,
    minsnr = 2.0,
    gaintype = 'G',
    calmode = 'p',
    gaintable = 'uid___A002_X439aba_X385.ms.split.bandpass_smooth20flat_ri')

os.system('rm -rf uid___A002_X439aba_X385.ms.split.ampli_short_inf') 
gaincal(vis = 'uid___A002_X439aba_X385.ms.split',
    caltable = 'uid___A002_X439aba_X385.ms.split.ampli_short_inf',
    field = '0~2', # J1924-292,Neptune,J1733-130
    selectdata = T,
    antenna = 'DA41,DA42,DV04,DV09&',
    solint = 'inf',
    refant = 'DV04',
    minblperant = 3,
    minsnr = 2.0,
    gaintype = 'T',
    calmode = 'ap',
    gaintable = ['uid___A002_X439aba_X385.ms.split.bandpass_smooth20flat_ri', 'uid___A002_X439aba_X385.ms.split.phase_short_int'])

os.system('rm -rf uid___A002_X439aba_X385.ms.split.flux_short_inf') 
os.system('rm -rf uid___A002_X439aba_X385.ms.split.fluxscale') 
casalog.setlogfile('uid___A002_X439aba_X385.ms.split.fluxscale')


fluxscale(vis = 'uid___A002_X439aba_X385.ms.split',
    caltable = 'uid___A002_X439aba_X385.ms.split.ampli_short_inf',
    fluxtable = 'uid___A002_X439aba_X385.ms.split.flux_short_inf',
    reference = '1',
    refspwmap=[0,1,1,1]) # Neptune

casalog.setlogfile('')

f = open('uid___A002_X439aba_X385.ms.split.fluxscale')
fc = f.readlines()
f.close()


for phaseCalName in ['J1733-130']:
    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 = 'uid___A002_X439aba_X385.ms.split',
          field = phaseCalName.replace(';','*;').split(';')[0],
          spw = spwId,
          fluxdensity = [flux,0,0,0])


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

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


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


#----- Save flags before applycal cal
flagmanager(vis = 'uid___A002_X439aba_X385.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_X439aba_X385.ms.split',
      field = i,
      gaintable = ['uid___A002_X439aba_X385.ms.split.bandpass_smooth20flat_ri', 'uid___A002_X439aba_X385.ms.split.phase_int', 'uid___A002_X439aba_X385.ms.split.flux_inf'],
      gainfield = ['', i, i],
      interp = 'linear,linear',
      calwt = F,
      flagbackup = F)


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

split(vis = 'uid___A002_X439aba_X385.ms.split',
    outputvis = 'uid___A002_X439aba_X385.corrected.ms',
    datacolumn = 'corrected',
    keepflags = T)

split(vis = 'uid___A002_X439aba_X385.ms.split',
    outputvis = 'uid___A002_X439aba_X385.HD163296.calibrated',
    field = '3',
    datacolumn = 'corrected',
    keepflags = T)


#----------------------------------------------------------------------------------
#----- Calibration of uid___A002_X4518c5_X3c3.ms ----------------------------------
#----------------------------------------------------------------------------------

# Using reference antenna = DV04

print "# A priori calibration"

#----- Running fixplanets on fields with 0,0 coordinates
fixplanets(vis = 'uid___A002_X4518c5_X3c3.ms',
    field = '1', # Neptune
    fixuvw = T)


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

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

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

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


#----- Generation of the antenna position cal table

os.system('rm -rf uid___A002_X4518c5_X3c3.ms.antpos') 
gencal(vis = 'uid___A002_X4518c5_X3c3.ms',
    caltable = 'uid___A002_X4518c5_X3c3.ms.antpos',
    caltype = 'antpos',
    antenna = 'DV05,DV08,DV09,DV11,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,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.00023328837383,0.000377916222577,-0.000436685084871])


#----- Application of the WVR, Tsys and antpos cal tables

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

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

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

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

applycal(vis = 'uid___A002_X4518c5_X3c3.ms',
    field = '3',
    spw = '17,19,21,23',
    gaintable = ['uid___A002_X4518c5_X3c3.ms.tsys', 'uid___A002_X4518c5_X3c3.ms.wvr.smooth', 'uid___A002_X4518c5_X3c3.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_X4518c5_X3c3.ms.split') 
split(vis = 'uid___A002_X4518c5_X3c3.ms',
    outputvis = 'uid___A002_X4518c5_X3c3.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_X4518c5_X3c3.ms.split.listobs')
listobs(vis = 'uid___A002_X4518c5_X3c3.ms.split',
    listfile = 'uid___A002_X4518c5_X3c3.ms.split.listobs')

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

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


#----- Initial flagging

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

# Flagging atmospheric line(s)
tflagdata(vis = 'uid___A002_X4518c5_X3c3.ms.split',
    mode = 'manual',
    spw = '2:0~3839',
    field = '1',
    flagbackup = F)

# Flagging atmospheric line(s)
tflagdata(vis = 'uid___A002_X4518c5_X3c3.ms.split',
    mode = 'manual',
    spw = '3:0~3839',
    field = '1',
    flagbackup = F)


#----- Putting a model for the flux calibrator(s)
setjy(vis = 'uid___A002_X4518c5_X3c3.ms.split',
    field = '1', # Neptune
    spw = '0,1,2,3',
    standard = 'Butler-JPL-Horizons 2010')

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

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


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


#----- Bandpass calibration

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

gaincal(vis = 'uid___A002_X4518c5_X3c3.ms.split',
    caltable = 'uid___A002_X4518c5_X3c3.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_X4518c5_X3c3.ms.split.bandpass') 
bandpass(vis = 'uid___A002_X4518c5_X3c3.ms.split',
    caltable = 'uid___A002_X4518c5_X3c3.ms.split.bandpass',
    field = '0', # J1924-292
    solint = 'inf',
    combine = 'scan',
    refant = 'DV04',
    solnorm = T,
    bandtype = 'B',
    gaintable = 'uid___A002_X4518c5_X3c3.ms.split.ap_pre_bandpass')
os.system('rm -rf uid___A002_X4518c5_X3c3.ms.split.bandpass_smooth20flat_ri') 

os.system('cp -r uid___A002_X4518c5_X3c3.ms.split.bandpass uid___A002_X4518c5_X3c3.ms.split.bandpass_smooth20flat_ri')


#----- Save flags before gain cal
flagmanager(vis = 'uid___A002_X4518c5_X3c3.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_X4518c5_X3c3.ms.split',field='J1733-130')

os.system('rm -rf uid___A002_X4518c5_X3c3.ms.split.phase_short_int') 
gaincal(vis = 'uid___A002_X4518c5_X3c3.ms.split',
    caltable = 'uid___A002_X4518c5_X3c3.ms.split.phase_short_int',
    field = '0~2', # J1924-292,Neptune,J1733-130
    selectdata = T,
    antenna = 'DA41,DV04,DV08,DV19&',
    solint = 'int',
    refant = 'DV04',
    minblperant = 3,
    minsnr = 2.0,
    gaintype = 'G',
    calmode = 'p',
    gaintable = 'uid___A002_X4518c5_X3c3.ms.split.bandpass_smooth20flat_ri')

os.system('rm -rf uid___A002_X4518c5_X3c3.ms.split.ampli_short_inf') 
gaincal(vis = 'uid___A002_X4518c5_X3c3.ms.split',
    caltable = 'uid___A002_X4518c5_X3c3.ms.split.ampli_short_inf',
    field = '0~2', # J1924-292,Neptune,J1733-130
    selectdata = T,
    antenna = 'DA41,DV04,DV08,DV19&',
    solint = 'inf',
    refant = 'DV04',
    minblperant = 3,
    minsnr = 2.0,
    gaintype = 'T',
    calmode = 'ap',
    gaintable = ['uid___A002_X4518c5_X3c3.ms.split.bandpass_smooth20flat_ri', 'uid___A002_X4518c5_X3c3.ms.split.phase_short_int'])

os.system('rm -rf uid___A002_X4518c5_X3c3.ms.split.flux_short_inf') 
os.system('rm -rf uid___A002_X4518c5_X3c3.ms.split.fluxscale') 
casalog.setlogfile('uid___A002_X4518c5_X3c3.ms.split.fluxscale')

fluxscale(vis = 'uid___A002_X4518c5_X3c3.ms.split',
    caltable = 'uid___A002_X4518c5_X3c3.ms.split.ampli_short_inf',
    fluxtable = 'uid___A002_X4518c5_X3c3.ms.split.flux_short_inf',
    reference = '1',
    refspwmap=[0,1,1,1]) # Neptune

casalog.setlogfile('')

f = open('uid___A002_X4518c5_X3c3.ms.split.fluxscale')
fc = f.readlines()
f.close()

for phaseCalName in ['J1733-130']:
    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 = 'uid___A002_X4518c5_X3c3.ms.split',
          field = phaseCalName.replace(';','*;').split(';')[0],
          spw = spwId,
          fluxdensity = [flux,0,0,0])

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

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

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


#----- Save flags before applycal cal
flagmanager(vis = 'uid___A002_X4518c5_X3c3.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_X4518c5_X3c3.ms.split',
      field = i,
      gaintable = ['uid___A002_X4518c5_X3c3.ms.split.bandpass_smooth20flat_ri', 'uid___A002_X4518c5_X3c3.ms.split.phase_int', 'uid___A002_X4518c5_X3c3.ms.split.flux_inf'],
      gainfield = ['', i, i],
      interp = 'linear,linear',
      calwt = F,
      flagbackup = F)

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

split(vis = 'uid___A002_X4518c5_X3c3.ms.split',
    outputvis = 'uid___A002_X4518c5_X3c3.corrected.ms',
    datacolumn = 'corrected',
    keepflags = T)

split(vis = 'uid___A002_X4518c5_X3c3.ms.split',
    outputvis = 'uid___A002_X4518c5_X3c3.HD163296.calibrated',
    field = '3',
    datacolumn = 'corrected',
    keepflags = T)


#----------------------------------------------------------------------------------
#----- Calibration of uid___A002_X4518c5_X221.ms ----------------------------------
#----------------------------------------------------------------------------------

# Using reference antenna = DV04

print "# A priori calibration"

#----- Running fixplanets on fields with 0,0 coordinates
fixplanets(vis = 'uid___A002_X4518c5_X221.ms',
    field = '1', # Neptune
    fixuvw = T)


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

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

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

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


#----- Generation of the antenna position cal table

os.system('rm -rf uid___A002_X4518c5_X221.ms.antpos') 
gencal(vis = 'uid___A002_X4518c5_X221.ms',
    caltable = 'uid___A002_X4518c5_X221.ms.antpos',
    caltype = 'antpos',
    antenna = 'DV05,DV08,DV09,DV11,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,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.00023328837383,0.000377916222577,-0.000436685084871])


#----- Application of the WVR, Tsys and antpos cal tables

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

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

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

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

applycal(vis = 'uid___A002_X4518c5_X221.ms',
    field = '3',
    spw = '17,19,21,23',
    gaintable = ['uid___A002_X4518c5_X221.ms.tsys', 'uid___A002_X4518c5_X221.ms.wvr.smooth', 'uid___A002_X4518c5_X221.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_X4518c5_X221.ms.split') 
split(vis = 'uid___A002_X4518c5_X221.ms',
    outputvis = 'uid___A002_X4518c5_X221.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_X4518c5_X221.ms.split.listobs')
listobs(vis = 'uid___A002_X4518c5_X221.ms.split',
    listfile = 'uid___A002_X4518c5_X221.ms.split.listobs')

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

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


#----- Initial flagging

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

# Flagging atmospheric line(s)
tflagdata(vis = 'uid___A002_X4518c5_X221.ms.split',
    mode = 'manual',
    spw = '2:0~3839',
    field = '1',
    flagbackup = F)

# Flagging atmospheric line(s)
tflagdata(vis = 'uid___A002_X4518c5_X221.ms.split',
    mode = 'manual',
    spw = '3:0~3839',
    field = '1',
    flagbackup = F)


#----- Putting a model for the flux calibrator(s)
setjy(vis = 'uid___A002_X4518c5_X221.ms.split',
    field = '1', # Neptune
    spw = '0,1,2,3',
    standard = 'Butler-JPL-Horizons 2010')

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

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


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


#----- Bandpass calibration
os.system('rm -rf uid___A002_X4518c5_X221.ms.split.ap_pre_bandpass') 

gaincal(vis = 'uid___A002_X4518c5_X221.ms.split',
    caltable = 'uid___A002_X4518c5_X221.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_X4518c5_X221.ms.split.bandpass') 
bandpass(vis = 'uid___A002_X4518c5_X221.ms.split',
    caltable = 'uid___A002_X4518c5_X221.ms.split.bandpass',
    field = '0', # J1924-292
    solint = 'inf',
    combine = 'scan',
    refant = 'DV04',
    solnorm = T,
    bandtype = 'B',
    gaintable = 'uid___A002_X4518c5_X221.ms.split.ap_pre_bandpass')

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

os.system('cp -r uid___A002_X4518c5_X221.ms.split.bandpass  uid___A002_X4518c5_X221.ms.split.bandpass_smooth20flat_ri')


#----- Save flags before gain cal
flagmanager(vis = 'uid___A002_X4518c5_X221.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_X4518c5_X221.ms.split',field='J1733-130')

os.system('rm -rf uid___A002_X4518c5_X221.ms.split.phase_short_int') 
gaincal(vis = 'uid___A002_X4518c5_X221.ms.split',
    caltable = 'uid___A002_X4518c5_X221.ms.split.phase_short_int',
    field = '0~2', # J1924-292,Neptune,J1733-130
    selectdata = T,
    antenna = 'DA41,DV04,DV08,DV19&',
    solint = 'int',
    refant = 'DV04',
    minblperant = 3,
    minsnr = 2.0,
    gaintype = 'G',
    calmode = 'p',
    gaintable = 'uid___A002_X4518c5_X221.ms.split.bandpass_smooth20flat_ri')

os.system('rm -rf uid___A002_X4518c5_X221.ms.split.ampli_short_inf') 
gaincal(vis = 'uid___A002_X4518c5_X221.ms.split',
    caltable = 'uid___A002_X4518c5_X221.ms.split.ampli_short_inf',
    field = '0~2', # J1924-292,Neptune,J1733-130
    selectdata = T,
    antenna = 'DA41,DV04,DV08,DV19&',
    solint = 'inf',
    refant = 'DV04',
    minblperant = 3,
    minsnr = 2.0,
    gaintype = 'T',
    calmode = 'ap',
    gaintable = ['uid___A002_X4518c5_X221.ms.split.bandpass_smooth20flat_ri', 'uid___A002_X4518c5_X221.ms.split.phase_short_int'])

os.system('rm -rf uid___A002_X4518c5_X221.ms.split.flux_short_inf') 
os.system('rm -rf uid___A002_X4518c5_X221.ms.split.fluxscale') 
casalog.setlogfile('uid___A002_X4518c5_X221.ms.split.fluxscale')

fluxscale(vis = 'uid___A002_X4518c5_X221.ms.split',
    caltable = 'uid___A002_X4518c5_X221.ms.split.ampli_short_inf',
    fluxtable = 'uid___A002_X4518c5_X221.ms.split.flux_short_inf',
    reference = '1',
    refspwmap=[0,1,1,1]) # Neptune

casalog.setlogfile('')

f = open('uid___A002_X4518c5_X221.ms.split.fluxscale')
fc = f.readlines()
f.close()

for phaseCalName in ['J1733-130']:
    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 = 'uid___A002_X4518c5_X221.ms.split',
          field = phaseCalName.replace(';','*;').split(';')[0],
          spw = spwId,
          fluxdensity = [flux,0,0,0])

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

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

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


#----- Save flags before applycal cal
flagmanager(vis = 'uid___A002_X4518c5_X221.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_X4518c5_X221.ms.split',
      field = i,
      gaintable = ['uid___A002_X4518c5_X221.ms.split.bandpass_smooth20flat_ri', 'uid___A002_X4518c5_X221.ms.split.phase_int', 'uid___A002_X4518c5_X221.ms.split.flux_inf'],
      gainfield = ['', i, i],
      interp = 'linear,linear',
      calwt = F,
      flagbackup = F)

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


split(vis = 'uid___A002_X4518c5_X221.ms.split',
    outputvis = 'uid___A002_X4518c5_X221.corrected.ms',
    datacolumn = 'corrected',
    keepflags = T)

#----- Split out science target 
split(vis = 'uid___A002_X4518c5_X221.ms.split',
    outputvis = 'uid___A002_X4518c5_X221.HD163296.calibrated',
    field = '3',
    datacolumn = 'corrected',
    keepflags = T)


#----- Concat all with a position tolerance of 5" 
concat (vis=['uid___A002_X42aee8_X6a.HD163296.calibrated','uid___A002_X4518c5_X3c3.HD163296.calibrated','uid___A002_X439aba_X385.HD163296.calibrated','uid___A002_X4267ef_X358.HD163296_spw123.calibrated','uid___A002_X4518c5_X221.HD163296.calibrated'],dirtol = '5.arcsec', concatvis='HD163296_Band7_concat.ms')

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


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