# v190_pipeline.py
# script and notes for reducing v190 pulsar astrometry data in ParselTongue
# 2007 January 30 Adam Deller

# This script is intended to produce the basic calibration and flagging
# necessary for looking at averaged pulsar data in DIFMAP

################################################################################
# AIPS imports
################################################################################
from AIPS import AIPS, AIPSDisk
from AIPSTask import AIPSTask, AIPSList
from AIPSData import AIPSUVData, AIPSImage, AIPSCat
from AIPSTV import AIPSTV
#import Wizardry.AIPSData

################################################################################
# General python imports
################################################################################

import sys, os, string, math, warnings, subprocess
import v190_utils, logs2antab, wizardry_scintamp, wizardry_correctmodel, wizardry_selfcalsnfix
import wizardry_createdifmapflagfile, wizardry_correctpropermotion_oneshot, wizardry_sensweight
import wizardry_weightfix
from optparse import OptionParser
warnings.defaultaction = "always"

################################################################################
# Skip class
################################################################################
class SkipManager:
    skipload = skipmodelcorrect = skipuvsrt = skipflag = skipquack = skipfring = False
    skiptecor = skiptsys = skipsnedt = skipcal_t = skipcal_f = skipselfcal = False
    skipwritecals = skipscintweight = skipuvshift = skipwriterefs = False
    skipwritetargets = skippmcor = skipantabgen = skipmccheck = False
    skipplot = skipfixweights = False
    skipcreatedflags = skipusedflags = skipuvedt = skipsensweight = True
    
    def __init__(self):
        skipload = skipmodelcorrect = skipuvsrt = skipflag = skipquack = skipfring = False
        skiptecor = skiptsys = skipsnedt = skipcal_t = skipcal_f = skipselfcal = False
        skipwritecals = skipscintweight = skipuvshift = skipwriterefs = False
        skipwritetargets = skippmcor = skipantabgen = skipmccheck = False
        skipplot = skipfixweights = False
        skipcreatedflags = skipusedflags = skipuvedt = skipsensweight = True

################################################################################
# some global variables, some mine, some set from command line arguments
################################################################################
# Command line arguments
usage = "usage: %prog -e <experiment name> [options] <rpfits file 1> [rpfitsfile 2] ... [rpfits file n]"
tasks = ["load","uvsrt","modelcorrect","flag","quack","fring","tecor","tsys","snedt","cal_t",\
         "cal_f","uvedt","selfcal","writecals","scintweight","uvshift","writerefs",\
         "writetargets","sensweight","antabgen","mccheck","plot","fixweights","pmcor",\
         "skipcreatedflags","skipusedflags"]
avoidhelpstring = "Execution steps to skip (one or more of "
dohelpstring = "Execution steps to do, skipping all others (one or more of "
for task in tasks:
    avoidhelpstring = avoidhelpstring + task + ","
    dohelpstring = dohelpstring + task + ","
parser = OptionParser(usage)

parser.add_option("-e", "--experiment", dest="experiment", default="",
                  help="experiment name eg v190x", metavar="EXPERIMENT")
parser.add_option("-o", "--overwrite", dest="autostomp", action="store_true",default=False,
                  help="Overwrite existing AIPS catalogue files when required")
parser.add_option("-s", "--sortuv", dest="douvsrt", action="store_true",default=False,
                  help="Sort uv files after loading")
parser.add_option("-a", "--avoid", dest="avoid", default="",
                  help=avoidhelpstring,metavar="avoid1[,avoid2,...avoidN]")
parser.add_option("-c", "--clversion", dest="clversion", default=1,
                  help ="Starting cl version")
parser.add_option("-d", "--doonly", dest="doonly", default="",
                  help=dohelpstring,metavar='doonly1[,doonly2...doonlyN]')
parser.add_option("-t", "--notv", dest="notv", action="store_true", default=False,
                  help="Don't start the AIPS TV")
parser.add_option("-b", "--bandpasson", dest="doband", action="store_true", default=False,
                  help="Do bandpass correction - table must exist")

(options, rpfitsfilenames) = parser.parse_args()
experiment = options.experiment
douvsrt = options.douvsrt
autostomp = options.autostomp
doband = options.doband
print "Doband is " + str(doband)
clversion = int(options.clversion)
avoidkeys = options.avoid.split(',')
doonlykeys = options.doonly.split(',')
skipmanager = SkipManager()
if options.doonly != "":
    if options.avoid != "":
        print "Can only specify one of -a/--avoid= and -d/--doonly= - aborting!!!"
        sys.exit()
    for task in tasks:
        setattr(skipmanager, 'skip' + task, True)
    for task in doonlykeys:
        setattr(skipmanager, 'skip' + task, False)
elif options.avoid != "":
    for skipkey in avoidkeys:
        setattr(skipmanager, 'skip' + skipkey, True)
nfits = len(rpfitsfilenames)       # Number of RPFITS files to load
if nfits <= 0 and not skipmanager.skipload:
    parser.error("You must supply at least one RPFITS file!!!")
if experiment=="":
    parser.error("You must supply an experiment name with -e xxx or --experiment=xxx!!!")
  
# Other global variables
userno = 715
rootdir = os.path.dirname('/nfs/cluster/ska/adeller/v190/')
v190_utils.setexperiment(rootdir, experiment)
solint = float(180.0)          # main SOLINT in seconds !!!
t_int = float(2.0)            # correlator integration time in seconds !!!
selfcal_interval = float(1.5) # Solution interval for CALIB, in minutes
delay_window_ns = 400
rate_window_mhz = 100
cl_table_interval = 15.       # desired CL table interval in seconds !!!
flagprefix = experiment + '/logs/flag'
tsysprefix = experiment + '/logs/tsys'
modelsprefix = 'calmodels'
onlineflagsuffix = '.oflg'
skdflagsuffix = '.flag'
tsyssuffix = '.antab'
aips_version_std = '31DEC05'
calcdirectory = '/home/ssi/adeller/correlator/phase1fx/corcalc/'

# this stuff is used to find the instrumental phase in the absence of PCAL info
calibrator_instrument = '0537-441'
# put all FRING antennas here, with the desired refant first
calibrator_instrument_search = ['CAT*', 'PKS', 'MOPRA', 'TID', 'HOB', 'CED']
twochar_antenna_names = ['At','Pk','Mp','Ti','Ho','Cd']
# Source names, divided into fringe finders, phase refs and targets
fringefinders = ['0537-441','1921-293','0823-500','1718-649',
                 '1424-418','2255-282']
primary_phaserefs = ['0111-1317','0439-4522','0628-2805','0736-303',
                     '1604-4441','2047-1639','2141-3729','2142-0437']
target_pulsars = ['0108-1431','0437-4715','0630-2834','0737-3039',
                  '1559-4438','2048-1616','2144-3933','2145-0750']
exp_start_hours = [22,9,0,11,9,1,1,16,23,2,9,1]
#target_shifts = [[1.7,-1.4],[1.1,-0.7],[-0.75,1.05],[0,0],[-0.1,0.2],[2.7,-1.6],[-0.4,-1.6],[-0.1,0.4]]
pulsar_best_mar2006_shifts =  [[2.03,-1.81],[0.86,-0.5],[-0.055,0.72],[0.08,-0.11],[-0.1,0.2],[2.7,-1.55],[-0.4,-1.54],[-0.04,-0.07]]
pulsar_best_pm = [[103,-150],[121.55,-71.8],[-45,20],[0,0],[3,13],[105,-4],[-63,-153],[-7,28]]
#if experiment > 'v190h':
#    pulsar_best_mar2006_shifts[0][0] =
#    pulsar_best_mar2006_shifts[0][1] =
#    pulsar_best_mar2006_shifts[2][0] =
#    pulsar_best_mar2006_shifts[2][1] =
#    pulsar_best_mar2006_shifts[3][0] =
#    pulsar_best_mar2006_shifts[3][1] =
    
pulsar_scintillation_time = [120, 45, 5, 15, 1.5, 1.5, 60, 60]
calibrator_selfcaltimes = [1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5]
#bad_data_stop = [0, 0, 0, 0, 0, 0, 0, 0]

################################################################################
# Set up the reminaing global variable information
################################################################################

AIPS.userno = userno
#AIPSTask.isbatch = 0    # AIPS behaves differently in batch mode, so
                        # even though this is a "batch" mode script, don't
                        # tell that to AIPS.  If you really want batch mode,
                        # set this to 32000 (stupid AIPS magic numbers).

# declare a global uvdataset to work on
uvdata = AIPSUVData(experiment.upper(), 'UVDATA', 1, 1)
if uvdata.exists() and not skipmanager.skipload:
    if v190_utils.yesno(autostomp, 'Delete existing ' + experiment + ' data from catalogue and proceed (no aborts pipeline)?'):
        uvdata.zap()
    else:
        sys.exit()
if not options.notv:
    tv = AIPSTV()
    tv.start()

################################################################################
# open the log file
AIPS.log = open(rootdir + '/' + experiment + '.aipslog', "a")

################################################################################
# replace the warnings module's showwarning function with my own to print
# to the AIPS.log file
def showwarning(message, category, filename, lineno, file=None):
    """Hook to write a warning to AIPS.log file"""
    if file is None:
        file = AIPS.log
    try:
        file.write(warnings.formatwarning(message, category, filename, lineno))
        file.flush()
        if(file != sys.stderr):
            sys.stderr.write(warnings.formatwarning(message, category,
                                                    filename, lineno))
            sys.stderr.flush()
    except IOError:
        pass # the file (probably stderr) is invalid - this warning gets lost.

# Set the warnings function to use the AIPS log
warnings.showwarning = showwarning


################################################################################
################################################################################
################################################################################
################ SUBROUTINES ###################################################
################################################################################
################################################################################

################################################################################
# Grab the UT1-UTC and UTC-TAI values from the CALC files
def get_calc_utcvals(uvdata):
    obsdate = uvdata.header.date_obs
    year = string.atoi(obsdate[0:4])
    month = string.atoi(obsdate[5:7])
    day = string.atoi(obsdate[8:10])
    mjd = v190_utils.header_mjd(obsdate)
    utcvals = [0.0,0.0]

    eopin = open(calcdirectory + '/usno_finals.erp','r')
    eoplines = eopin.readlines()
    eopin.close()

    #find the line with our date
    for line in eoplines[1:]:
        splitline = line.split()
        #if len(splitline) > 1 and matchingmonth(splitline[0], month) and string.atoi(splitline[1]) == day:
        if len(splitline) > 1 and splitline[0] != '#' and splitline[0][2:7] == mjd:
            utcvals[0] = string.atof(splitline[3])/1000000.0

    taiin = open(calcdirectory + '/UTC-TAI.dat', 'r')
    tailines = taiin.readlines()[23:] #Ignore the header and some early years where things were tricky
    taiin.close()

    #Find the line with our year
    for line in tailines:
        startyear = string.atoi(line[1:5])
        startmonth = line[7:10].upper()
        if line[17:21] == '    ':
            endyear = 3000;
            endmonth = 'DEC'
        else:
            endyear = string.atoi(line[17:21])
            endmonth = line[23:26].upper()
        if dateincluded(year,month,startyear,getmonthint(startmonth),endyear,getmonthint(endmonth)):
            utcvals[1] = string.atof(line[34:36])
            break
    utcvals[0] = utcvals[0] + utcvals[1]

    return utcvals

################################################################################
# Gets the 1-based int for a month
def getmonthint(monthname):
    if monthname == 'JAN':
        return 1
    if monthname == 'FEB':
        return 2
    if monthname == 'MAR':
        return 3
    if monthname == 'APR':
        return 4
    if monthname == 'MAY':
        return 5
    if monthname == 'JUN':
        return 6
    if monthname == 'JUL':
        return 7
    if monthname == 'AUG':
        return 8
    if monthname == 'SEP':
        return 9
    if monthname == 'OCT':
        return 10
    if monthname == 'NOV':
        return 11
    if monthname == 'DEC':
        return 12
    return 0

################################################################################
# Checks to see if a number corresponds to a month
def matchingmonth(monthname, monthint):
    if monthname == 'JAN' and monthint == 1:
        return 1
    if monthname == 'FEB' and monthint == 2:
        return 1
    if monthname == 'MAR' and monthint == 3:
        return 1
    if monthname == 'APR' and monthint == 4:
        return 1
    if monthname == 'MAY' and monthint == 5:
        return 1
    if monthname == 'JUN' and monthint == 6:
        return 1
    if monthname == 'JUL' and monthint == 7:
        return 1
    if monthname == 'AUG' and monthint == 8:
        return 1
    if monthname == 'SEP' and monthint == 9:
        return 1
    if monthname == 'OCT' and monthint == 10:
        return 1
    if monthname == 'NOV' and monthint == 11:
        return 1
    if monthname == 'DEC' and monthint == 12:
        return 1
    return 0

################################################################################
# Checks to see if a year and month fall within a pair of year/month values
# assumes starting on the first of each month (appropriate for UTC-TAI file)
def dateincluded(year, month, startyear, startmonth, endyear, endmonth):
    if year < startyear: return 0
    if year > endyear: return 0
    if year == startyear and month < startmonth: return 0
    if year == endyear and month >= endmonth: return 0
    return 1

