########################################################################
# Calibration Script for SgrA* Band 6 ALMA data                        #
# See also imaging script                                              #
########################################################################
import re
from glob import glob
"""

You need:

SgrA_Band6_UnCalibratedMSAndTablesForReduction.tgz

which contains;
Measurement sets
* uid___A002_X22afc4_X5b9.ms
* uid___A002_X22afc4_X716.ms
* uid___A002_X22afc4_X874.ms

wvr (water vapour radiometry) and tsys (system temperature) tables
(untar before use) contains
X5b9.tsys.fdm    X716.tsys.fdm    X874.tsys.fdm
X5b9.wvr  X716.wvr  (no table for 3rd dataset)
These tables are generated from instrumental measurements and have
been checked. 

Comments on the data:

From listobs:
12 antennas
spw 1,3,5,7 1.875 GHz in 3840 chans (FDM, used for most operations here)
spw 9,11,13,15  2GHz in 128 chans   (TDM, used to derive Tsys tables) 

spw0    1           3        5         7
WVR   217137.156 218939.506 230960.494 232762.844

X5b9.ms X716.ms
0           1       2      3       4~10
J1924-292, Titan, 3c279, nrao530, SgrA x7 
pointing    flux  BP     phref    target

X874.ms
0           1       2      3~9
J1924-292, 3c279, nrao530, SgrA x7 

DV05 was known to be bad and is to be flagged.

In the second and third data sets, 3C279 and Titan were observed at
elevations below 30 deg and were too noisy and/or had too deep
atmospheric lines to be useful and are not used. If they are retained,
note that none of these data have useful Tsys values for DV09 and in
the third data set 3c279 has no Tsys measurements.

"""

# Input MS
svis=['uid___A002_X22afc4_X5b9.ms','uid___A002_X22afc4_X716.ms','uid___A002_X22afc4_X874.ms']

roots=['X5b9','X716','X874']

#------------------------#
#------------------------#
#------ Script Sections ------#
#-- use value T to activate --#

#- TEST AND OPTIONAL STEPS -#

test= F			# Quick test to see if script looking in right place  
do_plots =F             # Triggers multiple instances of plotcal throughout script
                        # Add more instances as required
for_plotms =F           #applycal purely for checking in plotms (not included)


#- INITIAL INSPECTION/FLAGGING -#
do_listobs =F           # List observation data. 

initial_state =F        # Back up curent state of ms
initial_flag =F	        # Initial flags (shadowing, known bad data etc.)
flag_PtgAtm =F          # Flag pointing scans and WVR spectral windows
      
#- INSTRUMENTAL CALIBRATION -#

# Apply calibration to FDM and split out
FDM_apply_wvr_tsys =T	# Apply Tsys and wvr tables for FDM
split_FDM  =T		# Split out FDM SPWs  

correct_phaserefpos =T # Update phase reference source position

# Specific flagging
flag_all  =T            # Flag some atmospheric lines
flag_titan =T		# Flag lines seen worse in Titan 
flag_badpols  =T        # Flag data bad in one or both polarizations

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

#-TITAN-#
fix_titan =T		# Fix titan position ENSURE NO OTHER PROCESSES INTERRUPT 
set_jy_titan =T	        # Set model of Titan's flux distribution 

#- BASIC CALIBRATION USING CAL SOURCES -#

phase_bp =T	        # Generate phase correction table for 3c279
gen_BP  =T		# Generate BP calibration table from 3c279

first_phase_cals =T	# Generate phase cal table for calibrator sources 
first_phase_target =T	# Generate phase cal table to apply to target sources 
bootstrap =T		# Bootstrap Flux density scale for calibrators
set_phaseref_flux =T    # Set derived phase ref flux in all datasets
amp_cals =T             # Generate amplitude cal table for phase ref/target

# These can be omitted if secondary bandpass cal and blcal are made
apply_cal_cals =F	# Apply all calibrations to phase ref
apply_cal_target =F	# Apply all calibrations to target sources

