#!/usr/bin/env ParselTongue
# e-MERLIN pipeline script
# Original version: mkargo 20110805.  Comments/complaints/cake to <argo@astron.nl>

################################################################################
# e-MERLIN pipeline. Invoke as:
#    mkargo@jaglan> parseltongue eMERLIN_pipeline.py calib.inputs
# where "calib.inputs" is a text file containing inputs for various tasks. See
# the readme file and the example in /home/mkargo/eMERLIN/
#
# This script will carry out a basic calibration including: initial fringe fit,
# normal phase referencing, and a bandpass calibration.  It can also optionally
# run selfcal on the phase calibrator.  It also runs IMAGR, but it almost
# certainly won't produce a nice pretty map you could publish.
#
# You need to have just the file you want to calibrate in your AIPS number
# (or you can load the DBCON'd file directly from disk by specifying it in
# the inputs file).
#
# You are STRONGLY advised to back up your tables outside of AIPS first.
################################################################################

# Changelog
# v0.9 20111202 - Adjustments to some functions for neatness!
# v0.8 20111123 - Fixed indisk errors with imagr and setjy indexing issue.
# v0.7 20111027 - Checks the specified reference antenna is present.  If not, sets
#	it to something sensible based on a priority list.
# v0.6 20111022 - Implemented loading a fitsfile from disk and some error checking
#	of inputs.
# v0.5 20111021 - Add some more diagnostic plots, including contour plots of the 
#	sources once imaged.
# v0.4 20111019 - Add option to do n rounds of selfcal on phase calibrator.
# v0.3 20111018 - Take point source fluxes from the inputs file instead of being 
#	hard-wired here.
# v0.2 20111018 - Added optional ability to set cellsize and imsize via the inputs 
#	file, also some sanity checking and uses sensible defaults if the input is crazy.
# v0.1 20111013 - First version of the working script.  Will calibrate and image a 
#	simple continuum dataset which has already been flagged, combined and averaged.





import os, re, time, datetime, sys, math
from os.path import join, getsize
from datetime import date
from collections import deque
import Utilities
from AIPS import AIPS, AIPSDisk
from AIPSTask import AIPSTask, AIPSList
from AIPSData import AIPSUVData, AIPSImage, AIPSCat
from eMERLIN_tasks import *
import eMERLIN_tasks



#################################################################################################################
# Define functions needed in this file
#################################################################################################################

# Parse the list of inputs given in the specified file. (Modified from Cormac's evn_funcs.py)
def parse_inp(filename):
    # form a hash of parameter names and values parsed from an input file.
    INPUTFILE = open(filename, "r")
    control = dict()

    # a few useful regular expressions
    newline = re.compile(r'\n')
    space = re.compile(r'\s')
    char = re.compile(r'\w')
    comment = re.compile(r'#.*')

    # parse the input file assuming '=' is used to separate names from values
    for line in INPUTFILE:
        if char.match(line):
            line = comment.sub(r'', line)
            line = line.replace("'", '')
            (param, value) = line.split('=')

            param = newline.sub(r'', param)
            param = param.strip()
            param = space.sub(r'', param)

            value = newline.sub(r'', value)
            value = value.strip()
            valuelist = value.split(', ')
            control[param] = valuelist
    return control


