########################################################################
# Data reduction script for CenA Band 6:
# - Calibration script: "CenA_Band6_Calibration.py" -
# Tested in CASA Version 3.3.0 (r16856)
########################################################################
"""
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=T

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

#----- Input MS names 
#-----  (the data provided have already been coverted to CASA measurement set
#-----  format (.ms) using the task importasdm)
basename_ORIG=['uid___A002_X2762c7_X132','uid___A002_X2762c7_X26','uid___A002_X2762c7_X23e','uid___A002_X2762c7_X34a','uid___A002_X2762c7_X456','uid___A002_X2762c7_X562']

#----- Re-name the MS to something shorter
basename_all=['X132','X26','X23e','X34a','X456','X562']

for i in range(len(basename_ORIG)):
	os.system('mv '+ basename_ORIG[i]+'.ms '+ basename_all[i]+'.ms')
	os.system('mv '+ basename_ORIG[i]+'.ms.flagversions '+ basename_all[i]+'.ms.flagversions')

#----- List of field names to be used
field_names = ['Titan','3c279','J1427-421','Cent*']

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

#----- Create a summary text file for each dataset
for asdm in basename_all:
	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_all:
        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")

#----------------------------------------------------------------------------------
#----- Fix planet/satellite position ----------------------------------------------
#----------------------------------------------------------------------------------

for asdm in basename_all:
	fixplanets(vis=asdm+'.ms', field='Titan', fixuvw=True)

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

#----- Inspect amplitude versus time for each dataset
if(interact):
	for asdm in basename_all:
		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")

#----------------------------------------------------------------------------------
#----- Initial Flagging -----------------------------------------------------------
#----------------------------------------------------------------------------------

for asdm in basename_all:
#----- Remove all the previous flags
	flagdata(vis=asdm+'.ms',mode='manualflag', unflag=T, flagbackup=F)
#----- Flagging shadowed data
	flagdata(vis=asdm+'.ms',mode = 'shadow', diameter=12.0, flagbackup = F)
#----- Flagging pointing and atmosphere calibration scans
	flagdata(vis=asdm+'.ms', mode='manualflag', intent='*POINTING*',flagbackup = F)
	flagdata(vis=asdm+'.ms', mode='manualflag', intent='*ATMOSPHERE*',flagbackup = F)
#----- Flagging autocorrelation data
	flagdata(vis=asdm+'.ms',autocorr=True,flagbackup=F)
#----- Saving 'a priori' flags
	flagmanager(vis = asdm+'.ms', mode = 'save', versionname = 'Apriori')
#----- Restoring up 'a priori' flags"+asdm
#        flagmanager(vis = asdm+'.ms', mode = 'restore', versionname = 'Apriori')

#----------------------------------------------------------------------------------
#----- Examine and Apply Tsys and WVR calibration tables --------------------------
#----------------------------------------------------------------------------------

if(calplots):
	for asdm in basename_all:
#----- Plotting Tsys vs. time
		plotcal(caltable=asdm+'.tsys.cal.fdm', 
			xaxis="time",yaxis="amp",
			spw='1:1200~1200',plotsymbol=".", subplot=421,
			antenna='0~7',
			iteration='antenna', figfile=asdm+'.tsys_vs_time.page1.png',
			fontsize=6.0)    
		plotcal(caltable=asdm+'.tsys.cal.fdm', 
			xaxis="time",yaxis="amp",
			antenna='8~15',
			spw='1:1200~1200',plotsymbol=".", subplot=421,
			iteration='antenna', figfile=asdm+'.tsys_vs_time.page2.png',
			fontsize=6.0)    
#----- Plotting Tsys vs. frequency
		plotcal(caltable=asdm+'.tsys.cal.fdm',
			xaxis="freq",yaxis="amp",
			spw='', plotsymbol=".", subplot=421,
			iteration='antenna', figfile=asdm+'.tsys_vs_freq.page1.png',
			antenna='0~7', fontsize=6.0)    
		plotcal(caltable=asdm+'.tsys.cal.fdm',
			xaxis="freq",yaxis="amp",
			spw='', plotsymbol=".", subplot=421,
			iteration='antenna', figfile=asdm+'.tsys_vs_freq.page2.png',
			antenna='8~15', fontsize=6.0)

# In all datasets:  DV12 has low tsys values in polarization XX at spw = 7.
# DV08 has higher Tsys in spw =5, pol YY.
# DV02 has unrealistic high tsys measures in spw 5, both pol, for the last four datasets.
# See "Further Data Flagging" section below.
    
#----- Apply Tsys, WVR calibrations and correct antenna positions
for asdm in basename_all:
	for field in field_names:
		applycal(vis=asdm+".ms", spw='1,3,5,7', 
			 field=field, gainfield=[field,field],
			 interp='nearest', 
			 gaintable=[asdm+".tsys.cal.fdm", asdm+'.wvr.cal.smooth'],
			 flagbackup=F, calwt=T)

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

#----- Plot spectral plot in amplitude and phase
if(interact):
	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):
	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')

#----------------------------------------------------------------------------------
#----- Split out corrected data and save the flags ------------
#----------------------------------------------------------------------------------

#----- First save the flags
for asdm in basename_all:
	flagmanager(vis = asdm+'.ms', mode = 'save', versionname = 'Before_flagging')

#----- Now do the split and output another summary file
for asdm in basename_all:
	os.system('rm -rf '+ asdm+'.wvrtsys.ms')
	split(vis=asdm+'.ms', outputvis=asdm+'.wvrtsys.ms', 
	      datacolumn='corrected', spw='1,3,5,7', keepflags=F)
	os.system('rm '+asdm+'.wvrtsys.listobs.txt')
	listobs(vis=asdm+'.wvrtsys.ms', listfile=asdm+'.wvrtsys.listobs.txt', verbose=True)

#----- Remove the pointing tables
for asdm in basename_all:
	tb.open(asdm+'.wvrtsys.ms/POINTING',nomodify=False)
	a = tb.rownumbers()
	tb.removerows(a)
	tb.close()

#----------------------------------------------------------------------------------
#----- Further Data Flagging ------------------------------------------------------
#----------------------------------------------------------------------------------

#----- Flag shadowed data
for asdm in basename_all:
	flagdata(vis=name+'.wvrtsys.ms', flagbackup = F, mode = 'shadow')

#----- Flag for the Tsys issues identified above
for asdm in basename_all:
	flagdata2(flagbackup = F, vis=asdm+'.wvrtsys.ms',manualflag=T,mf_antenna='DV12',spw='3',correlation='XX', selectdata=True)
	flagdata2(flagbackup = F, vis=asdm+'.wvrtsys.ms',manualflag=T,mf_antenna='DV08',spw='2',correlation='YY', selectdata=True)

for asdm in ['X23e','X34a','X456','X562']:
	flagdata2(flagbackup = F, vis=asdm+'.wvrtsys.ms',manualflag=T,mf_antenna='DV02', selectdata=True)

#----- Flag the less sensitive edge channels of every data set
for asdm in basename_all:
	flagdata2(vis = asdm+'.wvrtsys.ms',manualflag=T, mf_spw = ['*:0~7;3831~3839'],  selectdata=True, flagbackup = F)

#----- Flag some birdies
for asdm in basename_all:
	flagdata2(vis = asdm+'.wvrtsys.ms', manualflag=T, mf_spw='0~3:320~322;447~448;720~721;1555~1556;2275~2276;2847~2848', selectdata=True, flagbackup = F)

#----- Flag possible contamination of Titan by Saturn
for asdm in basename_all:
	flagdata2(vis=asdm+'.wvrtsys.ms', manualflag=T, field='Titan', mf_uvrange='0~50',  selectdata=True, flagbackup = F)

for asdm in basename_all:
	flagmanager(vis=asdm+'.wvrtsys.ms',mode='save',versionname ='User')

#----------------------------------------------------------------------------------
#----- Bin to lower spectral res (~10 km/s) ---------------------------------------
#----------------------------------------------------------------------------------

for asdm in basename_all:
	os.system('rm -rf '+asdm+'.wvrtsys.res10.ms')
	split(vis = asdm+'.wvrtsys.ms',outputvis = asdm+'.wvrtsys.res10.ms', width=15, datacolumn='data')

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

#----------------------------------------------------------------------------------
#----- Calibration Steps ----------------------------------------------------------
#----------------------------------------------------------------------------------

#----- Running a short solution interval phase calibration for bandpass
for asdm in basename_all:
	os.system('rm -rf '+asdm+'.bpphase.gcal')
	gaincal(vis = asdm+'.wvrtsys.res10.ms',
		selectdata=T,field = '3c279',spw = '*:73~87',
		caltable = asdm+'.bpphase.gcal',
		solint = 'int',refant = 'DV11',calmode='p')

#----- Running a bandpass calibration
for asdm in basename_all:
	os.system('rm -rf '+asdm+'.bandpass.bcal')
	bandpass(vis = asdm+'.wvrtsys.res10.ms',
		 field = '3c279',
		 gaintable = asdm+'.bpphase.gcal',
		 caltable = asdm+'.bandpass.bcal',
		 bandtype='B',
		 solint = 'inf',combine = 'scan,field', solnorm=T,refant = 'DV11',
		 minblperant=3,minsnr=2,fillgaps=1)

#----- Plotting bandpass solutions
if(calplots):
	for asdm in basename_all:
		plotcal(caltable = asdm+'.bpphase.gcal',
			xaxis = 'time', yaxis = 'phase',
			iteration = 'antenna', plotrange=[0,0,-180,180],
			showgui=False, subplot=421, figfile=asdm+'.bpphase.page1.png',
			antenna='0~7')
		plotcal(caltable = asdm+'.bpphase.gcal',
			xaxis = 'time', yaxis = 'phase',
			iteration = 'antenna', plotrange=[0,0,-180,180],
			showgui=False, subplot=421, figfile=asdm+'.bpphase.page2.png',
			antenna='8~15')
		plotcal(caltable = asdm+'.bandpass.bcal', 
			xaxis = 'freq',yaxis = 'amp',
			plotrange = [0,0,0.8,1.2],
			antenna='0~7', iteration='antenna',
			showgui=False, subplot=421, figfile=asdm+'.bcal_amp.page1.png')
		plotcal(caltable = asdm+'.bandpass.bcal', 
			xaxis = 'freq',yaxis = 'amp',
			plotrange = [0,0,0.8,1.2],
			antenna='8~15', iteration='antenna',
			showgui=False, subplot=421, figfile=asdm+'.bcal_amp.page2.png')

if(calplots):
	for asdm in basename_all:
		plotcal(caltable = asdm+'.bandpass.bcal', 
			xaxis = 'freq',yaxis = 'phase', iteration='antenna',
			antenna='0~7', subplot=421, figfile=asdm+'.bcal_phase.page1.png',
			plotrange = [0,0,-180,180], showgui=False)
		plotcal(caltable = asdm+'.bandpass.bcal', 
			xaxis = 'freq',yaxis = 'phase', iteration='antenna',
			antenna='8~15', subplot=421, figfile=asdm+'.bcal_phase.page2.png',
			plotrange = [0,0,-180,180], showgui=False)

#----- Reading model for Titan (Setjy)
for asdm in basename_all:
	setjy(vis = asdm+'.wvrtsys.res10.ms',field = 'Titan',
	      standard = 'Butler-JPL-Horizons 2010')

#----- Carrying out short timescale (integration, ~ 8 sec) phase solution
for asdm in basename_all:
	os.system('rm -rf '+asdm+'.intphase.gcal')
	gaincal(vis=asdm+'.wvrtsys.res10.ms',
		gaintable=asdm+'.bandpass.bcal', 
                caltable=asdm+'.intphase.gcal',
                calmode='p',
                field='Titan,3c279,J1427*',
                spw='*:5~250',
                refant='DV11', solint='int',minsnr=2.0,minblperant=4)

#----- Carrying out longer timescale (by scan) phase solution
for asdm in basename_all:
	os.system('rm -rf '+asdm+'.scanphase.gcal')
	gaincal(vis=asdm+'.wvrtsys.res10.ms',
                gaintable=asdm+'.bandpass.bcal', 
                caltable=asdm+'.scanphase.gcal',
                calmode='p',
                field='Titan,3c279,J1427*',
                spw='*:5~250',
                refant='DV11', solint='Inf',minsnr=2.0,minblperant=4)

#----- Solving for longer (scan) interval amplitude solution
for asdm in basename_all:
	os.system('rm -rf '+asdm+'.amp.cal')
	gaincal(vis = asdm+'.wvrtsys.res10.ms',
               gaintable =[asdm+'.bandpass.bcal',asdm+'.intphase.gcal'],
               caltable = asdm+'.amp.cal',
               calmode='a',
               field = 'Titan,3c279,J1427*',
               spw='*:5~250',
               refant = 'DV11',solint = 'inf')

#----- Scaling amplitude calibration to match Titan (except for spw=2,
#-----  where the co line is. spw2 is referenced to spw=3)
for asdm in basename_all:
	os.system('rm -rf '+asdm+'.flux.cal')
	fluxscale(vis = asdm+'.wvrtsys.res10.ms',
                caltable = asdm+'.amp.cal',
                fluxtable = asdm+'.flux.cal',
                reference = 'Titan',
	        refspwmap=[0,1,3,3],
                transfer = '3c279,J1427*')

#----- Plotting solutions
if(calplots):
	for asdm in basename_all:
		plotcal(caltable = asdm+'.scanphase.gcal',
			xaxis = 'time', yaxis = 'phase',
			iteration = 'antenna', plotrange=[0,0,-180,180],
			showgui=True, subplot=421, figfile=asdm+'.scanphase.page1.png',
			antenna='0~7', fontsize=6.0)
		plotcal(caltable = asdm+'.scanphase.gcal',
			xaxis = 'time', yaxis = 'phase',
			iteration = 'antenna', plotrange=[0,0,-180,180],
			showgui=False, subplot=421, figfile=asdm+'.scanphase.page2.png',
			antenna='8~15', fontsize=6.0)
		plotcal(caltable = asdm+'.flux.cal', 
			xaxis = 'time',yaxis = 'amp',
			plotrange = [0,0,0,0],
			antenna='0~7', iteration='antenna',
			showgui=False, subplot=421, figfile=asdm+'.flux.page1.png',
			fontsize=6.0)
		plotcal(caltable = asdm+'.flux.cal', 
			xaxis = 'time',yaxis = 'amp',
			plotrange = [0,0,0,0],
			antenna='8~15', iteration='antenna', fontsize=6.0,
			showgui=False, subplot=421, figfile=asdm+'.flux.page2.png')

#----- Applying calibrations to the different sources:

for asdm in basename_all:
	applycal(vis=asdm+'.wvrtsys.res10.ms',field='3c279',
		gaintable=[asdm+'.bandpass.bcal',asdm+'.intphase.gcal',asdm+'.flux.cal'],
		interp=['nearest','nearest','nearest'],
		gainfield=['3c279','3c279','3c279'],flagbackup=F, calwt=F)
	applycal(vis=asdm+'.wvrtsys.res10.ms',field='Titan',
		gaintable=[asdm+'.bandpass.bcal',asdm+'.intphase.gcal',asdm+'.flux.cal'],
		interp=['nearest','nearest','nearest'],
		gainfield=['3c279','Titan','Titan'],flagbackup=F, calwt=F)
	applycal(vis=asdm+'.wvrtsys.res10.ms',field='J1427*',
		gaintable=[asdm+'.bandpass.bcal',asdm+'.intphase.gcal',asdm+'.flux.cal'],
		interp=['nearest','nearest','nearest'],
		gainfield=['3c279','J1427*','J1427*'],flagbackup=F, calwt=F)
	applycal(vis=asdm+'.wvrtsys.res10.ms',field='Cent*',
		interp=['nearest','linear','linear'],
		gaintable=[asdm+'.bandpass.bcal',asdm+'.scanphase.gcal',asdm+'.flux.cal'],
		gainfield=['3c279','J1427*','J1427*'],flagbackup=F, calwt=F)

#----- Inspecting calibrated data
if(interact):
	for asdm in basename_all:
		plotms(vis = asdm+'.wvrtsys.res10.ms', 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+'.wvrtsys.res10.ms', 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+'.wvrtsys.res10.ms', 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+'.wvrtsys.res10.ms', xaxis='time', yaxis='phase',
		        ydatacolumn='corrected', field='',
			avgchannel='3840', avgscan=F, avgbaseline=F, coloraxis='corr')

#----------------------------------------------------------------------------------
#----- Split and concatenate calibrated data on CenA ------------------------------
#----------------------------------------------------------------------------------

#----- Split out the calibrated data on CenA
for asdm in basename_all:
	os.system('rm -rf '+asdm+'.cal.ms')
	split(vis = asdm+'.wvrtsys.res10.ms',outputvis = asdm+'.cal.ms',
	      datacolumn='corrected', field = 'Cent*',spw='', keepflags=False)

#----- Create a summary file for the calibrated datasets
for asdm in basename_all:
	os.system('rm '+asdm+'.cal.listobs.txt')
	listobs(asdm+'.cal.ms',listfile=asdm+'.cal.listobs.txt')

#----- Concatenate into one dataset
cal_vis = [vis+'.cal.ms' for vis in basename_all]
os.system('rm -rf CenA.cal.ms')
concat(vis=cal_vis, concatvis='CenA.cal.ms', timesort=T)

#----- rebin the data with time averaging
split(vis='CenA.cal.ms', outputvis="CenA.cal.rebin.ms",datacolumn="data",timebin="60s")

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