################################################################################
# Loads the RPFITS files into AIPS using ATLOD
def load_rpfits(nfiles, filenames, uvdata, doband, douvsrt):
    fillm = AIPSTask('fillm', version = aips_version_std)
    fillm.aparm[1:11] = [0,0,0,0,0,0,0,0,0,0]
    fillm.bparm[1:11] = [0,0,0,0,0,0,0,0,0,0]
    fillm.cparm[1:11] = [0,0,0,0,0,0,0,0,0,0]
    fillm.dparm[1:11] = [0,0,0,0,0,0,0,0,0,0]
    atlod = AIPSTask('atlod', version = aips_version_std)
    atlod.outdata = uvdata
    atlod.xyphase[1:][1:] = [0]
    atlod.optype = ''                       # Load data
    atlod.sources[1:] = ''                  # for all sources
    atlod.timerang[1:] = [0]                # for all times
    atlod.freqsel = 0                       # for all frequencies
    atlod.ifsel = 0                         # for all IFs
    atlod.chansel[1:][1:] = [0]             # for all channels
    atlod.ifmap = 0
    atlod.nif = 0
    atlod.douvcomp = -1
        
    for f in filenames:
        atlod.infile = v190_utils.check_root_file(experiment + '/' + f)
        for i in range(10):
            atlod.dparm[10-i] = 0.0
            atlod.cparm[10-i] = 0.0
            atlod.bparm[10-i] = 0.0
            atlod.aparm[10-i] = 0.0
        atlod.aparm[7] = cl_table_interval/60.0
        atlod.dparm[3] = 10.0
        print "Atlod.aparm[7] is " + str(atlod.aparm[7])
        if(atlod.infile):
            atlod()

    if doband:
        #Load the bandpass table
        tbin = AIPSTask('tbin',  version = '31DEC05')
        tbin.outdata = uvdata
        tbin.infile = rootdir + '/' + experiment + '/' + experiment + '.bpass.table'
        tbin.bcount = 1
        tbin.ecount = 0
        if not os.path.exists(tbin.infile):
            print "Can't find bandpass table file " + tbin.infile + " - turning off doband"
            doband = False
        else:
            tbin()
    if not douvsrt:
        uvdata.table('CL', 1).zap()
        indxr = AIPSTask('indxr', version = aips_version_std)
        indxr.indata = uvdata
        indxr.infile = ''
        indxr.prtlev = 1
        indxr.cparm[1:] = [0]
        indxr.cparm[3] = 0.25
        indxr.bparm[1:] = [0]
        indxr.in2file = ''
        indxr()
    
################################################################################
# Loads the flag file and attaches it to the AIPS uvdata
def load_flagfile(uvdata, dflagfile):
    uvflg = AIPSTask('uvflg', version = aips_version_std)
    uvflg.indata = uvdata
    uvflg.sources[1:] = ''
    uvflg.subarray = 0
    uvflg.selband = -1
    uvflg.selfreq = -1
    uvflg.freqid = -1
    uvflg.timer[1:] = [0,0,0,0,0,0,0,0]
    uvflg.bchan = 0
    uvflg.echan = 0
    uvflg.bif = 0
    uvflg.eif = 0
    uvflg.antennas[1:] = [0]
    uvflg.baseline[1:] = [0]
    uvflg.stokes = ''
    uvflg.flagver = 1
    uvflg.aparm[1:] = [0]
    uvflg.bdrop = 0
    uvflg.edrop = 0
    uvflg.dohist = -1
    uvflg.opcode = 'FLAG'
    if not dflagfile:
        uvflg.infile = v190_utils.check_data_file(flagprefix, '', skdflagsuffix)
        uvflg.reason = 'SLEWING'
        uvflg()
    else:
        s = ['RR','LL']
        uvflg.reason = 'DIFMAP EDIT'
        for i in range(4):
            for j in range(2):
                uvflg.bif = i
                uvflg.eif = i
                uvflg.stokes = s[j]
                try:
                    uvflg.infile = v190_utils.check_data_file(flagprefix, '', '.f' + str(i) + '.p' + str(j) + '.dflg')
                    uvflg()
                except AssertionError:
                    print 'No Difmap-generated flag file for IF ' + str(i) + ', polarisation ' + s[j]

################################################################################
# Sorts a loaded file into TB order 
def tb_sort():
    global uvdata
    
    uvsrt = AIPSTask('uvsrt', version = aips_version_std)
    uvsrt.indata = uvdata
    sorteddata = AIPSUVData(experiment.upper(), 'UVSRT', 1, 1)
    if sorteddata.exists():
        if v190_utils.yesno(autostomp, 'Delete existing sorted ' + experiment + ' data from catalogue and proceed (no aborts pipeline)?'):
            sorteddata.zap()
        else:
            sys.exit()
    uvsrt.outdata = sorteddata
    uvsrt.baddisk[1:] = [0]
    uvsrt.sort = 'TB'
    uvsrt.rotate = 0
    uvsrt.defer = 0
    uvsrt()
    indxr = AIPSTask('indxr', version = aips_version_std)
    indxr.indata = sorteddata
    indxr.infile = ''
    indxr.prtlev = 1
    indxr.cparm[1:] = [0]
    indxr.cparm[3] = 0.25
    indxr.bparm[1:] = [0]
    indxr.in2file = ''
    indxr()

################################################################################
# Loads the flag file and attaches it to the AIPS uvdata
def default_quack(uvdata, atcaname):
    quack = AIPSTask('quack', version = aips_version_std)
    quack.indata = uvdata
    quack.sources[1:] = ''
    quack.subarray = 0
    quack.selband = -1
    quack.selfreq = -1
    quack.freqid = -1
    quack.timerang[1:] = [0,0,0,0,0,0,0,0]
    quack.antennas[1:] = [0]
    quack.flagver = 1
    quack.opcode = 'BEG'
    quack.reason = 'QUACK'
    quack.aparm[1] = 0
    quack.aparm[2] = 0.167
    quack.aparm[3] = 0.5
    quack() #Do the first 10 seconds of each scan
    quack.opcode = 'ENDB'
    quack.aparm[2] = 0.083
    quack() #Do the last 5 seconds of each scan
    quack.antennas[1] = v190_utils.get_antenna_number(atcaname, uvdata)
    quack.opcode = 'BEG'
    quack.aparm[2] = 0.5
    quack() #Do the first 30 seconds at the ATCA
    
################################################################################
# Corrects for ionospheric dispersion using IONEX format Total Electron Content
def tecor_ionosphere(uvdata):
    tecor = AIPSTask('tecor', version = aips_version_std)
    tecor.indata = uvdata
    files = os.listdir(rootdir + '/' + experiment + '/logs/ionex/')
    numfiles = 0
    for filename in files:
        if filename[-1:] == 'i':
            if numfiles == 0 or (rootdir + '/' + experiment + '/logs/ionex/' + filename) < tecor.infile:
                tecor.infile = rootdir + '/' + experiment + '/logs/ionex/' + filename
            numfiles = numfiles + 1
    tecor.nfiles = numfiles
    tecor.subarr = 0
    tecor.antennas[1:] = [0]
    tecor.gainver = clversion
    tecor.gainuse = (clversion+1)
    tecor.aparm[2:] = [0]
    tecor.aparm[1] = 1
    tecor.aparm[2] = 0.2
    if numfiles > 0:
        print 'Running TECOR with ' + str(numfiles) + ' files - first one is ' + tecor.infile
        tecor() #Do it!
    else:
        print 'Could not find any suitable ionex files - aborting!!!'
        sys.exit()
        
################################################################################
# Loads the tsys data when already in sntable format
def load_tsys_sn(uvdata, path):
    tbin = AIPSTask('tbin', version = aips_version_std)
    tbin.outdata = uvdata
    tbin.infile = path
    tbin.bcount = 1
    tbin.ecount = 0
    tbin()
    #tacop = AIPSTask('tacop', version = aips_version_std)
    #tacop.indata = uvdata
    #tacop.inext = 'SN'
    #tacop.ncount = 1
    #tacop.invers = 1
    #tacop.outdata = uvdata
    #tacop.outver = 3
    #tacop.keyword = ''
    #tacop.keyvalue[1:] = [0]
    #tacop.keystrng = ''
    #tacop()

################################################################################
# Loads the tsys data from a single experiment.antab file and creates SN table 1
def load_tsys(uvdata):
    antab = AIPSTask('antab', version = aips_version_std)
    antab.indata = uvdata
    antab.infile = v190_utils.check_data_file(tsysprefix, '', tsyssuffix)
    antab.subarray = 0
    antab.tyver = 1
    antab.gcver = 1
    antab.blver = 1
    antab.sparm[1:] = ''
    antab.prtlev = 1
    antab.offset = 1
    antab()
    apcal = AIPSTask('apcal', version = aips_version_std)
    apcal.indata = uvdata
    apcal.antennas[1:] = [0]
    apcal.subarray = 0
    apcal.stokes = ''
    apcal.bif = 0
    apcal.eif = 0
    apcal.freqid = 0
    apcal.sources[1:] = ''
    apcal.timerang[1:] = [0,0,0,0,0,0,0,0]
    apcal.tyver = 1
    apcal.gcver = 1
    apcal.snver = 1
    apcal.opcode = ''
    apcal.aparm[1:] = [0]
    apcal()
    
################################################################################
# Launches interactive editing for the given SN table
def tvedit_sndata(uvdata, tablenumber):
    snedt = AIPSTask('snedt', version = aips_version_std)
    snedt.indata = uvdata
    snedt.inext = 'sn'
    snedt.invers = tablenumber
    snedt.dodelay = 1
    snedt.timerang[1:] = [0,0,0,0,0,0,0,0]
    snedt.bif = 0
    snedt.eif = 0
    snedt.antennas[1:] = [0]
    snedt.freqid = -1
    snedt.subarray = 0
    snedt.solint = 0.05
    snedt.detime = 0
    snedt.dotwo = 1
    snedt.expert = 0
    snedt.crowded = 0
    snedt.antuse[1:] = [0]
    snedt.baddisk[1:] = [0]
    snedt()
    
################################################################################
# Launches interactive editing for flagging UV data
def tvedit_uvdata(uvdata, cltableversion):
    editr = AIPSTask('editr', version = aips_version_std)
    editr.indata = uvdata
    editr.sources[1:] = ''
    editr.qual = -1
    editr.calcode = ''
    editr.timerang[1:] = [0,0,0,0,0,0,0,0]
    editr.selband = -1
    editr.selfreq = -1
    editr.bif = 0
    editr.eif = 0
    editr.bchan = 0
    editr.echan = 0
    editr.subarray = 0
    editr.docalib = 1
    editr.gainuse = cltableversion
    editr.dopol = -1
    editr.blver = -1
    editr.bpver = -1
    fqtable = uvdata.table('FQ', 1)[0]
    editr.smooth[3:] = [0]
    editr.smooth[1] = 3
    editr.smooth[2] = fqtable.total_bandwidth[0]/fqtable.ch_width[0]
    editr.stokes = ''
    editr.flagver = 1
    editr.antennas[1:] = [0]
    editr.freqid = -1
    editr.uvrange[1:] = [0]
    editr.dohist = -1
    editr.solint = 0.1
    editr.detime = 0
    editr.doweight = 1
    editr.dotwo = 1
    editr.expert = 0
    editr.crowded = 0
    editr.reason = 'EDITR'
    editr.antuse[1:] = [0]
    editr.baddisk[1:] = [0]
    editr()
    
################################################################################
# Fringe-fits one source only using that source's model, creating a SN table
def fringe_fit(uvdata, source, clversion, usemodel):
    if usemodel:
        #Load the clean component file using fitld
        modeldata = AIPSImage(source, 'UVDATA', 1, 1)
        if modeldata.exists():
            if v190_utils.yesno(autostomp, 'Delete existing ' + source + ' data from catalogue and proceed \
                                (no aborts pipeline)?'):
                modeldata.zap()
            else:
                sys.exit()
        fitld = AIPSTask('fitld', version = aips_version_std)
        fitld.infile = v190_utils.check_root_file(modelsprefix + '/' + source + '.model.fits')
        fitld.outdata = modeldata
        fitld.optype = ''
        fitld.ncount = 0
        fitld.dotable = 1
        fitld.douvcomp = 0
        fitld.doconcat = -1
        fitld()
    
    #Fringe fit into SN table 2
    print "Going to fringe-fit using clversion " + str(clversion)
    fring = AIPSTask('fring', version = aips_version_std)
    fring.indata = uvdata
    fring.calsour[1] = source
    fring.calsour[2:] = ''
    fring.qual = -1
    fring.calcode = ''
    fring.selband = 0
    fring.selfreq = 0
    fring.freqid = -1
    fring.timer[1:] = [0,0,0,0,0,0,0,0]
    fring.bchan = 0
    fring.echan = 0
    fring.antennas[1:] = [0]
    fring.dofit[1:] = [0]
    fring.subarray = 0
    fring.uvrange[1:] = [0]
    fring.wtuv = 0
    fring.weightit = 0
    fring.docalib = 2
    fring.gainuse = clversion
    fring.dopol = -1
    fring.blver = -1
    fring.flagver = 1
    fring.doband = -1
    fring.bpver = -1
    fring.smooth[1:] = [0]
    if usemodel:
        fring.in2data = modeldata
    else:
        fring.in2name = ''
        fring.in2class = ''
        fring.in2seq = 0
        fring.in2disk = 0
    fring.invers = 1
    fring.ncomp[1:] = [0]
    fring.flux = 0
    fring.nmaps = 0
    fring.cmethod = 'DFT'
    fring.smodel[1:] = [0]
    fring.refant = v190_utils.get_antenna_number(calibrator_instrument_search[0], uvdata)
    fring.search[1:] = [0]
    fring.aparm[1:] = [3,0,0,0,0,2,7.5,0,0,0]
    fring.dparm[1:] = [0,delay_window_ns,rate_window_mhz,t_int,0,0,0,0]
    if source == "0736-303":
        fring.solint = 9
        fring.solsub = 3
        fring.aparm[9] = 1
    else:
        fring.solint = solint/60.0
        fring.solsub = solint/120.0
    fring.solmin = 0
    fring.snver = 2
    fring.antwt[1:] = [0]
    fring.baddisk[1:] = [0]
    fring()
    