def checkin(control):
	# convert the control hash to a list of global variables
	global userno, indisk, fitsfile, toload, plotdir, msgkill
	global targets, phsrefs, fluxcals, pointcals, pointflux, bpasscals
	global refant, refantnames, band, imsize, cellsize
	global interactive, mode, solint, snver, fring_snr, setjy_fluxes
	global nchan, bchan, echan, nif, stokes, plotstokes, polarizations
	global dosort, doselfcal

	AIPS.userno = int(control.get('userno',[0])[0])
	indisk = int(control.get('indisk', [0])[0])
	dosort = int(control.get('dosort',[1])[0])
	doselfcal = int(control.get('doselfcal',[0])[0])

	refant = control.get('refant',['none'])[0]

	# Options relating to loading from disk
	fitsdir = control.get('fitsdir', [False])[0]
	fitsfile = control.get('fitsfile', [False])[0]
	toload = 0
	if (fitsfile and not fitsdir):
		fitsdir = os.getcwd()
	if (fitsdir and not fitsfile):
		print "Error: fitsdir specified but don't know which fitsfile to load. Check your inputs!"
		sys.exit()
	if (fitsdir and fitsfile):
		toload = os.path.join(fitsdir, fitsfile)
		if not os.path.isfile(toload):
			print "Error: file", toload, "does not exist. Check your inputs."
			sys.exit()

	# plotdir defaults to the cwd if not specified
	plotdir = control.get('plotdir', [False])[0]
	if (not plotdir):
		plotdir = os.getcwd()
	if (not os.path.isdir(plotdir)):
		print "Error:", plotdir, "does not exist. Check your inputs."
		sys.exit()

	solint = float(control.get('solint', [10])[0])

	bpasscals = control.get('bpasscals', [])
	for i in range(len(bpasscals)):
		bpasscals[i] = bpasscals[i].upper()
	nbp_table = 0
	if bpasscals:
		nbp_table = 1

	phsrefs = control.get('phsrefs', [])
	for i in range(len(phsrefs)):
		phsrefs[i] = phsrefs[i].upper()

	targets = control.get('targets', [])
	for i in range(len(targets)):
		targets[i] = targets[i].upper()

	pointcals = control.get('pointcals', [])
	for i in range(len(pointcals)):
		pointcals[i] = pointcals[i].upper()

	fluxcals = control.get('fluxcals', [])
	for i in range(len(fluxcals)):
		fluxcals[i] = fluxcals[i].upper()

	pointflux = control.get('pointflux', [])

	fring_snr = float(control.get('fring_snr', [9])[0])

	interactive = int(control.get('interactive', [False])[0])

	imsize = int(control.get('imsize', [1024])[0])
	if imsize:
		imsize = imsize,imsize

	cellsize = control.get('cellsize',[False])[0]

	msgkill = int(control.get('msgkill', [-5])[0])

	assert (len(phsrefs) == len(targets)), """unmatched number of phase reference and target sources! """

	if (not AIPS.userno) or (not indisk):
		print 'Error: You must set both your AIPS userno and indisk.'
		sys.exit()
	if (not targets) or (not phsrefs) or (not pointcals) or (not bpasscals):
		print 'Error: You must specify targets, phaserefs and bpasscals.'
		sys.exit()
	if not pointflux:
		print 'Error: You must specify the fluxes for your pointcal!'
		sys.exit()


def get_tab(uvdata, table):
	for i in range(len(uvdata.tables)) :
		if table in uvdata.tables[i][1] :
			ver = uvdata.tables[i][0]
	return ver


def get_ant_num(uvdata, refant_name):
	# convert antenna name to number
	antab = uvdata.table('AN',1)
	for row in antab :
		if refant_name in row.anname :
			return row.nosta


def set_refant(uvdata, refant):
	# Check the reference antenna is present, set to sensible default if not
	searchants = ['Mk2', 'Pi', 'Da', 'Kn', 'Cm', 'Lo', 'De']
	antab = uvdata.table('AN',1)
	refantn = 0
	refantn = get_ant_num(uvdata, refant)
	if not refantn :
		print "Warning: No refant specified or requested reference antenna not present."
		for item in searchants:
			if refantn:
				break
			for row in antab :
				if item in row.anname :
					print "Warning: Using reference antenna", item
					refant = item
					refantn = get_ant_num(uvdata, refant)
					break

	refantlist = []
	for antenna in searchants:
		antnum = get_ant_num(uvdata, antenna)
		if antnum:
			refantlist.append(antnum)

	return (refant, refantn, refantlist)


