##############################################################################
# Data reduction script for VY CMa Band 7 658-GHz 12-m data
# Tested using CASA Version 4.2.2
#
# Datasets used:
#
# uid___A002_X6d34fd_Xb3d - observed Aug 18 2013 with 18 antennas
# Flux scale calibrator J0522-364
# Bandpass calibrator J0522-364
# Phase reference source J0648-3044
# Target Vy CMa
#
#
##############################################################################
#----------------------------------------------------------------------------------
#----- Optional Steps -------------------------------------------------------------
#----------------------------------------------------------------------------------
#----- Set this option to true if you want to make diagnostic plots along the way.
#----- Note that this will slow down the reduction somewhat.
calplots=T
#----------------------------------------------------------------------------------
#----- Some setup steps -----------------------------------------------------------
#----------------------------------------------------------------------------------
version = casalog.version()
print "You are using " + version
if (version < '4.2.0'):
print "\033[91m YOUR VERSION OF CASA IS TOO OLD FOR THIS GUIDE."
print "\033[91m PLEASE UPDATE IT BEFORE PROCEEDING."
else:
print "Your version of CASA is appropriate for this guide."
# Repeat imports if you exit and restart CASA if using copy and paste
#import re
#import os
thesteps=[]
step_title = {0: 'Listing the data and plotting antenna table',
1: 'Flag autocorrelations and pointing scans; back up flags',
2: 'Generation and time averaging of the WVR cal table',
3: 'Generation of the Tsys cal table',
4: 'Generation of the antenna position cal table',
5: 'Plot raw data',
6: 'Application of the WVR, Tsys and antpos cal tables',
7: 'Split out science SPWs and time average',
8: 'Listobs, clear pointing table, and save original flags',
9: 'Plot data with instrumental calibration',
10: 'Initial flagging',
11: 'Bandpass calibration',
12: 'Calibrate the flux scale',
13: 'Gain calibration',
14: 'Application of the bandpass and gain cal tables',
15: 'Plot calibrated phase ref data',
16: 'Split out the science data'}
# The Python variable 'mysteps' will control which steps
# are executed when you start the script using
# execfile('scriptForCalibration.py')
# e.g. setting
# mysteps = [2,3,4]# before starting the script will make the script execute
# only steps 2, 3, and 4
# Setting mysteps = [] will make it execute all steps.
try:
print 'List of steps to be executed ...', mysteps
thesteps = mysteps
except:
print 'global variable mysteps not set.'
if (thesteps==[]):
thesteps = range(0,len(step_title))
print 'Executing all steps: ', thesteps
######################################################
# Name of the asdm
basename='uid___A002_X6d34fd_Xb3d'
# We have already imported the asdm, so we skip the importasdm
# step that you would normally have to do
# We list the data and plot the antenna table
mystep = 0
if(mystep in thesteps):
casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
print 'Step ', mystep, step_title[mystep]
#os.system('rm -rf '+basename+'.ms')
#os.system('rm -rf '+basename+'.ms.flagversions')
#importasdm(asdm=basename, vis=basename+'.ms', process_flags=F,
# asis='Antenna Station CalWVR Receiver Source')
# Do listobs and write info to a text file
os.system('rm '+basename+'_listobs.txt')
listobs(vis=basename+'.ms', listfile=basename+'_listobs.txt', verbose=True)
print "<< Printed listobs output to "+basename+'.listobs.txt'
# Plot antenna configuration to a file
os.system('rm -rf plotants_'+basename+'.png')
plotants(basename+'.ms', figfile='plotants_'+basename+'.png')
print 'Choose refant - with lots of short baselines but not so close other antennas to\nthat it will be shadowed. Definition is at top of script, change if desired.'
# PROBLEM 1 #
##########################################################
############## CHOOSE YOUR REFERENCE ANTENNA
##########################################################
refantenna='DV05'
#################################
# Flag autocorrelations and pointing scans; back up flags
mystep = 1
if(mystep in thesteps):
casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
print 'Step ', mystep, step_title[mystep]
# Remove all previous flags
flagdata(vis=basename+'.ms',mode='unflag', flagbackup=F)
flagdata(vis = basename+'.ms',
mode = 'manual',
spw = '1~12',
autocorr = T,
flagbackup = F)
flagdata(vis = basename+'.ms',
mode = 'manual',
intent = '*POINTING*,*SIDEBAND_RATIO*,*ATMOSPHERE*',
flagbackup = F)
###################################################
# SAVE THE FLAGGING STATE
###################################################
# Store current flagging state
flagmanager(vis = basename+'.ms', mode = 'save',
versionname = 'Apriori')
###############
# Generation and time averaging of the WVR cal table
mystep = 2
if(mystep in thesteps):
casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
print 'Step ', mystep, step_title[mystep]
os.system('rm -rf ' +basename+'.ms.wvr')
mylogfile = casalog.logfile()
casalog.setlogfile(basename+'.ms.wvrgcal')
wvrgcal(vis = basename+'.ms',
caltable = basename+'.ms.wvr',
toffset = 0,
tie = ['Vy CMa,J0648-3044'],
statsource = 'Vy CMa')
casalog.setlogfile(mylogfile)
os.system('rm -rf ' +basename+'.ms.wvr.smooth')
smoothcal(vis = basename+'.ms',
tablein = basename+'.ms.wvr',
caltable = basename+'.ms.wvr.smooth',
smoothtype = 'mean',
smoothtime = 6.048)
if(calplots):
os.system('rm -rf '+basename+'.ms.wvr.smooth*.png')
for ant in range(3):
plotcal(caltable = basename+'.ms.wvr.smooth',
xaxis = 'time', yaxis = 'phase',
plotrange = [0,0,-180,180],
antenna=str(ant*8)+'~'+str(ant*8+7),
iteration = 'antenna',
markersize=3.0,
figfile=basename+'.ms.wvr.smooth_ant'+str(ant*8)+'-'+str(ant*8+7)+'.png',
subplot = 421)
# Flag WVR autocorrelations
flagdata(vis = basename+'.ms',
mode = 'manual',
spw = '0',
autocorr = T,
flagbackup = F)
#####################################
# Generation of the Tsys cal table
mystep = 3
if(mystep in thesteps):
casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
print 'Step ', mystep, step_title[mystep]
os.system('rm -rf '+basename+'.ms.tsys')
gencal(vis = basename+'.ms',
caltable = basename+'.ms.tsys',
caltype = 'tsys')
if(calplots):
os.system('rm -rf '+basename+'.ms.tsys*.png')
for ant in range(3):
plotcal(caltable = basename+'.ms.tsys',
xaxis = 'freq', yaxis = 'tsys',
spw='9:8~120',
markersize=3.0,
antenna=str(ant*8)+'~'+str(ant*8+7),
iteration = 'antenna',
figfile=basename+'.ms.tsys_ant'+str(ant*8)+'-'+str(ant*8+7)+'.png',
subplot = 421)
# PROBLEM 2 #
##########################################################
########## CAN YOU IDENTIFY THE OUTLYING TSYS POINTS (~3000K)
########## Do you want to flag them?
##########################################################
######################################
# Generation of the antenna position cal table
mystep = 4
if(mystep in thesteps):
casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
print 'Step ', mystep, step_title[mystep]
os.system('rm -rf '+ basename+'.ms.antpos')
os.system('rm -rf uid___A002_X6d34fd_Xb3d.ms.antpos')
gencal(vis = 'uid___A002_X6d34fd_Xb3d.ms',
caltable = 'uid___A002_X6d34fd_Xb3d.ms.antpos',
caltype = 'antpos',
antenna = 'DA41,DA46,DA52,DA56,DV03,DV04,DV05,DV07,DV11,DV13,DV14,DV15,DV17,DV18,DV19,DV22,DV24,DV25',
parameter = [-0.018216, 0.067084, 0.060448, 0.000006, 0.000228, -0.000297,0.000046, 0.000004, 0.000165,0.000065, 0.000035, 0.000115, -0.000196, -0.000061, 0.000178, -0.000132, 0.000341, -0.000190, 0.000068, 0.000267, -0.000494, -0.000524, 0.000838, 0.000626, 0.000112, -0.000116, 0.000112, -0.000032, 0.000188, -0.000261, 0.000043, -0.000173, 0.000008, 0.000037, 0.000468, -0.000484, 0.000167, -0.000474, 0.000167, -0.000047, 0.000177, -0.000151, -0.000282, -0.000304, -0.000094, 0.000009, 0.000227, -0.000283, -0.000404, 0.003222, 0.001466, -0.000296, 0.000305, 0.000940] )
mystep = 5
if(mystep in thesteps):
casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
print 'Step ', mystep, step_title[mystep]
# PROBLEM 3 #
#########################################################################
# Plot visibilities of bandpass calibrator before applying instrumental calibration
# Plot amplitudes vs frequency
# Use long averaging time
# Find the right spectral window
# Do it for phase and amplitude
# Use some channel averaging
# Output to a png file
#########################################################################
####################################
# Apply WVR, Tsys and antenna position corrections
mystep = 6
if(mystep in thesteps):
casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
print 'Step ', mystep, step_title[mystep]
from recipes.almahelpers import tsysspwmap
tsysmap = tsysspwmap(vis = basename+'.ms', tsystable = basename+'.ms.tsys')
# PROBLEM 4 #
##########################################################
########## WHAT IS THE LOGIC OF tsysmap?
##########################################################
# PROBLEM 5 #
##########################################################
########## WHY DOES interp HAVE TWO VALUES, AND gainfield HAS THREE?
##########################################################
for field in ['J0522-364','J0648-3044','Vy CMa']: #bandpass, gain, target
applycal(vis = 'uid___A002_X6d34fd_Xb3d.ms',
spw='11',
field=field,
gainfield=[field,'',''],
interp=[PROBLEM,PROBLEM],
gaintable=[basename+'.ms.tsys',basename+'.ms.wvr.smooth',basename+'.ms.antpos'],
flagbackup=F, calwt=T,
spwmap = PROBLEM)
# PROBLEM 6 #
##########################################################
########## DETERMINE THE LINE FREE CHANNELS
########## (THERE IS MORE THAN ONE EMISSION LINE)
########## USE PLOTMS TO DO THIS
##########################################################
contchans='PROBLEM'
###############################################
# Split out science spw with time averaging
mystep = 7
if(mystep in thesteps):
casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
print 'Step ', mystep, step_title[mystep]
# PROBLEM 7 #
##########################################################
########## WRITE THE SPLIT COMMAND THAT SPLITS OUT THE CORRECTED COLUMN
##########################################################
os.system('rm -rf '+ basename+'.ms.split')
split(PROBLEM)
#########################################
# Listobs, clear pointing table, and save original flags
mystep=8
if(mystep in thesteps):
casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
print 'Step ', mystep, step_title[mystep]
os.system('rm -rf '+basename+'.ms.split.listobs.txt')
listobs(vis = basename+'.ms.split',
listfile = basename+'.ms.split.listobs.txt')
# Remove POINTING tables
tb.open(basename+'.ms.split/POINTING', nomodify = False)
a = tb.rownumbers()
tb.removerows(a)
tb.close()
# Store pre-flagging state
flagmanager(vis = basename+'.ms.split',
mode = 'save',
versionname = 'Pre-flagging')
mystep = 9
if(mystep in thesteps):
casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
print 'Step ', mystep, step_title[mystep]
# PROBLEM 8 #
#########################################################################
# Plot visibilities of bandpass calibrator AFTER applying instrumental calibration
# Plot amplitudes vs frequency
# Use long averaging time
# Find the right spectral window
# Do it for phase and amplitude
# Use some channel averaging
#########################################################################
print 'You should see that the scatter in the data is slightly less.'
#################
# Flagging shadowed data and any other bad data and back up flags
mystep = 10
if(mystep in thesteps):
casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
print 'Step ', mystep, step_title[mystep]
flagdata(vis = basename+'.ms.split',
mode = 'shadow',
flagbackup = F)
# DA41 is bad in the first data set
#flagdata(vis = 'uid___A002_X6d5bd2_X31.ms.split',antenna='DA41',
#mode = 'manual', scan='0~24')
flagmanager(vis = basename+'.ms.split',
mode = 'save',
versionname = 'BeforeBandpassCalibration')
#################################
# Bandpass calibration
# First perform time-dependent calibration, averaging the central channels
mystep = 11
if(mystep in thesteps):
casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
print 'Step ', mystep, step_title[mystep]
# PROBLEM 9 #
##########################################################
########## FIND A SUITABLE CHANNEL RANGE TO AVERAGE
########## FOR THE PRE-BANDPASS PHASE-UP
########## USE A SENSIBLE SOLUTION INTERVAL
##########################################################
os.system('rm -rf cal-uid___A002_X6d34fd_Xb3d-BPint.Gap')
gaincal(vis='uid___A002_X6d34fd_Xb3d.ms.split',
caltable='cal-uid___A002_X6d34fd_Xb3d-BPint.Gap',
spw='PROBLEM',
field='J0522-364',
selectdata=F, solint='PROBLEM',
refant=refantenna, calmode='ap')
if(calplots): # make some plots if calplots=True
os.system('rm -rf cal-uid___A002_X6d34fd_Xb3d-time-BPint-Gap.png')
plotcal(caltable='cal-uid___A002_X6d34fd_Xb3d-BPint.Gap',
xaxis = 'time', yaxis = 'phase',
plotsymbol='.:', plotrange = [0,0,-180,180],
figfile='cal-uid___A002_X6d34fd_Xb3d-time-BPint-Gap.png',
subplot = 211)
plotcal(caltable='cal-'+'uid___A002_X6d34fd_Xb3d-BPint.Gap',
xaxis = 'time', yaxis = 'amp',
plotsymbol='.:',
figfile='cal-uid___A002_X6d34fd_Xb3d-time-BPint-Gap.png',
subplot = 212)
# PROBLEM 10 #
##########################################################
########## FOR THE BANDPASS,
########## FIND A SUITABLE SOLUTION INTERVAL, IN TIME AND FREQUENCY
########## WHICH GAINTABLE TO APPLY ON THE FLY?
##########################################################
os.system('rm -rf cal-uid___A002_X6d34fd_Xb3d.B1')
bandpass(vis='uid___A002_X6d34fd_Xb3d.ms.split',
caltable='cal-uid___A002_X6d34fd_Xb3d.B1',
field='J0522-364',
bandtype='B',
fillgaps=10, solnorm = T,
solint='PROBLEM,PROBLEM',
refant=refantenna,
gaintable='PROBLEM')
if(calplots): # make some plots if calplots=True
os.system('rm -rf cal-uid___A002_X6d34fd_Xb3d-phase*-B1.png')
os.system('rm -rf cal-uid___A002_X6d34fd_Xb3d-amp*-B1.png')
for ant in range(3):
plotcal(caltable = 'cal-uid___A002_X6d34fd_Xb3d.B1',
xaxis='freq', yaxis='phase',
antenna=str(ant*8)+'~'+str(ant*8+7),
iteration='antenna',subplot=421,
plotrange = [0,0,-180,180],
figfile='cal-uid___A002_X6d34fd_Xb3d-phase_ant'+str(ant*8)+'-'+str(ant*8+7)+'-B1.png')
plotcal(caltable = 'cal-uid___A002_X6d34fd_Xb3d.B1',
xaxis='freq', yaxis='amp',
antenna=str(ant*8)+'~'+str(ant*8+7),
iteration='antenna',subplot=421,
plotrange = [0,0,0,1.5],
figfile='cal-uid___A002_X6d34fd_Xb3d-amp_ant'+str(ant*8)+'-'+str(ant*8+7)+'-B1.png')
# PROBLEM 11 #
##########################################################
########## FIND BANDPASS SOLUTIONS AGAIN, BUT NOW USE POLYNOMIAL
########## SOLUTIONS. WRITE THE CORRECT BANDPASS COMMAND
##########################################################
# PROBLEM 11A - OPTIONAL#
##########################################################
########## TRY TO OVERPLOT THE B AND BPOLY SOLUTIONS
#########################################################
#################################
# Setting the flux scale
mystep = 12
if(mystep in thesteps):
casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
print 'Step ', mystep, step_title[mystep]
# Remove any previous models for QSOs
delmod(vis=basename+'.ms.split', field='J*')
# There is no flux calibrator in this asdm, but we know the scaling
# from other observations. These are the numbers for J0522:
J0522_flux= 4.77004 # Flux density in Jy
J0522_rfreq=316.093 # The reference frequency for this flux, in GHz
J0522_spidx=-0.53027 # The spectral index
# PROBLEM 12 #
##########################################################
########## WRITE THE SETJY COMMAND THAT WRITES THE CORRECT FLUX VALUES
########## TO THE MODEL COLUMN
#########################################################
#################
# Time-dependent calibration
mystep = 13
if(mystep in thesteps):
casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
print 'Step ', mystep, step_title[mystep]
# PROBLEM 13 #
##########################################################
########## USE AN APPROPRIATE SHORT TIME SOLUTION INTERVAL FOR
########## THE 'FAST' PHASE SOLUTIONS
########## TRY LONGER SOLINT IF 'INT' IS TOO NOISY
########## WHAT CALMODE TO USE HERE?
########## WHAT IS AN APPROPRIATE MINSNR VALUE
#########################################################
# phase calibration per integration
os.system('rm -rf cal-'+basename+'.Gp_int')
gaincal(vis=basename+'.ms.split',
caltable = 'cal-'+basename+'.Gp_int',
field='J*',
minsnr=PROBLEM,
spw=contchans,
solint='PROBLEM',
refant=refantenna,
calmode='PROBLEM',
gaintable = 'cal-uid___A002_X6d34fd_Xb3d.B1')
if(calplots): # make some plots if calplots=True
os.system('rm -rf cal-'+basename+'-phase-time-Gp_int.png')
plotcal(caltable='cal-'+basename+'.Gp_int',
xaxis = 'time', yaxis = 'phase',poln='X',
plotsymbol='o',
figfile='cal-'+basename+'-phase-time-Gp_int.png',
subplot = 211)
plotcal(caltable='cal-'+basename+'.Gp_int',
xaxis = 'time', yaxis = 'phase',poln='Y',
plotsymbol='o',
figfile='cal-'+basename+'-phase-time-Gp_int.png',
subplot = 212)
# PROBLEM 14 #
##########################################################
########## USE AN APPROPRIATE TIME SOLUTION INTERVAL FOR
########## THE 'SLOW' PHASE SOLUTIONS
########## WHAT CALMODE TO USE HERE?
########## WHAT IS AN APPROPRIATE MINSNR VALUE
#########################################################
# phase calibration per scan to apply to target
os.system('rm -rf cal-'+basename+'.Gp_inf')
gaincal(vis=basename+'.ms.split',
caltable = 'cal-'+basename+'.Gp_inf',
field='J*',
minsnr=PROBLEM,
spw=contchans,
solint='PROBLEM',
refant=refantenna,
calmode='PROBLEM',
gaintable = 'cal-uid___A002_X6d34fd_Xb3d.B1')
if(calplots): # make some plots if calplots=True
os.system('rm -rf cal-'+basename+'-phase-time-Gp_inf.png')
plotcal(caltable='cal-'+basename+'.Gp_inf',
xaxis = 'time', yaxis = 'phase',poln='X',
plotsymbol='o:',
figfile='cal-'+basename+'-phase-time-Gp_inf.png',
subplot = 211)
plotcal(caltable='cal-'+basename+'.Gp_inf',
xaxis = 'time', yaxis = 'phase',poln='Y',
plotsymbol='o:',
figfile='cal-'+basename+'-phase-time-Gp_inf.png',
subplot = 212)
# PROBLEM 15 #
##########################################################
########## USE AN APPROPRIATE TIME SOLUTION INTERVAL FOR
########## THE 'SLOW' AMPLITUDE SOLUTIONS
########## WHAT CALMODE TO USE HERE?
########## WHAT IS AN APPROPRIATE MINSNR VALUE
########## WHICH CALIBRATION TABLE TO APPLY 'ON-THE-FLY'?
#########################################################
# amp cal per scan
os.system('rm -rf cal-'+basename+'.Gap_inf')
gaincal(vis=basename+'.ms.split',
caltable = 'cal-'+basename+'.Gap_inf',
field='J*',
minsnr=PROBLEM,
spw=contchans,
solint='PROBLEM',
refant=refantenna,
gaintype='T',
calmode='PROBLEM',
gaintable = ['cal-uid___A002_X6d34fd_Xb3d.B1','PROBLEM')
if(calplots): # make some plots if calplots=True
os.system('rm -rf cal-'+basename+'-time-Gap_inf.png')
plotcal(caltable='cal-'+basename+'.Gap_inf',
xaxis = 'time', yaxis = 'phase',
plotsymbol='o:',
figfile='cal-'+basename+'-time-Gap_inf.png',
subplot = 211)
plotcal(caltable='cal-'+basename+'.Gap_inf',
xaxis = 'time', yaxis = 'amp',
plotsymbol='o:',
figfile='cal-'+basename+'-time-Gap_inf.png',
subplot = 212)
# PROBLEM 16 #
##########################################################
########### TRANSFER THE FLUX SCALE TO THE GAIN CALIBRATOR
########### WRITE THE COMMAND
########### THE OUTPUT TABLE SHOULD BE CALLED cal-'+basename+'.Gap_flux_inf
########################################################
#################
# Apply calibration
mystep = 14
if(mystep in thesteps):
casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
print 'Step ', mystep, step_title[mystep]
applycal(vis='uid___A002_X6d34fd_Xb3d.ms.split',
field='J0522-364',
gaintable=['cal-uid___A002_X6d34fd_Xb3d.B1','cal-uid___A002_X6d34fd_Xb3d.Gp_inf','cal-uid___A002_X6d34fd_Xb3d.Gap_flux_inf'],
gainfield=['','J0522-364','J0522-364'],
interp='linear,linear',
calwt=F,flagbackup=F)
applycal(vis=basename+'.ms.split',
field='J0648-3044',
gaintable=['cal-uid___A002_X6d34fd_Xb3d.B1','cal-'+basename+'.Gp_inf','cal-'+basename+'.Gap_flux_inf'],
gainfield=['','J0648-3044','J0648-3044'],
interp='linear,linear',
calwt=F,flagbackup=F)
# Many phase ref solutions may have failed. If these would be applied to the
# target, the corresponding data would be flagged. As target is bright enough
# to self-cal, we do not want to flag target data with failed phase-ref solutions
#
# PROBLEM 17 #
#############################################################
####### WHICH APPLYMODE TO USE HERE
#############################################################
applycal(vis=basename+'.ms.split',
field='V*',
gaintable=['cal-uid___A002_X6d34fd_Xb3d.B1','cal-'+basename+'.Gp_inf','cal-'+basename+'.Gap_flux_inf'],
gainfield=['','J0648-3044','J0648-3044'],
interp='linear,linear',
applymode='PROBLEM',
calwt=F,flagbackup=F)
#######
# Plot phase-ref visibilities after applying instrumental calibration
mystep = 15
if(mystep in thesteps):
casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
print 'Step ', mystep, step_title[mystep]
os.system('rm -rf '+basename+'.ms.split.J0648-3044_*.png')
plotms(vis=basename+'.ms.split',
field='J0648-3044',
xaxis='time', yaxis='amp',
selectdata=T, spw='0',
avgtime='120s',avgchannel='3840',
ydatacolumn='corrected',
antenna=refantenna+'&*',
showgui=F,
plotfile=basename+'.ms.split.J0648-3044_amp.png')
plotms(vis=basename+'.ms.split',
field='J0648-3044',
xaxis='time', yaxis='phase',
selectdata=T, spw='0',
avgtime='120s',avgchannel='3840',
ydatacolumn='corrected',
antenna=refantenna+'&*',
showgui=F,
plotfile=basename+'.ms.split.J0648-3044_phase.png')
###########################
# Split out the science data, concatenate and transform to lsrk to allow self-calibration
mystep = 16
if(mystep in thesteps):
casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
print 'Step ', mystep, step_title[mystep]
os.system('rm -rf '+basename+'_VYCMa_658.ms')
mstransform(vis=basename+'.ms.split',
outputvis='VYCMa_658_origspw.ms',datacolumn='corrected',
field='V*')
# Regrid to LSRK
os.system('rm -rf VYCMa_658.ms')
mstransform(vis='VYCMa_658_origspw.ms',
outputvis='VYCMa_658.ms',datacolumn='data',
spw='0',
field='V*',combinespws=T,regridms=T,
mode='frequency',
outframe='lsrk')
# Do listobs and write info to a text file
os.system('rm VYCMa_658_listobs.txt')
listobs(vis='VYCMa_658.ms', listfile='VYCMa_658_listobs.txt', verbose=True)
print '<< Printed listobs output to VYCMa_658_listobs.txt'
#----------------------------------------------------------------------------------
#----- End of calibration script.
#----- To continue, see imaging script "VYCMa_Band9_Imaging.py"
#----------------------------------------------------------------------------------