################################################################################
# Creates CL tables 3 and 4 from the Tsys and fringe-fit SN tables
def calibrate_tsys_fringefit(uvdata, skipcal_t, skipcal_f, skipplot):
    global clversion
    if not skipcal_t or not skipcal_f:
        clcal = AIPSTask('clcal', version = aips_version_std)
        clcal.indata = uvdata
        clcal.sources[1:] = ''
        clcal.soucode = ''
        clcal.calsour[1:] = ''
        clcal.qual = -1
        clcal.calcode = ''
        clcal.timerang[1:] = [0,0,0,0,0,0,0,0]
        clcal.subarray = 0
        clcal.antennas[1:] = [0]
        clcal.selband = -1
        clcal.selfreq = -1
        clcal.freqid = -1
        clcal.opcode = 'CALI'
        clcal.interpol = '2PT '
        clcal.cutoff = 0
        clcal.samptype = ''
        clcal.bparm[1:] = [0]
        clcal.doblank = 0
        clcal.dobtween = -1
        clcal.smotype = ''
    if not skipcal_t:
        try:
            table = uvdata.table('SN', 3)
            for row in table: break #Basically a test to see if it exists...
            clcal.snver = 3
            clcal.invers = 3
        except IOError:
            clcal.snver = 1
            clcal.invers = 1
        print "About to calibrate with snver and inver" + str(clcal.snver) + ", " + str(clcal.invers)
        clcal.gainver = clversion
        clcal.gainuse = (clversion+1)
        clcal.refant = v190_utils.get_antenna_number(calibrator_instrument_search[0], uvdata)
        clcal.baddisk[1:] = [0]
        for i in range(len(fringefinders)):
            if v190_utils.wasobserved(uvdata, fringefinders[i]):
                print "Calibrating tsys for " + fringefinders[i]
                clcal.sources[1] = fringefinders[i]
                clcal.calsour[1] = fringefinders[i]
                clcal() # Do the tsys values for this fringe finder
        for i in range(len(primary_phaserefs)):
            if v190_utils.wasobserved(uvdata, primary_phaserefs[i]):
                print "Calibrating tsys for " + primary_phaserefs[i] + " and " + target_pulsars[i]
                clcal.sources[1] = primary_phaserefs[i]
                clcal.sources[2] = target_pulsars[i]
                clcal.calsour[1] = primary_phaserefs[i]
                clcal.calsour[2] = target_pulsars[i]
                clcal() # Do the tsys values for this target/reference pair
        clversion = clversion+1
    if not skipcal_f:
        clcal.interpol = 'CUBE'
        clcal.sources[2] = ""
        clcal.calsour[2] = ""
        try:
            table = uvdata.table('SN', 4)
            for row in table: break #Basically a test to see if it exists...
            clcal.snver = 4
            clcal.invers = 4
        except IOError:
            clcal.snver = 2
            clcal.invers = 2
        clcal.gainver = clversion
        clcal.gainuse = (clversion+1)
        for i in range(len(fringefinders)):
            if v190_utils.wasobserved(uvdata, fringefinders[i]):
                print "Calibrating fring for " + fringefinders[i]
                clcal.sources[1] = fringefinders[i]
                clcal.calsour[1] = fringefinders[i]
                clcal() # Do the fringe-fit results for this fringe finder
        for i in range(len(primary_phaserefs)):
            if v190_utils.wasobserved(uvdata, primary_phaserefs[i]):
                print "Calibrating fring for " + primary_phaserefs[i] + " and " + target_pulsars[i]
                #if target_pulsars[i] == "0737-3039": #Want to smooth the delay solutions for this guy
                #    clcal.samptype = 'BOX'
                #    clcal.bparm[1] = clcal.bparm[2] = clcal.bparm[5] = solint/3600.0
                #    clcal.bparm[6] = clcal.bparm[7] = clcal.bparm[10] = solint/3600.0
                #    clcal.bparm[3] = clcal.bparm[4] = 10*solint/3600.0
                #    clcal.bparm[8] = clcal.bparm[9] = 10*solint/3600.0
                #    clcal.icut = 1.5
                #else:
                #    clcal.samptype = ''
                #    clcal.bparm[1:] = [0]
                #    clcal.icut = 0.0
                clcal.samptype = ''
                clcal.bparm[1:] = [0]
                clcal.icut = 0.0
                clcal.sources[1] = primary_phaserefs[i]
                clcal.sources[2] = target_pulsars[i]
                clcal.calsour[1] = primary_phaserefs[i]
                #clcal.calsour[2] = target_pulsars[i]
                clcal() # Do the fringe-fit results for this target/reference pair
        clversion = clversion+1
    if (not skipplot) and (not skipcal_t or not skipcal_f):
        snplt = AIPSTask('snplt', version = aips_version_std)
        snplt.indata = uvdata
        snplt.inext = 'CL'
        snplt.invers = clversion
        snplt.sources[1:] = ''
        snplt.qual = -1
        snplt.timerang[1:] = [0,0,0,0,0,0,0,0]
        snplt.stokes = ''
        snplt.selband = -1
        snplt.selfreq = -1
        snplt.freqid = -1
        snplt.subarray = 0
        snplt.antennas[1:] = [0]
        snplt.bif = 0
        snplt.eif = 0
        snplt.pixrange[1:] = [0, 0]
        snplt.nplots = 8
        snplt.xinc = 1
        snplt.opcode = ''
        snplt.do3col  = 0
        snplt.bcount = 0
        snplt.xaxis = 0
        snplt.symbol = 0
        snplt.factor = 0
        snplt.cutoff = 0
        snplt.ltype = 3
        snplt.dotv = 0
        snplt.grch = 1
        lwpla = AIPSTask('lwpla', version = aips_version_std)
        lwpla.userid = 0
        lwpla.indata = uvdata
        lwpla.plver = 1
        lwpla.invers = 1
        lwpla.aspmm = 0
        lwpla.lpen = 3
        lwpla.rgbgamma[1:] = [0]
        lwpla.functype = ''
        lwpla.dparm[1:] = [0]
        lwpla.copies = 1
        lwpla.dodark = 0
        lwpla.ofmfile = ' '
        lwpla.docolor = 0
        #lwpla.plcolors[1:] = [0]
    
        for output in ['AMP ','DELA','PHAS']:
            snplt.optype = output
            lwpla.outfile = rootdir + '/' + experiment + '/plots/' + experiment + \
                            '.tsys.' + output + '.ps'
            for ant in uvdata.antennas:
                snplt.antennas[1] = v190_utils.get_antenna_number(ant, uvdata)
                try:
                    snplt()
                except RuntimeError:
                    print "Looks like there was no data for antenna " + ant
                else:
                    lwpla()
                    uvdata.zap_table('AIPS PL', -1)

################################################################################
# Runs a selfcal on all the calibrators using CALIB, creating SN table 3
def aips_selfcal(uvdata, source):
    calib = AIPSTask('calib', version = aips_version_std)
    modeldata = AIPSImage(source, 'UVDATA', 1, 1)
    
    if not modeldata.exists():
        fitld = AIPSTask('fitld', version = aips_version_std)
        fitld.infile = v190_utils.check_root_file(modelsprefix + '/' + source + '.model.fits')
        fitld.outdata = modeldata
        fitld.optype = ''
        fitld.ncount = 0
        fitld.dotable = 1
        fitld.douvcomp = 1
        fitld.doconcat = -1
        fitld()
        
    calib.indata = uvdata
    calib.calsour[1] = source
    calib.calsour[2:] = ''
    calib.qual = -1
    calib.calcode = ''
    calib.selband = -1
    calib.selfreq = -1
    calib.freqid = -1
    calib.timerang[1:] = [0,0,0,0,0,0,0,0]
    calib.bchan = 0
    calib.echan = 0
    calib.antennas[1:] = [0]
    calib.dofit[1:] = [0]
    calib.antuse[1:] = [0]
    calib.subarray = 0
    calib.uvrange[1:] = [0, 0]
    calib.wtuv = 0
    calib.weightit = 0
    calib.docalib = 1
    calib.gainuse = clversion
    calib.dopol = -1
    calib.blver = -1
    calib.flagver = 1
    calib.doband = -1
    calib.smooth[1:] = [0]
    calib.in2data = modeldata
    calib.invers = 1
    calib.ncomp[1:] = [0]
    calib.flux = 0
    calib.nmaps = 0
    calib.cmethod = 'DFT'
    calib.cmodel = 'COMP'
    calib.smodel[1:] = [0]
    calib.refant = v190_utils.get_antenna_number(calibrator_instrument_search[0], uvdata)
    count = 0
    #Default value for selfcal interval
    calib.solint = selfcal_interval
    #Put in specific one
    for ref in primary_phaserefs:
        if ref == source:
            calib.solint = calibrator_selfcaltimes[count]
        count = count + 1
    calib.solsub = 0
    calib.solmin = 0
    calib.aparm[1:] = [4,0,0,0,0,3,8,0,0]
    calib.doflag = -2.5
    calib.soltype = ' '
    calib.solmode = 'A&P'
    calib.solcon = 0
    calib.minamper = 15
    calib.minphser = 10
    calib.cparm[1:] = [0,0,1,1,0,0]
    calib.snver = 5
    calib.antwt[1:] = [0]
    calib.baddisk[1:] = [0]
    try:
        calib()
    except RuntimeError:
        scount = 0
        print "Looks like there were no good solutions for source " + source
        for s in primary_phaserefs:
            if s == source:
                if scount != 0:
                    primary_phaserefs[:len(primary_phaserefs)-1] = primary_phaserefs[:scount] + primary_phaserefs[scount+1:]
                    target_pulsars[:len(target_pulsars)] = target_pulsars[:scount] + target_pulsars[scount+1:]
                    primary_phaserefs[len(primary_phaserefs)-1] = []
                    target_pulsars[len(target_pulsars)-1] = []
                else:
                    primary_phaserefs[0] = []
                    target_pulsars[0] = []
                #bad_data_stop[scount] = 1
            scount = scount + 1
    modeldata.zap()
    