def setup(uvdata, dosort, cellsize):
	# extract useful information and set up a few useful parameters

	# Create a list of all cals and a list of all sources
	fringlist=[]
	for i in range(len(phsrefs)):
		fringlist.append(phsrefs[i])
	for i in range(len(fluxcals)):
		fringlist.append(fluxcals[i])
	for i in range(len(pointcals)):
		fringlist.append(pointcals[i])
	for i in range(len(bpasscals)):
		fringlist.append(bpasscals[i])
	# Remove duplicate source entries
	tmpfringlist = set(fringlist)
	fringlist=[]
	for source in tmpfringlist:
		fringlist.append(source)
	sourcelist=[]
	for i in range(len(fringlist)):
		sourcelist.append(fringlist[i])
	for i in range(len(targets)):
		sourcelist.append(targets[i])
	caliblist=[]
	for i in range(len(phsrefs)):
		caliblist.append(phsrefs[i])
	for i in range(len(pointcals)):
		caliblist.append(pointcals[i])

	flagver = 0
	for i in range(len(uvdata.tables)) :
		if 'FG' in uvdata.tables[i][1] :
			flagver = uvdata.tables[i][0]

	allsources = uvdata.sources

	# Set up uvdata and make sure it's in a fit state for use
	if dosort == 1 :
		print 'Tidying up your data before we start:'
		runmsort(uvdata)
		runindxr(uvdata)
	print 'Deleting old PL, SN, BP and CL tables.'
	uvdata.zap_table('SN', -1)
	uvdata.zap_table('PL', -1)
	uvdata.zap_table('BP', -1)
	for table in uvdata.tables :
		if 'CL' in table[1] :
			if table[0]>1 :
				uvdata.zap_table('CL', table[0])

	stokes = uvdata.stokes

	global plotstokes
	assert(len(stokes) > 0)
	if len(stokes) == 1:
		plotstokes = stokes[0]
	else:
		plotstokes = 'HALF'

	polarizations = uvdata.polarizations
	polarizations.sort()
	assert(len(polarizations) > 0)

	chanwidth = uvdata.header.cdelt[2]  # in Hz
	nchan = uvdata.header.naxis[2]
	# ignore the outer 10% of channels for certain tasks
	bchan = int(math.floor(nchan * 0.1))
	echan = nchan - bchan
	nif = uvdata.header.naxis[3]
	ifwidth = chanwidth*nchan

	freq = uvdata.header.crval[2] / 1000000000
	if (freq > 1.3) and (freq < 1.7):
		band = 'L'

	if (freq > 4) and (freq < 8):
		band = 'C'

	if (freq > 22) and (freq < 24):
		band = 'X'

	if (cellsize):
		cellsize = float(cellsize),float(cellsize)

	if (not cellsize):
		print "Warning: No image parameters specified, using default for this band."
		if (band == 'C'):
			cellsize = 0.01,0.01
		if (band == 'L'):
			cellsize = 0.14,0.14
		if (band == 'X'):
			cellsize = 0.01,0.01

	if ((band == 'C') and (cellsize[0] > 0.01)):
		print "Warning: Sub-optimal cellsize requested for C-band. Will image with cellsize = 0.01."
		cellsize = 0.01,0.01

	if ((band == 'L') and (cellsize[0] > 0.14)):
		print "Warning: Sub-optimal cellsize requested for L-band. Will image with cellsize = 0.14."
		cellsize = 0.14,0.14

	if ((band == 'X') and (cellsize[0] > 0.01)):
		print "Warning: Sub-optimal cellsize requested for X-band. Will image with cellsize = 0.01."
		cellsize = 0.01,0.01

	return (sourcelist, fringlist, caliblist, nchan, nif, stokes, plotstokes, polarizations, flagver, bchan, echan, freq, band, cellsize)