#- ADDITIONAL CALIBRRATION-#
phase_BP =T		# Create secondary bandpass correction from
			# the phase cal nrao530 
phref_blcal =T          #  Baseline-based calibration 

#- APPLY CALIBRATION, EXTRACT AND IMAGE PHASE REF -#

apply_cal_cals_2 =T	# Apply all calibration to calibrator sources
phsref_split=T          # Split out phase ref and concatenate
phsref_image=T          # Image phase reference source (non-interactive)

#- APPLY CALIBRATION, EXTRACT TARGET -#			
apply_cal_target_2 =T	# Apply all calibration to target data 

# These averaging intervals reduce data volume/processing times
# but do not degrade field of view and seem OK for self-calibration
averagetime = '30s'     # Time averaging for target
averagechan = 5         # Channels averaged for target
target_split=T          # Split out and concatenate target data


#-IMAGING-#  see separate script

#------------------------#

#------------------------#
#------------------------#
#--------test------------#

if (test):
	os.system('ls -lh uid*')

#------------------------#
#------------------------#
#-----List Obs Data------#

if (do_listobs):
	default(listobs)		#show MS information
	verbose= F
	for i in range(len(svis)):
    		vis=svis[i]
    		listobs()
	print "<< Check logger for details of what is in your data\n"
	
#------------------------#
#------------------------#
#--- initial flagging ---#

if (initial_state):
	print "<< Save current state of flagging as Original\n"
	default(flagmanager)	       
	mode='save'
	versionname='Original'
	for i in range(len(svis)):
    		vis=svis[i]
    		flagmanager()

if (initial_flag):
	print "<< Flag autocorrelations\n"	
	default(flagautocorr)		
	for i in range(len(svis)):
    		vis=svis[i]
    		flagautocorr()

	print "<< Flag DV05, known to have no good data\n" 	
	default(flagdata2)		
	selectdata = T
	antenna='3&*'
	manualflag= T
	mf_antenna='3&*'
	for i in range(len(svis)):
    		vis=svis[i]
    		flagdata2()

	print "<< Flag shadowed antennas\n"
	default(flagdata2)		
	shadow= T
	diameter=12.
	for i in range(len(svis)):
    		vis=svis[i]
    		flagdata2()

	print "<< Save flagged state of ms\n"	
	default(flagmanager)		
	mode='save'
	versionname='AC_shadow_DV05'
	for i in range(len(svis)):
    		vis=svis[i]
    		flagmanager()
				

#------------------------------#
#------------------------------#
#- Flag data no longer needed -#

if (flag_PtgAtm):
	print "<< Flag pointing and atmospheric measurements\n"
	default(flagdata2)		
	selectdata= T
	intent='*POINTING*'
	manualflag= T
	mf_intent='*POINTING*'
	for i in range(len(svis)):
    		vis=svis[i]
    		flagdata2()

	intent='*ATMOSPHERE*'
	mf_intent='*ATMOSPHERE*'
	for i in range(len(svis)):
    		vis=svis[i]
    		flagdata2()

	print "<< Save current state of flagging as BeforeSplit\n"
	default(flagmanager)	       
	mode='save'
	versionname='BeforeSplit'
	for i in range(len(svis)):
    		vis=svis[i]
    		flagmanager()

#----------------------------#
#----------------------------#
#-Apply WVR and Tsys for FDM-#

# apply WVR per field for fields 1~10 in first 2 data sets, 1~9 in 3rd