################################################################################
# Runs a selfcal on the specified calibrators using difmap and cordump,loads the
# resultant SN table into AIPS and fiddle it with wizardry to duplicate necessary lines
# with an intelligent stab at a scaling factor.  Deletes the calibrator from disk
# once its finished
def difmap_selfcal(uvdata, isll, source, sourceindex, calno, basesntabno, startdayfrac):
    sourceuvfile = os.path.join(rootdir + "/" + experiment, experiment + "_" + source + "_uv.fits")
    modelfile = os.path.join(rootdir + "/" + modelsprefix, source + ".mod")
    
    #Set up the defaults that don't change, run to run
    difmap = subprocess.Popen("difmap", stdin=subprocess.PIPE)
    difmap.stdin.write("obs " + sourceuvfile + "\n")
    difmap.stdin.write("uvaver 10\n")
    if source == '0439-4522':
        difmap.stdin.write("mapsize 1024,0.4\n")
    else:
        difmap.stdin.write("mapsize 1024,2\n")
    difmap.stdin.write("rmod " + modelfile + "\n")
    difmap.stdin.write("device /xs\n")
    difmap.stdin.write("uvweight 2,-1\n")
    
    #First get the model sorted
    if isll:
        difmap.stdin.write("select ll,1,2,3,4\n")
    else:
        difmap.stdin.write("select i,1,2\n")
    if experiment == 'v190k':
        difmap.stdin.write("select rr,1,2,3,4\n")
    difmap.stdin.write("vplot\n")
    difmap.stdin.write("radplot\n")
    difmap.stdin.write("uvaver 60\n")
    print "First make a model you are happy with - enter to dump that model.  I suggest selfcal'ing at 1.5 min solution intervals"
    response = raw_input("Enter a difmap command - enter to go on to wmod")
    if response == 'loop':
        response = ''
        for i in range(20):
            response = response + 'modelfit 20\nselfcal true,true,1.5\n'
    while response != "":
        if response[:6] == 'uvaver':
            print "Averaging is not allowed!!!"
        else:
            difmap.stdin.write(response + "\n")
        response = raw_input("Enter a difmap command - enter to go on to wmod")
        if response == 'loop':
            response = ''
            for i in range(20):
                response = response + 'modelfit 20\nselfcal true,true,1.5\n'
    difmap.stdin.write("wmod " + rootdir + "/" + experiment + "/" + source + ".mod\n")
    difmap.stdin.write("save " + rootdir + "/" + experiment + "/junk" + source + "\n")

    #Then do left if required
    #difmap.stdin.write("obs " + sourceuvfile + "\n")
    #difmap.stdin.write("uvaver 20\n")
    #difmap.stdin.write("rmod " + rootdir + "/" + experiment + "/" + source + ".mod\n")
    if isll:
        difmap.stdin.write("select ll,1,2,3,4\n")
    else:
        difmap.stdin.write("select ll,1,2\n")
    if experiment == 'v190k':
        difmap.stdin.write("select rr,1,2,3,4\n")
    print "Selfcal (I recommend solution intervals of 1.5 minutes) until chi-squared is ok. NO MODEL-FITTING!!!"
    response = raw_input("Enter a difmap command - enter to go on to cordump: ")
    if response == 'loop':
        response = ''
        for i in range(20):
            response = response + 'selfcal true,true,1.5\n'
    while response != "":
        if response[:6] == 'uvaver':
            print "Averaging is not allowed!!!"
        else:
            difmap.stdin.write(response + "\n")
        response = raw_input("Enter a difmap command - enter to go on to cordump: ")
        if response == 'loop':
            response = ''
            for i in range(20):
                response = response + 'selfcal true,true,1.5\n'
    difmap.stdin.write("cordump " + rootdir + "/" + experiment + "/ll.sn\n")
    if not isll and experiment != 'v190k':
        #difmap.stdin.write("obs " + sourceuvfile + "\n")
        difmap.stdin.write("clrmod true\n")
        #difmap.stdin.write("uvaver 20\n")
        difmap.stdin.write("rmod " + rootdir + "/" + experiment + "/" + source + ".mod\n")
        difmap.stdin.write("select rr,1,2\n")
        print "Now selfcal for RR (I recommend solution intervals of 1.5 minutes) until chi-squared is ok. NO MODEL-FITTING!!!"
        response = raw_input("Enter a difmap command - enter to go on to cordump: ")
        if response == 'loop':
            response = ''
            for i in range(20):
                response = response + 'selfcal true,true,1.5\n'
        while response != "":
            if response[:6] == 'uvaver':
                print "Averaging is not allowed!!!"
            else:
                difmap.stdin.write(response + "\n")
            response = raw_input("Enter a difmap command - enter to go on to cordump: ")
            if response == 'loop':
                response = ''
                for i in range(20):
                    response = response + 'selfcal true,true,1.5\n'
        difmap.stdin.write("cordump " + rootdir + "/"  + experiment + "/rr.sn\n")
    difmap.stdin.write("exit\n\n")
    difmap.wait()

    #Change the sn file times if need be
    os.system("~/correlator/utilities/cordump_datefix " + rootdir + "/" + \
              experiment + "/ll.sn " + str(startdayfrac))
    if not isll and not experiment == 'v190k':
        os.system("~/correlator/utilities/cordump_datefix " + rootdir + "/" + \
                  experiment + "/rr.sn " + str(startdayfrac))
    
    if experiment == 'v190k':
        os.system("dasncon " + rootdir + "/" + experiment + "/ll.sn " + rootdir + "/" + experiment + \
                  "/ll.sn " + rootdir + "/" + experiment + "/all.sn")
    elif isll:
        os.system("dasncon " + rootdir + "/" + experiment + "/ll.sn " + rootdir + "/" + experiment + \
                  "/ll.sn " + rootdir + "/" + experiment + "/all.sn")
    else:
        os.system("dasncon " + rootdir + "/" + experiment + "/rr.sn " + rootdir + "/" + experiment + \
                  "/ll.sn " + rootdir + "/" + experiment + "/all.sn")
    snfile = os.path.join(rootdir + "/" + experiment, "all.sn")

    #Load the generated table into AIPS
    tbin = AIPSTask('tbin',  version = '31DEC05')
    tbin.outname = source
    tbin.outclass = 'SPLIT'
    tbin.outseq = 1
    tbin.outdisk = 1
    tbin.infile = snfile
    tbin.bcount = 1
    tbin.ecount = 0
    tbin()

    #Use wizardry to creat a new table
    #wizardry_selfcalsnfix.wizardry_selfcalcopy(experiment, uvdata.klass, source, isll, calno, basesntabno)
    #Not anymore, now just copy into position
    #copy_table(source, 'SPLIT', 1, 0, 'SN', 1, basesntabno+calno, uvdata)

    #Clean up
    os.system("rm -f " + sourceuvfile)
    os.system("rm -f " + experiment + "/*.sn")

################################################################################
# Load all the difmap files in prep for generating flag table from the difmap edits    
def load_junkfiles(uvdata, sources):
    modeldata = AIPSUVData(experiment.upper(), 'DFLAG', 1, -1)
    for i in range(99):
        modeldata = AIPSUVData(experiment.upper(), 'DFLAG', 1, i+1)
        if modeldata.exists():
            modeldata.zap()
    modeldata = AIPSUVData(experiment.upper(), 'DFLAG', 1, 0)
    fitld = AIPSTask('fitld', version = aips_version_std)
    fitld.outdata = modeldata
    fitld.optype = ''
    fitld.ncount = 0
    fitld.dotable = 1
    fitld.douvcomp = 1
    fitld.doconcat = 1
    fitld.qual = -1
    fitld.bchan = 0
    fitld.echan = 0
    fitld.bif = 0
    fitld.eif = 0
    fitld.digicor = -1
    fitld.wtthresh = 0

    nfound = 0
    for source in sources:
        if v190_utils.wasobserved(uvdata, source):
            try:
                fitld.infile = v190_utils.check_root_file(experiment + '/junk' + source + '.uvf')
                nfound = nfound + 1
                fitld()
            except AssertionError:
                print "No file " + rootdir + "/" + experiment + "/junk" + source + ".uvf; skipping"
    return nfound

################################################################################
# Produce clean maps, residual maps, and some statistics for a phase ref calibrator
def difmap_mapcal(source, isll):
    sourceuvfile = os.path.join(rootdir + "/" + experiment, experiment + "_" + source + "_uv.fits")
    if skipmanager.skipselfcal:
        modelfile = os.path.join(rootdir + "/" + modelsprefix, source + ".mod")
    else:
        modelfile = os.path.join(rootdir + "/" + experiment, source + ".mod")
    statsfile =  os.path.join(rootdir + "/" + experiment, source + ".stats")

    statsout = open(statsfile, "w")
    statsout.close()
    
    difmapout = open("difmap_commands", "w")
    difmapout.write("float pkflux\n")
    difmapout.write("float rmsflux\n")
    difmapout.write("integer ilevs\n")
    difmapout.write("float lowlev\n")
    difmapout.write("obs " + sourceuvfile + "\n")
    if source == '0439-4522':
        pixsize = 0.4
    else:
        pixsize = 2.5
    difmapout.write("mapsize 1024," + str(pixsize) + "\n")
    difmapout.write("mapcolor none\n")
    difmapout.write("uvweight 2,-1\n")
    difmapout.write("device /xs\n")
    difmapout.write("select ll,1,2,3,4\n")
    difmapout.write("vplot\n")
    difmapout.write("radplot\n")
    difmapout.write("select rr,1,2,3,4\n")
    difmapout.write("vplot\n")
    difmapout.write("radplot\n")
            
    #Do each individual band, one at a time
    for i in range(4):
        difmapout.write("select rr," + str(i+1) + "\n")
        difmapout.write("clrmod true\n")
        difmapout.write("rmod " + modelfile + "\n")
        difmapout.write("modelfit 20\n")
        difmapout.write("device " + rootdir + "/" + experiment + "/" + source + ".rr." + str(i+1) + ".resid.ps/PS\n")
        difmapout.write("mappl\n")
        difmapout.write("rmsflux = imstat(rms)\n")
        difmapout.write("restore\n")
        difmapout.write("pkflux = peak(flux,max)\n")
        difmapout.write("ilevs = pkflux/rmsflux\n")
        difmapout.write("lowlev = 300.0/ilevs\n")
        difmapout.write("loglevs lowlev\n")
        difmapout.write("device " + rootdir + "/" + experiment + "/" + source + ".rr." + str(i+1) + ".clean.ps/PS\n")
        difmapout.write("mappl cln\n")
        imagefile = rootdir + '/' + experiment + '/' + source + '.rr.' + str(i+1) + '.image.fits'
        if os.path.exists(imagefile):
            os.system("rm -f " + imagefile)
        difmapout.write("wmap " + imagefile + "\n")
        difmapout.write("print \"# Freq " + str(i+1) + " pol RR, S/N:\", ilevs\n")
        
        difmapout.write("select ll," + str(i+1) + "\n")
        difmapout.write("clrmod true\n")
        difmapout.write("rmod " + modelfile + "\n")
        difmapout.write("modelfit 20\n")
        difmapout.write("device " + rootdir + "/" + experiment + "/" + source + ".ll." + str(i+1) + ".resid.ps/PS\n")
        difmapout.write("mappl\n")
        difmapout.write("rmsflux = imstat(rms)\n")
        difmapout.write("restore\n")
        difmapout.write("pkflux = peak(flux,max)\n")
        difmapout.write("ilevs = pkflux/rmsflux\n")
        difmapout.write("lowlev = 300.0/ilevs\n")
        difmapout.write("loglevs lowlev\n")
        difmapout.write("device " + rootdir + "/" + experiment + "/" + source + ".ll." + str(i+1) + ".clean.ps/PS\n")
        difmapout.write("mappl cln\n")
        imagefile = rootdir + '/' + experiment + '/' + source + '.ll.' + str(i+1) + '.image.fits'
        if os.path.exists(imagefile):
            os.system("rm -f " + imagefile)
        difmapout.write("wmap " + imagefile + "\n")
        difmapout.write("print \"# Freq " + str(i+1) + " pol LL, S/N:\", ilevs,\"\\n\"\n")
        
    #Now do the supposed combined one
    if experiment == 'v190k':
        difmapout.write("select rr,1,2,3,4\n")
    elif isll:
        difmapout.write("select ll,1,2,3,4\n")
    else:
        difmapout.write("select i,1,2\n")
    difmapout.write("clrmod true\n")
    difmapout.write("rmod " + modelfile + "\n")
    difmapout.write("modelfit 20\n")
    difmapout.write("rmsflux = imstat(rms)\n")
    difmapout.write("device " + rootdir + "/" + experiment + "/" + source + ".combined.resid.ps/PS\n")
    difmapout.write("mappl\n")
    difmapout.write("restore\n")
    difmapout.write("pkflux = peak(flux,max)\n")
    difmapout.write("ilevs = pkflux/rmsflux\n")
    difmapout.write("lowlev = 300.0/ilevs\n")
    difmapout.write("loglevs lowlev\n")
    difmapout.write("device " + rootdir + "/" + experiment + "/" + source + ".combined.clean.ps/PS\n")
    difmapout.write("mappl cln\n")
    imagefile = rootdir + '/' + experiment + '/' + source + '.combined.image.fits'
    if os.path.exists(imagefile):
        os.system("rm -f " + imagefile)
    difmapout.write("wmap " + imagefile + "\n")
    difmapout.write("print \"# COMBINED S/N:\", ilevs\n")
    difmapout.write("exit\n")
    difmapout.close()

    os.system("difmap < difmap_commands >> " + statsfile)

################################################################################
# Produce clean maps, residual maps, and some statistics for a pulsar
def difmap_mappsr(source, isll, centre_ra_deg, centre_dec_deg):
    sourceuvfile = os.path.join(rootdir + "/" + experiment, experiment + "_" + source + "_uv.fits")
    psrinfofile = rootdir + "/" + experiment + "/" + source + ".psrinfo"

    if source == '0437-4715':
        pixsize = 0.4
    else:
        pixsize = 2.5

    psrinfoout = open(psrinfofile, "w")
    
    difmap = subprocess.Popen("difmap", stdin=subprocess.PIPE, \
                              stdout=psrinfoout)
    #difmapout = open("difmap_commands", "w")
    difmap.stdin.write("float pkflux\n")
    difmap.stdin.write("float peakx\n")
    difmap.stdin.write("float peaky\n")
    difmap.stdin.write("float finepeakx\n")
    difmap.stdin.write("float finepeaky\n")
    difmap.stdin.write("float rmsflux\n")
    difmap.stdin.write("integer ilevs\n")
    difmap.stdin.write("float lowlev\n")
    difmap.stdin.write("obs " + sourceuvfile + "\n")
    difmap.stdin.write("mapsize 1024," + str(pixsize) + "\n")
    difmap.stdin.write("uvweight 2,-1\n")
    #if not source == '2144-3933':
    #    difmap.stdin.write("uvaver 10,true\n")
    difmap.stdin.write("uvaver 20\n")
    difmap.stdin.write("mapcolor none\n")
    difmap.stdin.write("device /xs\n")
    difmap.stdin.write("select ll,1,2,3,4\n")
    response = raw_input("Enter a difmap command for the LL data - enter to go on to fitting")
    while response != "":
        difmap.stdin.write(response + "\n")
        response = raw_input("Enter a difmap command for the LL data - enter to go on to fitting")
    difmap.stdin.write("select rr,1,2,3,4\n")
    response = raw_input("Enter a difmap command for the RR data - enter to go on to fitting")
    while response != "":
        difmap.stdin.write(response + "\n")
        response = raw_input("Enter a difmap command for the RR data - enter to go on to fitting")
    difmap.stdin.write("save " + rootdir + "/" + experiment + "/junk" + source + "\n")

    #Find the peak from the combined first
    if experiment == 'v190k':
        difmap.stdin.write("select rr,1,2,3,4\n")
    elif isll:
        difmap.stdin.write("select ll,1,2,3,4\n")
    else:
        difmap.stdin.write("select i,1,2\n")
    difmap.stdin.write("peakx = peak(x,max)\n")
    difmap.stdin.write("peaky = peak(y,max)\n")
    difmap.stdin.write("shift -peakx,-peaky\n")
    difmap.stdin.write("mapsize 1024,0.1\n")
    difmap.stdin.write("finepeakx = peak(x,max)\n")
    difmap.stdin.write("finepeaky = peak(y,max)\n")
    
    #Do each individual band, one at a time
    pols = ['rr','ll']
    for i in range(4):
        for p in pols:
            difmap.stdin.write("select " + p + "," + str(i+1) + "\n")
            write_difmappsrscript(source, p + "." + str(i+1), difmap, pixsize)

    #Now do the supposed combined one
    if experiment == 'v190k':
        difmap.stdin.write("select rr,1,2,3,4\n")
    elif isll:
        difmap.stdin.write("select ll,1,2,3,4\n")
    else:
        difmap.stdin.write("select i,1,2\n")
    write_difmappsrscript(source, "combined", difmap, pixsize)
    difmap.stdin.write("exit\n\n")
    difmap.wait()
    psrinfoout.close()
    
    #difmapout.close()
    #os.system("difmap < difmap_commands >> " + psrinfofile)

    #Reprocess the psrinfo file into something more human readable
    write_psr_summary(source, centre_ra_deg, centre_dec_deg)