def dovplot(uvdata, docalib, gainuse, flagver, doband, bpver, step, mode):
	if (mode == 'a&p'):
		# Plot amp on cals (UVPLT, VPLOT)
		aparm = 0,0
		bparm = 12,1,0
		for pol in stokes:
			runvplot(uvdata, fringlist, pol, timer, refantlist, refantlist, 1, nif, bchan, echan, -1, 0, 0, -1, 0, aparm, bparm, refantn, dotv, 9)

		# LWPLA to disk
		filename = 'VPLOT_' + str(step) + '_amp.ps'
		outfile = os.path.join(plotdir, filename)
		print 'Plotting VPLOT diagnostics:' + outfile
		runlwpla(uvdata, outfile)

	# Plot phase on cals (UVPLT, VPLOT)
	aparm = 0,0
	bparm = 12,2,0
	runvplot(uvdata, fringlist, 'I', timer, refantlist, refantlist, 1, nif, bchan, echan, -1, 0, 0, -1, 0, aparm, bparm, refantn, dotv, 9)

	# LWPLA to disk
	filename = 'VPLOT_' + str(step) + '_phase.ps'
	outfile = os.path.join(plotdir, filename)
	print 'Plotting VPLOT diagnostics:' + outfile
	runlwpla(uvdata, outfile)





#################################################################################################################
# Start doing stuff
#################################################################################################################

# Sanity checks before we start

print "-----------------------------------------------------------"
print "********** Data reduction pipeline for e-MERLIN. **********"
print "-----------------------------------------------------------"
print "Your data should already be loaded, flagged and DBCON'd. If"
print "this is not the case, please quit now and come back here "
print "later."
print "-----------------------------------------------------------"

# Check for an inputs file
if len(sys.argv)==2 :
	print "Using inputs specified in", sys.argv[1]
	afile = sys.argv[1]
	if os.path.isfile(afile) :
		control = parse_inp(afile)
	else :
		print "Error:" + afile + "does not exist, quitting."
		sys.exit()
else :
	print "Error: no parameter file specified, quitting."
	print "Usage: parseltongue pipeline.py inputs.txt"
	sys.exit()

prompt = "Continue? ([y]/n): "
retries=1
complaint='Yes or no, please!'
while True :
	ok = raw_input(prompt)
	if ok == '' :
		print "You wanted to continue but were too lazy to say so."
		break
	if ok in ('y', 'ye', 'yes','Y'):
		print "Thank you, continuing."
		break
	if ok in ('n', 'no', 'nop', 'nope'):
		print "OK, quitting."
		sys.exit()
	retries = retries - 1
	if retries < 0:
		print "Don't be daft."
		sys.exit()
	print complaint



# Check the inputs
checkin(control)

if interactive:
	dotv = 1
else:
	dotv = -1
AIPSTask.msgkill = msgkill
logfile = os.path.join(plotdir, 'eMERLIN_pipeline.log')
AIPS.log = open(logfile, 'w')
print "Logging AIPS messages to", logfile


pca = AIPSCat(indisk)
if toload:
	if len(pca[indisk]) != 0 :
		print "You appear to already have file(s) on your AIPS disk, quitting."
		sys.exit()
	print "Loading ", toload
	runfitld(toload,indisk)
else:
	if len(pca[indisk])== 0 :
		print "You appear to have an empty AIPS disk. Did you forget to load your data?"
		sys.exit()
	if len(pca[indisk]) > 1 :
		print "You appear to have more than one file on your AIPS disk, quitting."
		sys.exit()

pca = AIPSCat(indisk)
print "Your AIPS catalog looks like this:"
print AIPSCat(indisk)


uvdata = AIPSUVData(pca[indisk][0]["name"], pca[indisk][0]["klass"], indisk, pca[indisk][0]["seq"])

(sourcelist, fringlist, caliblist, nchan, nif, stokes, plotstokes, polarizations, flagver, bchan, echan, freq, band, cellsize) = setup(uvdata, dosort, cellsize)

if nif != len(pointflux):
	print "Number of IFs does not match your point cal info, check \'pointflux\'!"
	sys.exit()


(refant, refantn, refantlist) = set_refant(uvdata, refant)