if (FDM_apply_wvr_tsys):
	print "<< Apply WVR per field, interpolate Tsys calibration\n"
	default(applycal)
	vis=svis[0]
	interp = 'nearest','nearest'
	gaintable = roots[0]+'.wvr',roots[0]+'.tsys.fdm'
	flagbackup = F
	for f in range(1,11):
		f = str(f)
		field=f
		gainfield=[f,'']
		applycal()
	
	default(applycal)
	vis=svis[1]
	interp = 'nearest','nearest'
	gaintable = roots[1]+'.wvr',roots[1]+'.tsys.fdm'
	flagbackup = F
	for f in range(1,11):
		f = str(f)
		field=f
		gainfield=[f,'']
		applycal()
	
	default(applycal)
	vis=svis[2]
	interp = 'nearest','nearest'
	gaintable = roots[2]+'.tsys.fdm'
	flagbackup = F
	applycal()


	print "<< Applied Tsys and WVR calibration to FDM data\n"
	print "<< Ignor	errors about spw which are not in FDM mode\n" 	
#------------#
#------------#
#-SPLIT  FDM-#

if (split_FDM):	
	default(split)
	spw='1,3,5,7'
	for i in range(len(svis)):
    		vis=svis[i]
    		outputvis=roots[i]+'.fdm_split'
    		split()
	print "<< Created roots[i].fdm_split; spws now numbered 0,1,2,3\n"

#------------------------------#
#------------------------------#
#-- Correct nrao530 position --#

if (correct_phaserefpos):
	default(fixvis)
	field='nrao530'
	reuse = F
	phasecenter='J2000 17h33m02.705775s -13d04m49.54825s'
	for i in range(len(roots)):
    		vis=roots[i]+'.fdm_split'
                fixvis()
	print "<< Corrected position of phase reference source\n"
    	
#----------------------------------------------------------------------------------------------#
#--------------------------- END OF INSTRUMENTAL CALIBRATIONS ---------------------------------#
#----------------------------------------------------------------------------------------------#

# Flag worst atmospheric spikes in all data
# Flag Titan more severely, to ensure accurate flux scale.
# Phase-ref and target have slight remnants of ozone effects
# to be dealt with later
# Flag polarizations known to be bad (if partially flagged, other hand
# is also ignored so we flag both)  

#--------------------------------------#
#--------------------------------------#
#-Flag remaining atmospheric lines etc-#

if (flag_all):
	print "<< Flag worst atmospheric lines \n" 		
	default(flagdata2)
	selectdata = T
	spw='2:550~750,3:2400~2450'
	manualflag= T
	mf_spw='2:550~750,3:2400~2450'
	for i in range(len(roots)):
		vis=roots[i]+'.fdm_split'
		flagdata2()

if (flag_titan):
	print "<< Flag worse lines for Titan, we will only use first ms\n"
	default(flagdata2)
	selectdata = T
	field='Titan'
	spw='1:1307~1370,3:2291~2470'
	manualflag= T
	mf_field='Titan'
	mf_spw='1:1307~1370,3:2291~2470'
	for i in range(len(roots)):
		vis='X5b9.fdm_split'
		flagdata2()

if (flag_badpols):
	print "<< Flag data known to be bad in one or both pols\n"
	default(flagdata2)
	selectdata = T
	manualflag= T
	for i in range(len(roots)):
		vis=roots[i]+'.fdm_split'
		field=''
		antenna='DV02'
		spw='2,3'
		mf_antenna='DV02'
		mf_spw='2,3'
		flagdata2()
		antenna='DV09'
		spw='1'
		mf_antenna='DV09'
		mf_spw='1'
		flagdata2()
		antenna='DV07'
		spw='0'
		mf_antenna='DV07'
		mf_spw='0'
		flagdata2()

#---------TITAN-----------#

# We only use the first observation of Titan, in X5b9.fdm_split, since
# the other scans on Titan (and 3C279) are at low elevation.

#--------------------#
#--------------------#
#-Fix Titan Position-#
if (fix_titan):
	print "<< About to fix position of Titan. Do not interrupt, e.g. quit plotms\n"
	print "<< beforehand. If you abort task, you may have to recreate concatenated MS\n"
	fixplanets('X5b9.fdm_split', 'Titan', True)
	print "<< Titan Position Fixed\n"