################################################################################
# Write the necessary difmap commands to a file for one band selection for one pulsar
def write_difmappsrscript(source, bands, difmap, pixsize):
    difmap.stdin.write("clrmod true\n")
    difmap.stdin.write("unshift\n")
    #difmap.stdin.write("peakx = peak(x,max)\n")
    #difmap.stdin.write("peaky = peak(y,max)\n")
    difmap.stdin.write("shift -peakx,-peaky\n")
    difmap.stdin.write("mapsize 1024,0.1\n")
    #difmap.stdin.write("finepeakx = peak(x,max)\n")
    #difmap.stdin.write("finepeaky = peak(y,max)\n")
    difmap.stdin.write("pkflux = peak(flux,max)\n")
    difmap.stdin.write("addcmp pkflux, true, finepeakx, finepeaky, true, 0, false, 1, false, 0, false, 0, 0, 0\n")
    difmap.stdin.write("mapsize 1024," + str(pixsize) + "\n")
    difmap.stdin.write("modelfit 50\n")
    difmap.stdin.write("rmsflux = imstat(rms)\n")
    difmap.stdin.write("device " + rootdir + "/" + experiment + "/" + source + "." + bands + ".resid.ps/PS\n")
    difmap.stdin.write("mappl\n")
    difmap.stdin.write("restore\n")
    difmap.stdin.write("pkflux = peak(flux,max)\n")
    difmap.stdin.write("ilevs = pkflux/rmsflux\n")
    difmap.stdin.write("lowlev = 300.0/ilevs\n")
    difmap.stdin.write("loglevs lowlev\n")
    difmap.stdin.write("device " + rootdir + "/" + experiment + "/" + source + "." + bands + ".clean.ps/PS\n")
    difmap.stdin.write("mappl cln\n")
    imagefile = rootdir + '/' + experiment + '/' + source + '.' + bands + '.image.fits'
    if os.path.exists(imagefile):
        os.system("rm -f " + imagefile)
    difmap.stdin.write("wmap " + imagefile + "\n")
    difmap.stdin.write("unshift\n")
    difmap.stdin.write("wmod\n")
    difmap.stdin.write("print rmsflux\n")
    difmap.stdin.write("print imstat(bmin)\n")
    difmap.stdin.write("print imstat(bmaj)\n")
    difmap.stdin.write("print imstat(bpa)\n")
    
################################################################################
# Translates the dumped pulsar data from difmap into human readable form, deletes dump data
def write_psr_summary(source, centre_ra_deg, centre_dec_deg):
    psrinfoin = open(rootdir + "/" + experiment + "/" + source + ".psrinfo", "r")
    psrinfolines = psrinfoin.readlines()
    psrinfoin.close()
    psrstatout = open(rootdir + "/" + experiment + "/" + source + ".psrstats", "w")

    #Get the rest of the data and translate it
    at = 0
    pols = ['rr','ll']
    for i in range(4):
        for p in pols:
            rahh = int(centre_ra_deg/15.0)
            ramm = int((centre_ra_deg/15.0 - float(rahh))*60.0)
            rass = (centre_ra_deg/15.0 - (float(rahh) + float(ramm)/60.0))*3600.0
            decdd = int(centre_dec_deg)
            decmm = int((centre_dec_deg - float(decdd))*60.0)
            decss = (centre_dec_deg - (decdd + float(decmm)/60.0))*3600.0
            psrstatout.write("Pulsar " + source + ": Frequency " + str(i) + ", pol " + p.upper() + ":\n")
            while psrinfolines[at][0:44] != "Writing 1 model components to file: (stdout)":
                at = at + 1
            at = at + 1
            while psrinfolines[at][0] == '!':
                at = at + 1
            splitline = psrinfolines[at].split()
            flux = float(splitline[0][:-1])
            radius = float(splitline[1][:-1])
            thetadeg = float(splitline[2][:-1])
            rms = float(psrinfolines[at+1])
            bmin = float(psrinfolines[at+3])
            bmaj = float(psrinfolines[at+4])
            bpadeg =  float(psrinfolines[at+5])
            theta = thetadeg*math.pi/180.0
            bpa = bpadeg*math.pi/180.0
            at = at+6
            psrstatout.write("Centre RA:            %02d:%02d:%09.6f\n" % (rahh, ramm, rass))
            psrstatout.write("Centre Dec:           %02d:%02d:%09.6f\n" % (decdd, abs(decmm), abs(decss)))
            psrstatout.write("Flux (mJy):           " + str(flux*1000.0) + "\n")
            psrstatout.write("S/N:                  " + str(flux/rms) + "\n")
            psrstatout.write("RA offset (mas):      " + str(radius*math.sin(theta)) + "\n")
            psrstatout.write("Dec offset (mas):     " + str(radius*math.cos(theta)) + "\n")
            ra = centre_ra_deg + radius*math.sin(theta)/(math.cos(centre_dec_deg*math.pi/180.0)*3600000.0)
            dec = centre_dec_deg + radius*math.cos(theta)/3600000.0
            rahh = int(ra/15.0)
            ramm = int((ra/15.0 - float(rahh))*60.0)
            rass = (ra/15.0 - (float(rahh) + float(ramm)/60.0))*3600.0
            decdd = int(dec)
            decmm = int((dec - float(decdd))*60.0)
            decss = (dec - (decdd + float(decmm)/60.0))*3600.0
            psrstatout.write("Actual RA:            %02d:%02d:%09.6f\n" % (rahh, ramm, rass))
            psrstatout.write("Actual Dec:           %02d:%02d:%09.6f\n" % (decdd, abs(decmm), abs(decss)))
            psrstatout.write("Beam:                 " + str(bmin) + "x" + str(bmaj) + " at " + str(bpadeg) + " degrees\n")
            raerr = math.sqrt((bmaj*bmaj*math.sin(bpa)*math.sin(bpa) + bmin*bmin*math.cos(bpa)*math.cos(bpa))/
                              (1+bmin*bmin/(bmaj*bmaj)))/(flux/rms)
            decerr = math.sqrt((bmaj*bmaj*math.cos(bpa)*math.cos(bpa) + bmin*bmin*math.sin(bpa)*math.sin(bpa))/
                               (1+bmin*bmin/(bmaj*bmaj)))/(flux/rms)
            psrstatout.write("Est. RA error (mas):  " + str(raerr) + "\n")
            psrstatout.write("Est. RA error (hms):  " + str(raerr/(math.cos(centre_dec_deg*math.pi/180.0)*15.0)) + "\n")
            psrstatout.write("Est. Dec error (mas): " + str(decerr) + "\n")
                               
    psrstatout.write("Pulsar " + source + " combined results:\n")

    while psrinfolines[at][0:44] != "Writing 1 model components to file: (stdout)":
        at = at + 1
    at = at + 1
    while psrinfolines[at][0] == '!':
        at = at + 1
       
    rahh = int(centre_ra_deg/15.0)
    ramm = int((centre_ra_deg/15.0 - float(rahh))*60.0)
    rass = (centre_ra_deg/15.0 - (float(rahh) + float(ramm)/60.0))*3600.0
    decdd = int(centre_dec_deg)
    decmm = int((centre_dec_deg - float(decdd))*60.0)
    decss = (centre_dec_deg - (decdd + float(decmm)/60.0))*3600.0 
    splitline = psrinfolines[at].split()
    flux = float(splitline[0][:-1])
    radius = float(splitline[1][:-1])
    thetadeg = float(splitline[2][:-1])
    rms = float(psrinfolines[at+1])
    bmin = float(psrinfolines[at+3])
    bmaj = float(psrinfolines[at+4])
    bpadeg =  float(psrinfolines[at+5])
    theta = thetadeg*math.pi/180.0
    bpa = bpadeg*math.pi/180.0
    at = at+6
    psrstatout.write("Centre RA:            %02d:%02d:%09.6f\n" % (rahh, ramm, rass))
    psrstatout.write("Centre Dec:           %02d:%02d:%09.6f\n" % (decdd, abs(decmm), abs(decss)))
    psrstatout.write("Flux (mJy):           " + str(flux*1000.0) + "\n")
    psrstatout.write("S/N:                  " + str(flux/rms) + "\n")
    psrstatout.write("RA offset (mas):      " + str(radius*math.sin(theta)) + "\n")
    psrstatout.write("Dec offset (mas):     " + str(radius*math.cos(theta)) + "\n")
    ra = centre_ra_deg + radius*math.sin(theta)/(math.cos(centre_dec_deg*math.pi/180.0)*3600000.0)
    dec = centre_dec_deg + radius*math.cos(theta)/3600000.0
    rahh = int(ra/15.0)
    ramm = int((ra/15.0 - float(rahh))*60.0)
    rass = (ra/15.0 - (float(rahh) + float(ramm)/60.0))*3600.0
    decdd = int(dec)
    decmm = int((dec - float(decdd))*60.0)
    decss = (dec - (decdd + float(decmm)/60.0))*3600.0
    psrstatout.write("Actual RA:            %02d:%02d:%09.6f\n" % (rahh, ramm, rass))
    psrstatout.write("Actual Dec:           %02d:%02d:%09.6f\n" % (decdd, abs(decmm), abs(decss)))
    psrstatout.write("Beam:                 " + str(bmin) + "x" + str(bmaj) + " at " + str(bpadeg) + " degrees\n")
    raerr = math.sqrt((bmaj*bmaj*math.sin(bpa)*math.sin(bpa) + bmin*bmin*math.cos(bpa)*math.cos(bpa))/
                      (1+bmin*bmin/(bmaj*bmaj)))/(flux/rms)
    decerr = math.sqrt((bmaj*bmaj*math.cos(bpa)*math.cos(bpa) + bmin*bmin*math.sin(bpa)*math.sin(bpa))/
                       (1+bmin*bmin/(bmaj*bmaj)))/(flux/rms)
    psrstatout.write("Est. RA error (mas):  " + str(raerr) + "\n")
    psrstatout.write("Est. RA error (hms):  " + str(raerr/(math.cos(centre_dec_deg*math.pi/180.0)*15.0)) + "\n")
    psrstatout.write("Est. Dec error (mas): " + str(decerr) + "\n")
    psrstatout.close()
    
################################################################################
# Merges all the separate SN tables generated by selfcal'ing on each reference source,
# renumbers the merged table and then deletes all the intermediate tables
def mergesn(uvdata, phaserefs, numcals, basesntabno):
    clcal = AIPSTask('clcal', version = aips_version_std)
    clcal.indata = uvdata
    clcal.sources[1:] = ''
    clcal.soucode = ''
    clcal.calsour[1:] = ''
    clcal.qual = -1
    clcal.calcode = ''
    clcal.timerang[1:] = [0]
    clcal.subarray = 0
    clcal.antennas[1:] = [0]
    clcal.selband = -1
    clcal.selfreq = -1
    clcal.freqid = -1
    clcal.opcode = 'MERG'
    clcal.samptype = 'BOX'
    clcal.interpol = '2PT'
    clcal.cutoff = 0
    clcal.bparm[1:] = [0]
    clcal.bparm[1] = 1
    clcal.doblank = 0
    clcal.dobtween = -1
    clcal.smotype = 'AMPL'
    clcal.snver = basesntabno+1
    clcal.inver = basesntabno+numcals
    clcal.refant = v190_utils.get_antenna_number(calibrator_instrument_search[0], uvdata)
    clcal.baddisk[1:] = [0]
    clcal() # Merge the tables from endsn+1 to endsn+numcals+1 into endsn+numcals+2

    #tacop = AIPSTask('tacop', version = aips_version_std)
    #tacop.indata = uvdata
    #tacop.inext = 'SN'
    #tacop.ncount = 1
    #tacop.invers = 0
    #tacop.outdata = uvdata
    #tacop.outvers = basesntabno
    #tacop.keyword = ''
    #tacop.keyvalue[1:] = [0]
    #tacop.keystrng = ''
    #tacop() #Copy the generated table to a sensible version number

    #for i in range(numcals+1):
    #    uvdata.table('SN', basesntabno+i+1).zap() #Zap all the intermediate files