snver = 0
clver = 0


# Make some diagnostic plots before we start

timer = 0,0,0,0,0,0,0,0
aparm = 0,0,0,0,0,0,0,0,0,1
bparm = 0,0,0,0,0,0,0,0,0,0
for i in range(len(fringlist)) :
	sourceID = fringlist[i],
	print 'Running POSSM on', fringlist[i]
	runpossm(uvdata, sourceID, timer, refantlist, refantlist, aparm, bparm, 0, nif, -1, 0, flagver, plotstokes, -1, 0, 'A&P', 30, 9, -1, 1)

# LWPLA to disk
outfile = os.path.join(plotdir, 'POSSM_1.ps')
print 'Plotting POSSM diagnostics:' + outfile
runlwpla(uvdata, outfile)

# Plot VPLOT diagnostics
dovplot(uvdata, docalib=-1, gainuse=0, flagver=flagver, doband=-1, bpver=0, step=1, mode='a&p')


print '***********************************************************'
print '    Finished plotting pre-calibration diagnostic plots.'
print '         It would be wise to inspect them closely.'
print '***********************************************************'

# Initial fringe fit on everything *but* the target(s)
aparm = 3, 0, 0, 0, 0, 0, 2.5, 0, 0, 0
dparm = 3, 50, 0, 0, 0, 0, 0, 5, 0, 0 
timer = 0,0,0,0,0,0,0,0
solint = 10
print 'Running initial fringe fit (SN#1)'
runfring(uvdata, fringlist, timer, -1, 0, 0, -1, -1, refantn, refantlist, solint, aparm, dparm, 0, fring_snr, bchan, echan)

# Plot Delay solutions (SNPLT)
optype = 'DELA'
snver=0

snver = get_tab(uvdata, 'SN')
print 'SN version' + format(snver) + 'contains solutions for:', fringlist
runsnplt(uvdata, 'SN', 1, fringlist, 9, optype, dotv)

# LWPLA to disk
outfile = os.path.join(plotdir, 'SNPLT_1.ps')
print 'Plotting SNPLT diagnostics:' + outfile
runlwpla(uvdata, outfile)


# If in interactive mode: EDIT Delay Solutions (SNEDT) - remove extreme values
#	  ----> SN = SN + 1
#if interactive :
#	runsnedt()
#	snver=0
#	for i in range(len(uvdata.tables)) :
#		if 'SN' in uvdata.tables[i][1] :
#			snver = uvdata.tables[i][0]
#	runsnplt()
#	runlwpla()


# Apply solutions with CLCAL
runclcal(uvdata, sourcelist, fringlist, 'CALI', 'AMBG', snver, snver, 0, 0, refantn)



# N.B. Sanity checking that the number of fluxes given matches the number of IFs
# is done before processing starts.

fluxcal=[]
for i in range(len(pointcals)):
	print 'Running SETJY on', pointcals[i]
	for j in range(len(pointflux)):
		thing = float(pointflux[j])
		print "Setting flux for source " + str(pointcals[i]) + ", IF", str(j+1), '=', thing
		fluxcal.append(pointcals[i],)
		runsetjy(uvdata, fluxcal, j+1, j+1, thing, '')
		fluxcal.pop()


# Run calib - phase-only on phase and point source cals
aparm = 3, 0, 0, 0, 0, 0, 5, 0, 0, 0
cparm = 10, 0, 0, 0, 0, 0, 0
timer = 0,0,0,0,0,0,0,0
solint=10
uvrang = 0, 0
antwt = 0, 0
clver = get_tab(uvdata, 'CL')
print "Attempting phase-only calibration with CALIB", caliblist, "and clver", clver
runcalib(uvdata, caliblist, timer, uvrang, 100, clver, 0, -1, 0, 'DFT', refantn, solint, aparm, 1, 'L1', 'P', 10, 10, cparm, 0, antwt, 1)