#--------------------#
#--------------------#
#--Set TITAN model --#

if (set_jy_titan):
	print "<< About to insert Titan model with setjy\n"
	default(setjy)
	vis='X5b9.fdm_split'
	field='Titan'
	standard='Butler-JPL-Horizons'
	setjy()
	print "<< Titan model set\n"
	
#---------3c279-----------#
	
#---------------------#
#---------------------#
#--Phase Cal 3c279----#

if (phase_bp):
	print "<< About to calibrate 3c279 phase with time\n"
	default(gaincal)
	vis='X5b9.fdm_split'
	caltable='3c297.phcal0'
	field='3c279'
	refant='8'
	calmode='p'
	solint='int'
	minsnr=2.
	minblperant=4
	gaincal()
	print "<< Generated 3c297.phcal0 \n"

if (for_plotms):	
	default(applycal)   # just to check, inspect in plotms
	vis='X5b9.fdm_split'
	gaintable=caltable
	field='3C279'
	interp='nearest'
	applycal()
	print "<< 3c297.phcal0 applied\n"
	
#---------------------------#
#---------------------------#
#--Generate BP Cal 3c279----#

if (gen_BP):
	print"<< About to generate Bandpass table\n"
	default(bandpass)
	vis='X5b9.fdm_split'
	caltable='3c297.bp' 
	field='3c279'
	refant='8'
	gaintable='3c297.phcal0'
	solint='inf'
	combine=''
	solnorm= T
	minblperant=4
	fillgaps=100
	bandpass()
	print "<< Generated bandpass table 3c297.bp\n"

if (for_plotms):		
	default(applycal)   # if checking in plotms wanted
	vis='X5b9.fdm_split'
	gaintable='3c297.phcal0','3c297.bp'
	interp='linear','nearest'
	applycal()
	
#-------------------------------------------------#
#-------------------------------------------------#
#-Generate phase cal table for calibrator sources-#

# First derive calibration for all calibrators in first data set,
# then for the phase-reference source in all data sets.

if (first_phase_cals):
	print "<< About to generate phase cal table for cal sources\n"
	default(gaincal)
	vis='X5b9.fdm_split'
	caltable='X5b9.phcal1'  #Phase cal table.
	field='Titan,3c279,nrao530'
	refant='8'
	calmode='p'
	solint='int'                      # per integration Bright sources 
	gaintable='3c297.bp'		  # BP cal table
	interp='nearest'                  # to extrapolate BP table
	minsnr=2.
	minblperant=4
	gaincal()

	default(gaincal)
	for i in range(len(roots)):
		vis=roots[i]+'.fdm_split'
		caltable=roots[i]+'_cals.phcal1'  #Phase cal table.
		field='nrao530'
		refant='8'
		calmode='p'
		solint='int'                      # per integration Bright sources 
		gaintable='3c297.bp'		  # BP cal table
		interp='nearest'                  # to extrapolate BP table
		minsnr=2.
		minblperant=4
		gaincal()
	print "<< Generated roots[i]_cals.phcal1 for cal sources\n"

if (do_plots):
	default(plotcal)
	for i in range(len(roots)):
		caltable=roots[i]+'_cals.phcal1'
		xaxis='time'
		yaxis='phase'
		plotrange=[0,0,-180,180]
		field='nrao530'
		subplot=421
		plotcal()
	
#-----------------------------------------------------#
#-----------------------------------------------------#
#-Generate phase cal table to apply to target sources-#
# This uses a longer solution interval to improve interpolation
# across the target