################################################################################
# Creates CL table from the calibrator selfcal corrections
def calibrate_selfcalcorrections(uvdata, sntableno, skipplot, sourceindex):
    global clversion
    clcal = AIPSTask('clcal', version = aips_version_std)
    clcal.indata = uvdata
    clcal.sources[1:] = ''
    clcal.soucode = ''
    clcal.calsour[1:] = ''
    clcal.qual = -1
    clcal.calcode = ''
    clcal.timerang[1:] = [0]
    clcal.subarray = 0
    clcal.antennas[1:] = [0]
    clcal.selband = -1
    clcal.selfreq = -1
    clcal.freqid = -1
    clcal.opcode = 'CALI'
    clcal.interpol = '2PT'
    clcal.cutoff = 0
    clcal.samptype = ''
    clcal.bparm[1:] = [0]
    clcal.doblank = 1
    clcal.dobtween = -1
    clcal.smotype = ''
    clcal.snver = sntableno
    clcal.inver = sntableno
    clcal.gainver = clversion
    clcal.gainuse = (clversion+1)
    clcal.refant = v190_utils.get_antenna_number(calibrator_instrument_search[0], uvdata)
    clcal.baddisk[1:] = [0]
    if sourceindex < 0:
        for i in range(len(primary_phaserefs)):
            if v190_utils.wasobserved(uvdata, primary_phaserefs[i]):
                clcal.sources[1] = primary_phaserefs[i]
                clcal.sources[2] = target_pulsars[i]
                clcal.calsour[1] = primary_phaserefs[i]
                clcal() # Do the selfcal corrections
        clversion = clversion+1
    else:
        clcal.sources[1] = primary_phaserefs[sourceindex]
        clcal.sources[2] = target_pulsars[sourceindex]
        clcal.calsour[1] = primary_phaserefs[sourceindex]
        clcal() # Do the selfcal corrections

    if not skipplot and sourceindex == 0:
        snplt = AIPSTask('snplt', version = aips_version_std)
        snplt.indata = uvdata
        snplt.inext = 'CL'
        snplt.invers = clversion
        snplt.sources[1:] = ''
        snplt.qual = -1
        snplt.timerang[1:] = [0]
        snplt.stokes = ''
        snplt.selband = -1
        snplt.selfreq = -1
        snplt.freqid = -1
        snplt.subarray = 0
        snplt.bif = 0
        snplt.eif = 0
        snplt.pixrange[1:] = [0, 0]
        snplt.nplots = 8
        snplt.xinc = 1
        snplt.opcode = ''
        snplt.do3col  = 0
        snplt.bcount = 0
        snplt.xaxis = 0
        snplt.symbol = 0
        snplt.factor = 0
        snplt.cutoff = 0
        snplt.ltype = 3
        snplt.dotv = 0
        lwpla = AIPSTask('lwpla', version = aips_version_std)
        lwpla.userid = 0
        lwpla.indata = uvdata
        lwpla.plver = 1
        lwpla.invers = 1
        lwpla.aspmm = 0
        lwpla.lpen = 3
        lwpla.rgbgamma[1:] = [0]
        lwpla.functype = ''
        lwpla.dparm[1:] = [0]
        lwpla.copies = 1
        lwpla.dodark = 0
        lwpla.ofmfile = ''
        lwpla.docolor = 0
#        lwpla.plcolors[1:] = [0]  
    
        for output in ['AMP','DELA','PHAS']:
            snplt.optype = output
            lwpla.outfile = rootdir + '/' + experiment + '/plots/' + \
                            experiment + '.selfcal.' + output + '.ps'
            for ant in uvdata.antennas:
                snplt.antennas[1] =  v190_utils.get_antenna_number(ant, uvdata)
                try:
                    snplt()
                except RuntimeError:
                    print "Looks like there was no data for antenna " + ant
                else:
                    lwpla()
                    uvdata.zap_table('AIPS PL', -1)

################################################################################
# Splits a list of sources, optionally with averaging
def split(uvdata, sourcelist, doaverage, clversion, doband):
    split = AIPSTask('split', version = aips_version_std)
    count = 0
    for source in sourcelist:
        if v190_utils.wasobserved(uvdata, source):
            count = count + 1
            split.sources[count] = source
            tempdata =  AIPSUVData(source, 'SPLIT', 1, 1)
            if tempdata.exists():
                if v190_utils.yesno(autostomp, 'Delete existing ' + source + ' data from catalogue and proceed (no aborts pipeline)?'):
                    tempdata.zap()
                else:
                    sys.exit()
    split.indata = uvdata
    split.sources[count+1:] = ''
    split.qual = -1
    split.calcode = ''
    split.timerang[1:] = [0]
    split.stokes = ''
    split.selband = -1
    split.selfreq = -1
    split.freqid = -1
    split.bif = 0
    split.eif = 0
    if doaverage:
        nchan = uvdata.table('FQ', 1)[0].total_bandwidth[0]/uvdata.table('FQ', 1)[0].ch_width[0]
        split.bchan = nchan/10
        split.echan = nchan - nchan/10
        split.aparm[1] = 2
    else:
        split.bchan = 0
        split.echan = 0
        split.aparm[1] = 1
        split.nchav = 1
        split.chinc = 1
    split.aparm[2:] = [0]
    split.subarray = 0
    split.docalib = 2
    split.gainuse = clversion
    split.dopol = -1
    split.blver = -1
    split.flagver = 1
    split.doband = -1
    split.bpver = 0
    if doband:
        #8 GHz stuff, so do bandpass correction
        try:
            bptable = uvdata.table('BP', 1)
            for row in bptable: break #Basically a test to see if it exists...
            split.doband = 4
            split.bpver = 1
        except IOError:
            print "BP table 1 does not exist - doband will be set to false!"
            doband = False
    split.smooth[1:] = [0]
    split.outclass = 'SPLIT'
    split.outseq = 1
    split.outdisk = 1
    split.douvcomp = -1
    split.ichansel[1:][1:] = [0]
    split.baddisk[1:] = [0]
    split()
    
################################################################################
# Shifts a list of sources by the specified amount
def uvshift(sourcelist, shift_arcseconds, outversion):
    at = 0
    utcvals = get_calc_utcvals(uvdata)
    uvfix = AIPSTask('uvfix', version = 'ADAMAIPS')
    uvfix.uvfixprm[1:] = [0]
    uvfix.uvfixprm[11] = utcvals[0]
    uvfix.uvfixprm[12] = utcvals[1]
    uvfix.uvfixprm[17] = 2000.0
    uvfix.uvfixprm[20] = 1
    for source in sourcelist:
        if v190_utils.wasobserved(uvdata, source):
            uvfix.indata =  AIPSUVData(source, 'SPLIT', 1, 1)
            outputdata = AIPSUVData(source, 'SPLIT', 1, outversion)
            if outputdata.exists():
                if v190_utils.yesno(autostomp, 'Delete existing ' + source + ' data from catalogue and proceed (no aborts pipeline)?'):
                    outputdata.zap()
                else:
                    sys.exit()
            uvfix.outdata =  outputdata
            uvfix.shift[1] = shift_arcseconds[at][0]
            uvfix.shift[2] = shift_arcseconds[at][1]
            uvfix()
        at = at+1
#            wizardry_uvshift.uvshift(source, 'SPLIT', 1, 1, shift_arcseconds[at][0]*math.pi/648000, shift_arcseconds[at][1]*math.pi/648000)


################################################################################
# Averages the list of sources down to a single channel
def average_frequency(sourcelist, inversion, outversion):
    for source in sourcelist:
        if v190_utils.wasobserved(uvdata, source):
            sourcedata = AIPSUVData(source, 'SPLIT', 1, inversion)
            outputdata = AIPSUVData(source, 'SPLIT', 1, outversion)
            if outputdata.exists():
                if v190_utils.yesno(autostomp, 'Delete existing ' + source + ' data from catalogue and proceed (no aborts pipeline)?'):
                    outputdata.zap()
                else:
                    sys.exit()
            avspc = AIPSTask('avspc', version = aips_version_std)
            avspc.indata = sourcedata
            avspc.outdata = outputdata
            nchan = uvdata.table('FQ', 1)[0].total_bandwidth[0]/uvdata.table('FQ', 1)[0].ch_width[0]
            avspc.ichansel[1][1] = nchan/10
            avspc.ichansel[1][2] = nchan - nchan/10
            avspc.ichansel[1][3] = 1
            avspc.ichansel[2:][1:] = [0]
            avspc()

################################################################################
# Averages the list of sources in time to the specified time
def average_time(sourcelist, timescales, inversion, outversion):
    uvavg = AIPSTask('uvavg', version = aips_version_std)
    count = 0
    for source in sourcelist:
        if v190_utils.wasobserved(uvdata, source):
            sourcedata = AIPSUVData(source, 'SPLIT', 1, inversion)
            outputdata = AIPSUVData(source, 'SPLIT', 1, outversion)
            if outputdata.exists():
                if v190_utils.yesno(autostomp, 'Delete existing ' + source + ' data from catalogue and proceed (no aborts pipeline)?'):
                    outputdata.zap()
                else:
                    sys.exit()
            uvavg.indata = sourcedata
            uvavg.outdata = outputdata
            uvavg.xinc = 1
            uvavg.yinc = timescales[count]*60.0
            uvavg.zinc = 0
            uvavg.opcode = '    '
            uvavg.optype = '    '
            uvavg()
        count = count + 1
        
################################################################################
# Copies a table from one version of a uv dataset to another
def copy_table(pulsar, aipsclass, inseq, outseq, tabletype, tableversion, outtabversion, outdata):
    tacop = AIPSTask('tacop', version = aips_version_std)
    tacop.indata = AIPSUVData(pulsar, aipsclass, 1, inseq)
    tacop.inext = tabletype
    tacop.ncount = 1
    tacop.invers = tableversion
    if outdata == None:
        tacop.outdata = AIPSUVData(pulsar, aipsclass, 1, outseq)
    else:
        tacop.outdata = outdata
    tacop.outvers = outtabversion
    tacop.keyword = ''
    tacop.keyvalue[1:] = [0]
    tacop.keystrng = ''
    tacop()
    
################################################################################
# Applies an SN table using splat
def apply_sn_table(pulsar, aipsclass, snversion, inseq, outseq):
    splat = AIPSTask('splat', version = aips_version_std)
    splat.indata =  AIPSUVData(pulsar, aipsclass, 1, inseq)
    outputdata = AIPSUVData(pulsar, aipsclass, 1, outseq)
    splat.outdata = outputdata
    if outputdata.exists():
        if v190_utils.yesno(autostomp, 'Delete existing ' + pulsar + ' data from catalogue and proceed (no aborts pipeline)?'):
            outputdata.zap()
        else:
            sys.exit()
    splat.sources[1:] = ''
    splat.qual = -1
    splat.calcode = ''
    splat.timer[1:] = [0]
    splat.stokes = ''
    splat.selband = -1
    splat.selfreq = -1
    splat.freqid = -1
    splat.bif = 0
    splat.eif = 0
    splat.bchan = 0
    splat.echan = 0
    splat.subarray = 0
    splat.docalib = 2
    splat.gainuse = snversion
    splat.dopol = -1
    splat.blver = -1
    splat.flagver = 0
    splat.doband = -1
    splat.bpver = 0
    splat.smooth[1:] = [0]
    splat.douvcomp = -1
    splat.aparm[1] = 1
    splat.aparm[2:] = [0]
    splat.aparm[4] = 1
    splat.ichansel[1:][1:] = [0]
    splat.channel = 0
    splat.solint = 0
    splat.baddisk[1:] = [0]
    splat()    
    
################################################################################
# Writes a list of sources to disk in uvfits format
def write_to_disk(sourcelist, outclass, sequence, zapafter):
    for source in sourcelist:
        if v190_utils.wasobserved(uvdata, source):
            sourcedata = AIPSUVData(source, outclass, 1, sequence)
            fittp = AIPSTask('fittp', version = aips_version_std)
            fittp.indata = sourcedata
            fittp.outfile = rootdir + '/' + experiment + '/' + experiment + '_' + source + '_uv.fits'
            fittp.dostokes = -1
            fittp.donewtab = 1
            fittp.format = 0
            fittp()
            #if(zapafter):
            #    sourcedata.zap()

################################################################################
# Replaces the CAT* in the global antennas array with correct ATCA reference
def replace_atca_name(uvdata):
    global calibrator_instrument_search

    count = 0
    for antenna in calibrator_instrument_search:
        if antenna == 'CAT*':
            replaceindex = count
        count = count+1
    
    for antenna in uvdata.antennas:
        if antenna[0:3] == 'CAT':
            calibrator_instrument_search[replaceindex] = antenna
            return antenna
    return ''
            