# Run clcal
snver = get_tab(uvdata, 'SN')
print "Running clcal on snver", snver, "with", caliblist
runclcal(uvdata, '', caliblist, 'CALI', 'AMBG', snver, snver, 0, 0, refantn)





# Run calib - P&A on phase cal and point source cal
aparm = 3, 0, 0, 0, 0, 3, 5, 0, 0, 0
cparm = 10, 0, 0, 0, 0, 0, 0
timer = 0,0,0,0,0,0,0,0
solint=10
uvrang = 0, 0
antwt = 0, 0
clver = get_tab(uvdata, 'CL')
print "Running A&P calib on", caliblist, "with clver", clver
runcalib(uvdata, caliblist, timer, uvrang, 100, clver, 0, -1, 0, 'DFT', refantn, solint, aparm, 1, 'L1', 'A&P', 10, 10, cparm, 0, antwt, 1)


# If interactive, SNEDT
# If in interactive mode: EDIT Delay Solutions (SNEDT) - remove extreme values
#	  ----> SN = SN + 1
#if interactive :
#	runsnedt()
#	for i in range(len(uvdata.tables)) :
#		if 'SN' in uvdata.tables[i][1] :
#			snver = uvdata.tables[i][0]
#	runsnplt()
#	runlwpla()



# Run getjy
snver = get_tab(uvdata, 'SN')
#for phscal in phsrefs :
for i in range(len(phsrefs)) :
	sourceID = phsrefs[i],
	print "Running GETJY on phscal", sourceID[0]
	rungetjy(uvdata, sourceID, pointcals, 1, nif, snver)


# Run clcal
snver = get_tab(uvdata, 'SN')
print "Running clcal on snver", snver, "with", caliblist
runclcal(uvdata, '', caliblist, 'CALI', 'AMBG', snver, snver, 0, 0, refantn)




# More diagnostics:
# runpossm() using solint ~30mins on all cals.
timer = 0,0,0,0,0,0,0,0
aparm = 0,0,0,0,0,0,0,0,0,1
bparm = 0,0,0,0,0,0,0,0,0,0
for i in range(len(caliblist)) :
	sourceID = caliblist[i],
	print 'Running POSSM on', caliblist[i]
	runpossm(uvdata, sourceID, timer, refantlist, refantlist, aparm, bparm, 0, nif, 1, 0, 0, plotstokes, -1, 0, 'A&P', 30, 9, -1, 1)

# LWPLA to disk
outfile = os.path.join(plotdir, 'POSSM_2.ps')
print 'Plotting POSSM diagnostics:' + outfile
runlwpla(uvdata, outfile)


# Plot SNPLT diagnostics
clver = get_tab(uvdata, 'CL')
#dovplot(uvdata, docalib=1, gainuse=gainuse, flagver=0, doband=-1, bpver=0, step=2, mode='a&p')
runsnplt(uvdata, 'CL', clver, caliblist, nplots=9, optype='PHAS', dotv=dotv)
outfile = os.path.join(plotdir, 'SNPLT_2.ps')
print 'Plotting SNPLT diagnostics:' + outfile
runlwpla(uvdata, outfile)


# Bandpass solutions
soltyp = 'L1'
solint = -1
docal = 100
bpassprm = 0, 0, 0, 0, 1, 0, 0, 0, 1, 3, 1 
print 'Running BPASS.'
bpasssrc = bpasscals[0],
runbpass(uvdata, bpasssrc, refantn, bpassprm, 1, nif, soltyp, solint, docal)

timer = 0,0,0,0,0,0,0,0
aparm = 0,0,0,0,0,0,0,0,0,1
bparm = 0,0,0,0,0,0,0,0,0,0
for i in range(len(caliblist)) :
	sourceID = caliblist[i],
	print 'Running POSSM on', caliblist[i]
	runpossm(uvdata, sourceID, timer, refantlist, refantlist, aparm, bparm, 0, nif, 1, 0, 0, plotstokes, 1, 1, 'A&P', 30, 9, -1, 1)