if (first_phase_target):
	print "<< About to generate phase cal table for target sources\n"
	default(gaincal)
	for i in range(len(roots)):
		vis=roots[i]+'.fdm_split'
		caltable=roots[i]+'_forSgrA.phcal1'  #Phase cal table.
		field='nrao530'			  
		refant='8'
		calmode='p'
		solint='inf'                      # per-scan for interpolation to target 
		gaintable='3c297.bp'		  # BP cal table
		interp='nearest'                  # to extrapolate BP table
		minsnr=2.
		minblperant=4
		gaincal()
	print "<< Generated roots[i]_forSgrA.phcal1 to apply to target sources\n"

if (do_plots):
	default(plotcal)
	xaxis='time'
	yaxis='phase'
	subplot=421
	plotrange=[0,0,-180,180]
	field='nrao530'
	for i in range(len(roots)):
		caltable=roots[i]+'_forSgrA.phcal1' 
		plotcal()
	
#--------------------------------#
#--------------------------------#
#-Use Titan to derive flux scale-#

# Derive gain solutions for calibration sources and compare the
# solutions for the other sources with Titan to derive the other flux
# densities. 3c279 is not needed here but serves as a check. 

# spw 1 and 2 are badly affected by atmospheric lines so we extrapolate
# from the other spw

if (bootstrap):
	default(gaincal)
	vis='X5b9.fdm_split'
	caltable='X5b9.ampcal0'
	field='1,2,3'
	solint='inf'
	refant='8'
	calmode='ap'
	gaintable='3c297.bp','X5b9.phcal1'
	interp='nearest','linear'
	minsnr=2.
	minblperant=4
	gaincal()
	print "<< Generated X5b9.ampcal0 with amp solutions for calibrators\n"
	
	print "<< About to derive flux scale using X5b9.ampcal0\n"
	default(fluxscale)
	vis='X5b9.fdm_split'
	caltable='X5b9.ampcal0'
	reference='1'
	refspwmap=[0,0,3,3]		#to avoid using spw 1,2 
	fluxtable='X5b9.fluxcal'
	fluxscale()
	print "<< Check logger for derived fluxes of nrao530 and 3c279\n"
	print "<< Values for spw 1,2 are less accurate due to atmospheric lines\n"

# extract derived flux densities for nrao530 from logger. Use spw
# frequencies from listobs to calculate spectral index and interpolate
# values for spw 1 and 2.
        print "<< Extract nrao530 flux density and derive spectral index\n" 
        print "<< from spw 0, 3 to interplate into spw 1, 2\n"	
        freqs=[217137.156, 218939.506, 230960.494, 232762.844]
        w0='nrao530 in SpW=0'
	w3='nrao530 in SpW=3'
	for l in open('casapy.log'):
		if (re.match(w0, l[57:73])):
			f0=float(l[78:84])
		if (re.match(w3, l[57:73])):
			f3=float(l[78:84])
	alpha=log(f0/f3)/log(freqs[0]/freqs[3])		
	f1=f0*(freqs[1]/freqs[0])**alpha
	f2=f0*(freqs[2]/freqs[0])**alpha
        print "<< Flux density of nrao530 in spw 0, 1, 2, 3 (expect ~2.7 Jy)\n"
	print f0, f1, f2, f3
	print "<< Spectral index", alpha

#--------------------------------#
#--------------------------------#
#-Set flux densities for nrao530-#

if (set_phaseref_flux):      
	print "<< Set phase ref flux density model in all datasets\n"

        default(setjy)
	field='nrao530'
	for i in range(len(roots)):
		vis=roots[i]+'.fdm_split'
		spw='0'
		fluxdensity=[float(f0),0,0,0]
		setjy()
		spw='1'
		fluxdensity=[float(f1),0,0,0]
		setjy()
		spw='2'
		fluxdensity=[float(f2),0,0,0]
		setjy()
		spw='3'
		fluxdensity=[float(f3),0,0,0]
		setjy()
	print "nrao530 flux densities set\n"			

#----------------------------------#
#----------------------------------#
#-Generate phase & amp cal tables -#