################################################################################
# Update the predicted flag file produced by sched with more accurate online
# info where available
def update_skdflag():
    skdflagfile = v190_utils.check_data_file(flagprefix, '', skdflagsuffix)
    skdflagin = open(skdflagfile,'r')
    skdlines = skdflagin.readlines()
    skdflagin.close()
    edited = 0
    
    acount = 0
    # test the availability of an online flag file for each antenna
    for ant in twochar_antenna_names:
        try:
            onlinefile = v190_utils.check_data_file(flagprefix, ant, onlineflagsuffix)
            print 'Merging online flagging file ' + onlinefile
            edited = 1
            count = 0
            mincount = 9999999
            maxcount = 0
            
            #Grab the contents of the online flag file
            onlineflagin = open(onlinefile,'r')
            onlinelines = onlineflagin.readlines()
            #if ant != 'Ti' and ant != 'Ho':
            #    for line in onlinelines:
            #        line = line.replace('reason=slewing', 'reason=\'slewing\'')
            #        #line = line.replace('/', '\'/')
            #        onlinelines[count] = line
            #        count = count + 1

            count = 0
            #Skip through skd file lines finding lines to omit
            for line in skdlines:
                if line[10:12].upper() == ant.upper() or line[10:13].upper() == calibrator_instrument_search[acount][:3].upper():
                    if mincount == 9999999: mincount = count
                    if count > maxcount: maxcount = count
                count = count+1

            #Create a new array of lines for the flag file based on
            #online flag file + other telescopes
            if mincount > maxcount:
                print 'Didn\'t find any predicted flagging for ' + ant
                skdlines = skdlines + onlinelines
            else:
                skdlines = skdlines[:mincount] + onlinelines + skdlines[maxcount+1:]

            #Close the online flag file
            onlineflagin.close()

            #Change the name of the online flag file to reflect fact
            #its been inserted
            os.rename(onlinefile, onlinefile + '.merged')
        except AssertionError:
            print 'No online flag file for ' + ant
        acount = acount + 1

    if(edited):
        #Dump out the contents of memory into new, updated flag file
        skdflagout = open(skdflagfile, 'w')
        skdflagout.writelines(skdlines)
    

################################################################################
# Writes an SN table out to disk and reads it back in again
def safeguard_sntable(uvdata, snver):
    tempfile = rootdir + '/tempjunk.sn'
    if os.path.exists(tempfile):
        os.system("rm -f " + tempfile)
    tbout = AIPSTask('tbout', version = aips_version_std)
    tbout.indata = uvdata
    tbout.inext = 'SN'
    tbout.inver = snver
    tbout.outfile = tempfile
    tbout.docrt = 132
    tbout.bcount = 1
    tbout.ecount = 0
    tbout()
    sntable = uvdata.table('SN', snver)
    sntable.zap()
    tbin = AIPSTask('tbin', version = aips_version_std)
    tbin.outdata = uvdata
    tbin.infile = tempfile
    tbin.bcount = 1
    tbin.ecount = 0
    tbin()
    #os.system("rm -f " + tbout.outfile)

################################################################################
# Fits a single model component using JMFIT and write to the jmfit summary file
def jmfit(pulsar, isll):
    #Do for each IF/pol, then the combined
    todo = ['rr.1', 'll.1', 'rr.2', 'll.2', 'rr.3', 'll.3', 'rr.4', 'll.4', 'combined']
    imagedata = AIPSImage('JUNK', 'IMG', 1, 1)
    outdata = AIPSImage('JUNK', 'IMG', 1, 2)
    psrstatout = open(rootdir + "/" + experiment + "/" + pulsar + ".jmfitstats", "w")
    for ifpol in todo:
        #First load the image data from difmap
        if imagedata.exists():
            imagedata.zap()
        if outdata.exists():
            outdata.zap()
        fitld = AIPSTask('fitld', version = aips_version_std)
        fitld.infile = rootdir + '/' + experiment + '/' + pulsar + '.' + ifpol + '.image.fits'
        fitld.outdata = imagedata
        fitld.optype = ''
        fitld.ncount = 0
        fitld.dotable = 1
        fitld.douvcomp = -1
        fitld.doconcat = -1
        fitld()
        print "Imagedata.header[1] is " + str(imagedata.header.crval[1]) + ", " + str(imagedata.header.crpix[1])
        print "Imagedata.header[2] is " + str(imagedata.header.crval[2]) + ", " + str(imagedata.header.crpix[2])
    
        jmfit = AIPSTask('jmfit')
        jmfit.indata = imagedata
        jmfit.outdata = outdata
        jmfit.niter = 4000
        jmfit.ngauss = 1
        jmfit.fwidth[1][1] = 0
        jmfit.fwidth[1][2] = 0
        jmfit.fwidth[1][3] = 0
        jmfit.domax[1:] = [1]
        jmfit.dopos[1][1] = 1
        jmfit.dopos[1][2] = 1
        jmfit.dowidth[1][1] = 1
        jmfit.dowidth[1][2] = 1
        jmfit.dowidth[1][3] = 1
        jmfit.blc[1:] = [int(imagedata.header.crpix[0])-40, int(imagedata.header.crpix[1])-40]
        jmfit.trc[1:] = [int(imagedata.header.crpix[0])+40, int(imagedata.header.crpix[1])+40]
        jmfit.go()

        jmfitmessage = jmfit.message()
        msgindex = 0
        exciselinenos = []
        for i in range(len(jmfitmessage)-1):
            if jmfitmessage[len(jmfitmessage)-(i+1)][:5] != 'JMFIT':
                exciselinenos.append(len(jmfitmessage)-(i+1))
        for e in exciselinenos:
            jmfitmessage = jmfitmessage[:e] + jmfitmessage[e+1:]
        while jmfitmessage[msgindex].find('X-ref') == -1 and msgindex < len(jmfitmessage):
            msgindex = msgindex + 1
        if msgindex >= len(jmfitmessage):
            print "Did not find solution for " + pulsar + ", band " + ifpol
            continue
        centrerasplit = jmfitmessage[msgindex].split()
	centredecsplit = jmfitmessage[msgindex+1].split()
        msgindex = msgindex + 2
        while jmfitmessage[msgindex].find('********* Solution from JMFIT') == -1 and msgindex < len(jmfitmessage):
            msgindex = msgindex + 1
        if msgindex >= len(jmfitmessage):
            print "Did not find solution for " + pulsar + ", band " + ifpol
            continue
        sourcerasplit = jmfitmessage[msgindex+7].split()
        sourcedecsplit = jmfitmessage[msgindex+8].split()
        fluxsplit = jmfitmessage[msgindex+3].split()
        beammajsplit = jmfitmessage[msgindex+12].split()
        beamminsplit = jmfitmessage[msgindex+13].split()
        beampasplit = jmfitmessage[msgindex+14].split()
        try:
            beammaj = float(beammajsplit[4])
            beammin = float(beamminsplit[4])
        except ValueError:
            print "Bad value for one or more of major axis " +  beammajsplit[4] + \
                  ", minor axis " + beamminsplit[4] + ", both will be set to zero"
            beammaj = 0.0
            beammin = 0.0
        
        #declination =  (float(jmfit.message()[4][29:31])*60 + float(jmfit.message()[4][32:34]))*60 + float(jmfit.message()[4][35:43])
        #ra_centre = float(jmfit.message()[4][35:43])
        #dec_centre = float(jmfit.message()[5][36:43])
        #ra_source = float(jmfit.message()[36][39:47]) #in seconds
        #dec_source = float(jmfit.message()[37][40:47]) #in arcseconds
        #delta_ra = ra_centre - ra_source
        #delta_dec = -dec_centre + dec_source
        if ifpol != 'combined':
            psrstatout.write("Pulsar " + pulsar + ": Frequency " + ifpol[-1] + ", pol " + ifpol[:2].upper() + ":\n")
        else:
            psrstatout.write("Pulsar " + pulsar + " combined results:\n")
        psrstatout.write("Centre RA:            " + centrerasplit[5] + ":" + centrerasplit[6] + ":" + \
                         centrerasplit[7] + "\n")
        
        psrstatout.write("Centre Dec:           " + centredecsplit[5] + ":" + centredecsplit[6] + ":" + \
                         centredecsplit[7] + "\n")
        if fluxsplit[4] == '+/-':
            flux = float(fluxsplit[3][1:])*1000.0
            rms = float(fluxsplit[5])*1000.0
        else:
            flux = float(fluxsplit[4])*1000.0
            rms = float(fluxsplit[6])*1000.0
        psrstatout.write("Flux (mJy):           " + str(flux) + "\n")
        psrstatout.write("S/N:                  " + str(flux/rms) + "\n")
        centrerahours = float(centrerasplit[5]) +  float(centrerasplit[6])/60.0 +  float(centrerasplit[7])/3600.0
        centredecdegs = float(centredecsplit[5]) +  float(centredecsplit[6])/60.0 +  float(centredecsplit[7])/3600.0
        srcrahours = float(sourcerasplit[2]) +  float(sourcerasplit[3])/60.0 +  float(sourcerasplit[4])/3600.0
        srcdecdegs = float(sourcedecsplit[2]) +  float(sourcedecsplit[3])/60.0 +  float(sourcedecsplit[4])/3600.0
        raoff = (srcrahours-centrerahours)*15*60*60*1000*math.cos(centredecdegs*math.pi/180.0)
        decoff = (srcdecdegs-centredecdegs)*60*60*1000
        psrstatout.write("RA offset (mas):      " + str(raoff) + "\n")
        psrstatout.write("Dec offset (mas):     " + str(decoff) + "\n")
        psrstatout.write("Actual RA:            " + sourcerasplit[2] + ":" + sourcerasplit[3] + ":" + \
                         sourcerasplit[4] + "\n")
        psrstatout.write("Actual Dec:           " + sourcedecsplit[2] + ":" + sourcedecsplit[3] + ":" + \
                         sourcedecsplit[4] + "\n")
        psrstatout.write("Beam:                 " + str(beammin*1000) + "x" + \
                         str(beammaj*1000) + " at " + beampasplit[4] + " degrees\n")
        raerr = float(sourcerasplit[6])*1000
        decerr = float(sourcedecsplit[6])*1000
        psrstatout.write("Est. RA error (mas):  " + str(raerr*math.cos(centredecdegs*math.pi/180.0)*15.0) + "\n")
        psrstatout.write("Est. RA error (hms):  " + str(raerr) + "\n")
        psrstatout.write("Est. Dec error (mas): " + str(decerr) + "\n")
    psrstatout.close()

################################################################################
################################################################################
################################################################################
################ MAIN CODE #####################################################
################################################################################
################################################################################

# Load the data into AIPS
if skipmanager.skipload and clversion == 1:
    raw_input("Remember to use --clversion if necessary - press any key to continue or ctrl-c if you stuffed up!!")
if not skipmanager.skipload:
    load_rpfits(nfits, rpfitsfilenames, uvdata, doband, douvsrt)

if douvsrt:
    #Sort into TB order (usually not necessary)
    if not skipmanager.skipuvsrt:
        tb_sort()
    uvdata = AIPSUVData(experiment.upper(), 'UVSRT', 1, 1)

# Find the reference antenna for the ATCA and replace CAT* with this name
atcaname = replace_atca_name(uvdata)

#Update the weights to reflect telescope sensitivities
if not skipmanager.skipsensweight:
    print "Overwriting CL table " + str(clversion) + " weights!!!"
    wizardry_sensweight.fix_weights(experiment, uvdata.klass, 1, 1)

if not skipmanager.skipmodelcorrect:
    if os.path.exists(rootdir + '/' + experiment + '/applied/'):
        print "Correcting model!"
        wizardry_correctmodel.correct_model(rootdir, experiment, uvdata.klass, True,
                                            skipmanager.skipmccheck)
        clversion = clversion+1
    else:
        print "Skipping model-correction, since there is no applied/ directory"

# Edit the skdflag file with more accurate online flagging for the telescopes
# where this is available
if not skipmanager.skipflag:
    update_skdflag()

# Flag the data
if not skipmanager.skipflag:
    load_flagfile(uvdata, False)
if not skipmanager.skipusedflags:
    load_flagfile(uvdata, True)
# Quack the first and last 10 seconds of every scan, and the first 30 seconds of every ATCA scan
if not skipmanager.skipquack:
    default_quack(uvdata, atcaname)

# Correct for ionosphere using Tecor (CL 1->2 or 2->3)
if not skipmanager.skiptecor:
    if experiment == 'v190a' or experiment == 'v190e' or experiment == 'v190g' or experiment == 'v190k':
        print "TECOR is permanently switched off for 8 GHz experiments!"
    else:
        tecor_ionosphere(uvdata)
        clversion = clversion + 1

# Load the tsys data into SN 1 - call another python module to do the translation/concatenation
if not skipmanager.skiptsys:
#    if (not os.path.isfile(tsysprefix + '/' + experiment + tsyssuffix)) or \
#       (not skipmanager.skipantabgen):
#        logs2antab.create_antab(rootdir, experiment)
#        if not skipmanager.skipantabgen:
#        if v190_utils.yesno(False, 'Recreate experiment antab file (yes) or use existing (no)'):
#            logs2antab.create_antab(rootdir, experiment)
#    else:
#        logs2antab.create_antab(rootdir, experiment)
    #load_tsys(uvdata)
    load_tsys_sn(uvdata, rootdir+'/'+tsysprefix+'/'+experiment+'.tsys.sn')

# Fringe fit using calibrator models into SN 2
for source in uvdata.sources:
    if v190_utils.isphaseref(source, fringefinders):
        if not skipmanager.skipfring:
            print "Fringe fitting " + source + " using cltable " + str(clversion)
            fringe_fit(uvdata, source, clversion, False)
    if v190_utils.isphaseref(source, primary_phaserefs):
        if not skipmanager.skipfring:
            print "Fringe fitting " + source + " using cltable " + str(clversion)
            fringe_fit(uvdata, source, clversion, False)
            
#If we will be snediting fring data, safeguard it
if not skipmanager.skipsnedt and (not skipmanager.skipcal_f):
    safeguard_sntable(uvdata, 2)