# LWPLA to disk
outfile = os.path.join(plotdir, 'POSSM_3+BP.ps')
print 'Plotting POSSM diagnostics:' + outfile
runlwpla(uvdata, outfile)


step = 2


# Do selfcal if requested
# To cope with multiple sources we are going to have to create a list containing the name of each phase reference source and the corresponding CL table which is associated with it?
# For now, just assume we have only one phase cal and one source.
imver = 1
if doselfcal :
	mode = 'P'
	print "Self-calibration requested. Repeats:", doselfcal
	# Setup the base table for the start of each calibrator
	clverstart = get_tab(uvdata, 'CL')

	for i in range(doselfcal):
		j = i
		step = step + 1
		if i == (doselfcal - 1):
			mode = 'a&p'

		# image phase ref
		print "Mapping phase cal", phsrefs[0]
		calimsize = 256,256
		#for i in range(len(phsrefs)) :
		sourceID = phsrefs[0],
		runimagr(uvdata, sourceID, 1, 0, 0, 1, 0, bchan, echan, nchan, 1, cellsize, calimsize, 500, -1, indisk)


		aparm = 3, 0, 0, 0, 0, 0, 5, 0, 0, 0
		cparm = 10, 0, 0, 0, 0, 0, 0
		timer = 0,0,0,0,0,0,0,0
		solint=10
		uvrang = 0, 0
		antwt = 0, 0
		clver = get_tab(uvdata, 'CL')
		loop = i + 1
		uvdata2 = AIPSImage(sourceID[0],'ICL001',indisk,loop)
		print "Attempting self-calibration with", sourceID[0], "clver", clver, " and selfcal mode:", mode
		runselfcalib(uvdata, uvdata2, sourceID, timer, uvrang, 100, clver, 0, 1, 0, 'DFT', refantn, solint, aparm, 1, 'L1', mode, 10, 10, cparm, 0, antwt, 1)

		# Run clcal
		snver = get_tab(uvdata, 'SN')
		print "Running clcal on snver", snver, "with", sourceID[0]
		runclcal(uvdata, sourceID, sourceID, 'CALI', 'AMBG', snver, snver, 0, 0, refantn)

		# Plot VPLOT diagnostics
		gainuse = get_tab(uvdata, 'CL')
		runsnplt(uvdata, 'CL', gainuse, caliblist, nplots=9, optype='PHAS', dotv=dotv)
		outfile = os.path.join(plotdir, 'SNPLT_' + str(step) + '.ps')
		print 'Plotting SNPLT diagnostics:' + outfile
		runlwpla(uvdata, outfile)


	imver = j



print '***********************************************************'
print 'End of calibration procedures. Mapping.'
print '***********************************************************'


print "Mapping phase cal(s)", phsrefs
calimsize = 512,512
for i in range(len(phsrefs)) :
	sourceID = phsrefs[i],
	runimagr(uvdata, sourceID, 1, 0, 0, 1, 0, bchan, echan, nchan, 1, cellsize, calimsize, 500, -1, indisk)
	image = AIPSImage(sourceID[0],'ICL001',indisk,imver)
	runkntr(image,1)
	outfile = os.path.join(plotdir, sourceID[0] + str(imver) + '.ps')
	runlwpla(image,outfile)

print "Mapping target(s)", targets
for i in range(len(targets)) :
	sourceID = targets[i],
	runimagr(uvdata, sourceID, 1, 0, 0, 1, 0, bchan, echan, nchan, 1, cellsize, imsize, 500, -1, indisk)
	image = AIPSImage(sourceID[0],'ICL001',indisk,1)
	runkntr(image,5)
	outfile = os.path.join(plotdir, sourceID[0] + '.ps')
	runlwpla(image,outfile)



print '***********************************************************'
print "Finished. Now admire the awesomeness of your e-MERLIN maps."
print "Your diagnostic plots and pipelined images can be found in:"
print plotdir
print '***********************************************************'