if (amp_cals):
	print "<< About to generate amp cal tables for phase ref\n"
	default(gaincal)
	field='nrao530'
	refant='8'
	calmode='ap'
	solint='inf'                      	
	interp='nearest','linear'
	minsnr=2.
	minblperant=4
	
	for i in range(len(roots)):
		vis=roots[i]+'.fdm_split'
		caltable=roots[i]+'_nrao530.acal1'
		gaintable='3c297.bp',roots[i]+'_cals.phcal1'  #BP & phase cal table
	       	gaincal()
	print "<< Generated roots[i]_nrao530.acal1 from phase ref\n"

if (do_plots):
	default(plotcal)
	xaxis='time'
	yaxis='amp'
	iteration='antenna'
	field='3'
	subplot=421
	plotrange=[0,0,0.15,0.35]		
	for i in range(len(roots)):
		caltable=roots[i]+'_nrao530.acal1'
		plotcal()
	
#----------------------------------#
#----------------------------------#
#-Apply all calibration to nrao530-#
# Can skip if you intend to do secondary bandpass and baseline
# calibration

if (apply_cal_cals):
	print "<< About to apply  Cal tables to phase reference\n"
	interp='nearest','linear','linear'
	field='nrao530'
	default(applycal)
	for i in range(len(roots)):
		vis=roots[i]+'.fdm_split'
		gaintable='3c297.bp',roots[i]+'_cals.phcal1',roots[i]+'_nrao530.acal1'
		applycal()
	print "<< Calibration applied to nrao530"

#--------------------------------------#
#--------------------------------------#
#---Apply all calibration to target----#
# Can skip if you intend to do secondary bandpass and baseline
# calibration

if (apply_cal_target):
	print "<< About to apply Cal tables to Target Sources\n"
	interp='nearest','linear','linear'
	field='SgrA'
	default(applycal)
	for i in range(len(roots)):
		vis=roots[i]+'.fdm_split'
		gaintable='3c297.bp',roots[i]+'_forSgrA.phcal1',roots[i]+'_nrao530.acal1'
		applycal()
	print "<< Calibration applied to SgrA*"
	
#------------------------------------------#
#------------------------------------------#
#-Generate BP table from phase cal nrao530-#
# Inspection shows that 3c279 was observed at lower elevation than the target/phase-ref
# Hence, the bandpass corrections from 3c279 over-compensate for atmospheric lines
# in the target and phase-ref
	
if (phase_BP):
	print "<< Generating BP table from phase ref due to atmospheric lines in 3c279\n"
	default(bandpass)
	field='nrao530'
	refant='8'
       	solint='inf'
	combine=''
	solnorm= T
	minblperant=4
	fillgaps=100
	for i in range(len(roots)):
		vis=roots[i]+'.fdm_split'
		caltable=roots[i]+'_nrao530.bp' 
		gaintable='3c297.bp',roots[i]+'_cals.phcal1',roots[i]+'_nrao530.acal1'
		bandpass()
	print "<< Generated roots[i]_nrao530.bp using phase ref source\n"

#-------------------------------------#
#-------------------------------------#
#---blcal on phase-reference source---#
# Experiment shows that this greatly improves snr on the phase ref.

if (phref_blcal):
	print "<< Generating baseline-dependent calibration using nrao530\n"
	default(blcal)
	field='nrao530'
	combine=''
	solnorm= T
	for i in range(len(roots)):
		vis=roots[i]+'.fdm_split'
		gaintable='3c297.bp',roots[i]+'_cals.phcal1',roots[i]+'_nrao530.acal1',roots[i]+'_nrao530.bp'
		caltable=roots[i]+'_nrao530.bl'
		blcal()
	print "<< Generated roots[i]_nrao530.bl using phase ref source\n"

#----------------------------------#
#----------------------------------#
#-Apply all calibration to nrao530-#