# Use SNEDT to get rid of any bad SN points
#if not skipmanager.skipsnedt and (not skipmanager.skipcal_t):
#    tvedit_sndata(uvdata, 1)
if not skipmanager.skipcal_t:
    tacop = AIPSTask('tacop', version = aips_version_std)
    tacop.indata = uvdata
    tacop.inext = 'SN'
    tacop.ncount = 1
    tacop.invers = 1
    tacop.outdata = uvdata
    tacop.outver = 3
    tacop.keyword = ''
    tacop.keyvalue[1:] = [0]
    tacop.keystrng = ''
    tacop()
if not skipmanager.skipsnedt and (not skipmanager.skipcal_f):
    tvedit_sndata(uvdata, 2)

# CLCAL for the tsys and fringe-fit SN tables, producing plots of the
# amplitude, phase, delay and rate corrections (CL 2->3->4 or 3->4->5)
calibrate_tsys_fringefit(uvdata, skipmanager.skipcal_t, skipmanager.skipcal_f, skipmanager.skipplot)

# Use EDITR to have a quick squiz for any obviously bad data (applying calibration & flagging)
if not skipmanager.skipuvedt:
    tvedit_uvdata(uvdata, 4)

# Self-cal on the calibrator models and generate corrections using calib (creates SN 3[5])
if not skipmanager.skipselfcal:
    expstartdayfrac = exp_start_hours[ord(experiment[-1]) - ord('a')]/24.0
    #Split the calibrators and write
    split(uvdata, primary_phaserefs, 1, clversion, doband)
    #Fix the weights if desired
    if not skipmanager.skipfixweights:
        for ref in primary_phaserefs:
            if v190_utils.wasobserved(uvdata, ref):
                wizardry_weightfix.fixweights(rootdir, experiment, ref, uvdata.klass, uvdata.seq, 1, 'SPLIT')
    
    write_to_disk(primary_phaserefs, 'SPLIT', 1, True)

    #basesntabno = 3
    #if not skipmanager.skipsnedt:
    #    basesntabno = 5
    basesntabno = 5 #ALWAYS! Now.
    
    #Go through and use difmap to selfcal/cordump
    calno = 0
    isll = (int((uvdata.table('FQ', 1)[0].if_freq[0] + uvdata.header.crval[2])/100000000)==16) #18cm, therefore ll only
    if experiment == "v190j":
        isll = False;
    print "IF Freq is " + str((uvdata.table('FQ', 1)[0].if_freq[0] + uvdata.header.crval[2])/100000000) + ", in 100s of MHz: " + \
          str(int((uvdata.table('FQ', 1)[0].if_freq[0] + uvdata.header.crval[2])/100000000)) + " so isll is " + str(isll)
    for source in uvdata.sources:
        if v190_utils.isphaseref(source, primary_phaserefs):
            print "About to selfcal source " + source
            sourceindex = 0
            for p in primary_phaserefs:
                if p == source:
                    break
                sourceindex = sourceindex + 1
            #aips_selfcal(uvdata, source)
            #selfcal, cordump, then load SN table to 4[6] + calno, then deletes from disk
            
            difmap_selfcal(uvdata, isll, source, sourceindex, calno, basesntabno, expstartdayfrac)
            if not skipmanager.skipsnedt:
                tvedit_sndata(AIPSUVData(source, 'SPLIT', 1, 1), 1)
                copy_table(source, 'SPLIT', 1, 0, 'SN', 2, basesntabno+calno, uvdata)
            else:
                copy_table(source, 'SPLIT', 1, 0, 'SN', 1, basesntabno+calno, uvdata)
            calibrate_selfcalcorrections(uvdata, basesntabno+calno, skipmanager.skipplot, sourceindex)
            calno = calno + 1
    clversion = clversion + 1

    #Merge the SN tables into SN 3[5], delete 4[6] to 4[6]+numcals-1
    #mergesn(uvdata, primary_phaserefs, calno, basesntabno)
    
    # Do SNEDT again for selfcal corrections if required
    #if not skipmanager.skipsnedt:
    #    safeguard_sntable(uvdata, basesntabno+calno+1)
    #    tvedit_sndata(uvdata, basesntabno+calno+1)
    #
    #    # CLCAL for the selfcal corrections from the phase reference calibrators
    #    # (CL 4->5or 5->6), producing plots of the amplitude and phase corrections in turn
    #    calibrate_selfcalcorrections(uvdata, basesntabno+calno+2, skipmanager.skipplot, -1)
    #else:
    #    calibrate_selfcalcorrections(uvdata, basesntabno+calno+1, skipmanager.skipplot, -1)

#isll = True
#tvedit_sndata(uvdata, 14)
#calibrate_selfcalcorrections(uvdata, 15, skipmanager.skipplot)

# Now create all the flag files, if required
if not skipmanager.skipcreatedflags:
    if load_junkfiles(uvdata, primary_phaserefs) > 0:
        #Work out and provide sequence numbers, zap afterwards
        seqnos = []
        for seqno in range(99):
            if AIPSUVData(experiment.upper(), 'DFLAG', 1, seqno).exists():
                seqnos.append(seqno)
        wizardry_createdifmapflagfile.create_difmapflagfile(experiment.upper(), 'DFLAG', flagprefix, False,
                                                            experiment, uvdata.klass, v190_utils.getdoy(uvdata.header.date_obs),
                                                            seqnos)
        load_flagfile(uvdata, True)
        for seqno in seqnos:
            AIPSUVData(experiment.upper(), 'DFLAG', 1, seqno).zap()
        #junkdata = AIPSUVData(experiment.upper(), 'DFLAG', 1, 99)
        #junkdata.zap()
difmapsuffixes = ['.par','.uvf','.fits','.mod','.win']
for source in primary_phaserefs:
    if v190_utils.wasobserved(uvdata, source):
        for suff in difmapsuffixes:
            if os.path.exists(rootdir + '/' + experiment + '/junk' + source + suff):
                os.remove(rootdir + '/' + experiment + '/junk' + source + suff)

# Split all the calibrators, with averaging
if not skipmanager.skipwriterefs:
    split(uvdata, primary_phaserefs, 1, clversion, doband)

    #Fix the weights if desired
    if not skipmanager.skipfixweights:
        for ref in primary_phaserefs:
            if v190_utils.wasobserved(uvdata, ref):
                wizardry_weightfix.fixweights(rootdir, experiment, ref, uvdata.klass, uvdata.seq, 1, 'SPLIT')

    # Image and write out uv fits for all the calibrators, zap them afterwards
    write_to_disk(primary_phaserefs, 'SPLIT', 1, True)
    isll = (int((uvdata.table('FQ', 1)[0].if_freq[0] + uvdata.header.crval[2])/100000000)==16) #18cm, therefore ll only
    if experiment == "v190j":
        isll = False;
    for source in primary_phaserefs:
        if v190_utils.wasobserved(uvdata, source):
            difmap_mapcal(source, isll)
            jmfit(source, isll)

final_target_seqno = 3
# Split all the targets, without averaging
if not skipmanager.skipuvshift or not skipmanager.skipscintweight:
    split(uvdata, target_pulsars, 0, clversion, doband)
else:
    if not skipmanager.skipwritetargets:
        split(uvdata, target_pulsars, 1, clversion, doband)
        final_target_seqno = 1

#epoch_offset is fractional years since jan 1, 2006
epoch_offset = float(v190_utils.header_mjd(uvdata.header.date_obs) - 53736)/365.25
for i in range(len(target_pulsars)):
    pulsar_best_mar2006_shifts[i][0] = pulsar_best_mar2006_shifts[i][0]+pulsar_best_pm[i][0]*epoch_offset/1000.0
    pulsar_best_mar2006_shifts[i][1] = pulsar_best_mar2006_shifts[i][1]+pulsar_best_pm[i][1]*epoch_offset/1000.0

# Shift to the best guess position and work out weights
if not skipmanager.skipscintweight:
    final_target_seqno = 5
    
    #Uvshift to the best guess position
    #epoch_offset is fractional years since jan 1, 2006
    uvshift(target_pulsars, pulsar_best_mar2006_shifts, 2)
        
    #Average in frequency, then in time
    average_frequency(target_pulsars, 2, 3)

    average_time(target_pulsars, pulsar_scintillation_time, 3, 4)
    
    #Create an SN table with gain based on the scintillating amplitudes
    pcount = 0
    exclude = [0,0,0] #antenna numbers of HOB, MOP, CED
    try:
        exclude[0] = v190_utils.get_antenna_number('HOB', uvdata)
    except RuntimeError:
        print 'Hobart not used in this experiment'
    try:
        exclude[1] = v190_utils.get_antenna_number('MOPRA', uvdata)
    except RuntimeError:
        print 'Mopra not used in this experiment'
    try:
        exclude[2] = v190_utils.get_antenna_number('CED', uvdata)
    except RuntimeError:
        print 'Ceduna not used in this experiment'

    for pulsar in target_pulsars:
        if v190_utils.wasobserved(uvdata, pulsar) and pulsar_scintillation_time[pcount] > 0:
            wizardry_scintamp.create_scint_sn_table(pulsar, 'SPLIT', 1, 4, experiment, uvdata.klass, 2, pulsar_scintillation_time[pcount], exclude)
        pcount = pcount + 1
        
if not skipmanager.skipuvshift:    
    # Average all the targets in frequency
    if not skipmanager.skipscintweight:
        # Uvshift all the targets to the desired position, then average
        #uvshift(target_pulsars, target_shifts, final_target_seqno-2)
        
        #average_frequency(target_pulsars, final_target_seqno-2, final_target_seqno-1)
        
        pcount = 0
        for pulsar in target_pulsars:
            if pulsar_scintillation_time[pcount] > 0 and v190_utils.wasobserved(uvdata, pulsar):
                #Copy the amplitude table from that copy to here
                copy_table(pulsar, 'SPLIT', 4, 3, 'SN', 1, 1, None)
                #use splat to apply the SN table created based on amplitudes
                apply_sn_table(pulsar, 'SPLIT', 1, final_target_seqno-2, final_target_seqno)
                pcount = pcount + 1
    else:
        # Uvshift all the targets to the desired position, then average
        uvshift(target_pulsars, pulsar_best_mar2006_shifts, final_target_seqno-1)
        average_frequency(target_pulsars, final_target_seqno-1, final_target_seqno)

if not skipmanager.skippmcor:
    for pulsar in target_pulsars:
        if v190_utils.wasobserved(uvdata, pulsar):
            finaldata = AIPSUVData(pulsar, 'SPLIT', 1, final_target_seqno+1)
            if finaldata.exists():
                if v190_utils.yesno(autostomp, 'Delete existing ' + pulsar + ',' + str(final_target_seqno+1) + ' data from catalogue and proceed (no aborts pipeline)?'):
                    finaldata.zap()
                else:
                    sys.exit()
            wizardry_correctpropermotion_oneshot.correct_propermotion(rootdir, experiment, pulsar,
                                                                      uvdata, final_target_seqno,
                                                                      'SPLIT')
    final_target_seqno = final_target_seqno + 1

if not skipmanager.skipfixweights:
    for pulsar in target_pulsars:
        if v190_utils.wasobserved(uvdata, pulsar):
            wizardry_weightfix.fixweights(rootdir, experiment, pulsar, uvdata.klass, uvdata.seq, final_target_seqno, 'SPLIT')

if not skipmanager.skipwritetargets:
    # Write out uv fits for all the targets
    write_to_disk(target_pulsars, 'SPLIT', final_target_seqno, False)
    isll = (int((uvdata.table('FQ', 1)[0].if_freq[0] + uvdata.header.crval[2])/100000000)==16) #18cm, therefore ll only
    if experiment == "v190j":
        isll = False;
    for source in target_pulsars:
        if v190_utils.wasobserved(uvdata, source):
            pdata = AIPSUVData(source, 'SPLIT', 1, final_target_seqno)
            difmap_mappsr(source, isll, pdata.header.crval[4], pdata.header.crval[5])
            jmfit(source, isll)

    # Now create all the flag files, if required
    if not skipmanager.skipcreatedflags:
        if load_junkfiles(uvdata, target_pulsars) > 0:
            #Work out and provide sequence numbers, zap afterwards
            seqnos = []
            for seqno in range(99):
                if AIPSUVData(experiment.upper(), 'DFLAG', 1, seqno).exists():
                    seqnos.append(seqno)
            wizardry_createdifmapflagfile.create_difmapflagfile(experiment.upper(), 'DFLAG', flagprefix, False,
                                                            experiment, uvdata.klass, v190_utils.getdoy(uvdata.header.date_obs),
                                                            seqnos)
            load_flagfile(uvdata, True)
            for seqno in seqnos:
                AIPSUVData(experiment.upper(), 'DFLAG', 1, seqno).zap()
            #junkdata = AIPSUVData(experiment.upper(), 'DFLAG', 1, 99)
            #junkdata.zap()
        difmapsuffixes = ['.par','.uvf','.fits','.mod','.win']
        for source in target_pulsars:
            if v190_utils.wasobserved(uvdata, source):
                for suff in difmapsuffixes:
                    if os.path.exists(rootdir + '/' + experiment + '/junk' + source + suff):
                        os.remove(rootdir + '/' + experiment + '/junk' + source + suff)

if not options.notv:
    tv.kill()
# Zap all the unaveraged targets and the unshifted targets
#for target in target_pulsars:
#    for i in range(final_target_seqno - 1):
#        zapdata = AIPSUVData(target, 'SPLIT', 1, i+1)
#        if zapdata.exists():
#            zapdata.zap()