if (apply_cal_cals_2):
	print "<< About to apply all cal tables to phase ref\n"
	default(applycal)
	field='nrao530'
	interp='nearest','linear','linear','nearest','linear'
	for i in range(len(roots)):
		vis=roots[i]+'.fdm_split'
		gaintable='3c297.bp',roots[i]+'_cals.phcal1',roots[i]+'_nrao530.acal1',roots[i]+'_nrao530.bp',roots[i]+'_nrao530.bl'
		applycal()
	print "<< Calibration applied to nrao530"

#---------------------------------------#
#---------------------------------------#
#--- Split and concatenate phase ref ---#

if (phsref_split):
	print "<< Split out nrao530, average data and drop pointing tables\n"
        default(split) 
        field='nrao530'
	timebin='30s'
	width=20
	copypointing = F
        for i in range(len(roots)):
	        vis=roots[i]+'.fdm_split'
		outputvis=roots[i]+'nrao530.split'
                split()

        default(concat)
	vis= glob('*nrao530.split')
	concatvis='nrao530_B6.ms'
	concat()
	print "<< Split out phase ref, averaged every 30s, 20 channels\n"
	print "<< Concatenated into nrao530_B6.ms\n"
	
#------------------------#
#------------------------#
#--- Image phase ref ---#

if (phsref_image):
	default(clean)
	vis='nrao530_B6.ms'
	field='nrao530'
	imagename = 'nrao530.clean'
	psfmode='hogbom'
	interactive = F
	cell=['0.25arcsec']
	threshold='0.15mJy'
	mask = [121, 121, 136, 136]
	clean()
        print "<< Imaged phase ref, created nrao530.clean.image\n"

	default(imstat)
	imagename  =  'nrao530.clean.image'
	box= '21,21,51,200'
	nraostats=imstat()
	nraorms=nraostats['rms']
	box= '101,101,156,156'
	nraostats=imstat()
	nraomax=nraostats['max']
	snr = nraomax/nraorms
	print 'nrao530 snr', snr[0], 'peak', nraomax[0], 'rms', nraorms[0]

# I get rms ~1 mJy/beam.  Sensitivity calculator gives 0.15 mJy
# for 220 GHz, 10 antennas, 7 GHz, 12 min on-source.  Ouch!
# Dynamic range limited? Dynamic range is ~2000
# also residual atmospheric lines

#----------------------------------------#
#----------------------------------------#
#--- Apply all calibration to target ----#

if (apply_cal_target_2):
	print "<< About to apply all Cal tables to Target Sources\n"
	interp='nearest','linear','linear','nearest','linear'
	field='SgrA'
	default(applycal)
	for i in range(len(roots)):
		vis=roots[i]+'.fdm_split'
		interp='nearest','linear','linear','nearest','linear'
		gaintable='3c297.bp',roots[i]+'_forSgrA.phcal1',roots[i]+'_nrao530.acal1',roots[i]+'_nrao530.bp',roots[i]+'_nrao530.bl'
		applycal()
	print "<< Calibration applied to Targets"

#------------------------------------#
#------------------------------------#
#--- Split and concatenate target ---#

if (target_split):
	print "<< Split out SgrA, average data and drop pointing tables\n"
        default(split) 
        field='SgrA'
	timebin=averagetime
	width=averagechan
        for i in range(len(roots)):
	        vis=roots[i]+'.fdm_split'
		outputvis=roots[i]+'SgrA.split'
                split()
	print "<< Split out target fields, averaged in time and channel by\n"
        print averagetime, averagechan

        default(concat)
	vis= glob('*SgrA.split')
	concatvis='SgrA_B6.ms'
	copypointing = F
	concat()
	print "<< Concatenated into SgrA_B6.ms\n"

	vis='SgrA_B6.ms'
	listobs()
	print "<< See logger for new spw (0~3) and fields (0~6)\n"


#--------------------------------------------------------------------------------#
#--------------------------END OF STANDARD CALIBRATIONS--------------------------#
#--------------------------------------------------------------------------------#

# See separate script for imaging