#!/usr/bin/env ParselTongue
# Pipeline script for processing EVN data
# Cormac Reynolds (reynolds@jive.nl): Jan 2006 original version
from AIPS import AIPS, AIPSDisk
from AIPSTask import AIPSTask
from AIPSData import AIPSUVData, AIPSImage
#from AIPSTV import AIPSTV
#import Wizardry.AIPSData
from evn_aips_tasks import *
from evn_funcs import *
import evn_funcs
import evn_aips_tasks
import crfuncs
import traceback

import LocalProxy
from xmlrpclib import ServerProxy

import copy, optparse, os, sys
import re, string, numarray, numarray.mlab, pprint, math
import time

#default_ver = '31DEC05'
#alt_ver = '31DEC06'

#AIPSTask.version = default_ver
#AIPSTask.isbatch = 0


c = 2.99792458e8

# define our local functions. More general purpose functions can be found in
# evn_funcs.py. The actual AIPS task calls can be found in evn_aips_tasks.py.

def getfile(extension):
    # get calibration files from indir assuming standard names.
    file = indir + '/' + experiment + '.' + extension
    if not os.path.isfile(file):
        file = indir + '/' + experiment + '.' + extension
        #assert(os.path.isfile(file)), file + ' does not exist!'
    return file



def checkin(control):
    # convert the control hash to a list of variables of the appropriate type.
    # Should probably add some more validity checking here. Make these global
    # for ease of reference - they should never change.
    global indir, refantnames, plotrefnames, fits, outdir
    global nfits, heads, solint, experiment
    global bpass_calibrators, phaseref_sources, target_sources, avg, plotavg
    global tmask, fitsdir, fits_file, nselfcal, doplot, msgkill
    global source_select, experiment_dir, output_prefix, nbp_table, fring_snr
    global disk, freqid, dopng

    evn_aips_tasks.AIPSTask.version = control['version']['default']
    evn_aips_tasks.task_versions = control['version']

    experiment = control['experiment'][0].lower()
    # remove the pass number from the experiment name if required
    experiment_dir = get_experiment_dir(experiment)

    refantnames = control['refant']
    for i in range(len(refantnames)):
        refantnames[i] = refantnames[i].upper()

    plotrefnames = control['plotref']
    for i in range(len(plotrefnames)):
        plotrefnames[i] = plotrefnames[i].upper()

    fitsdir = control.get('fitsdir', [False])[0]
    if (not fitsdir):
        fitsdir = archive_fitsdir(experiment_dir)

    outdir = control.get('outdir', [False])[0]
    if (not outdir):
        outdir = os.environ.get('OUT') + '/' + experiment_dir

    indir = control.get('indir', [False])[0]
    if (not indir):
        indir = os.environ.get('IN') + '/' + experiment_dir

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

    # heads defaults to 1
    heads = control.get('glue_pass', [1])
    for i in range(len(heads)):
        heads[i] = int(heads[i])
    if (len(heads) >= 4):
        raise """Too many passes specified (max = 4) - check 'glue_pass' in input file"""

    # solint defaults to the typical scan length on the phaseref source
    solint = float(control.get('solint', [0])[0])

    AIPS.userno = int(control['userno'][0])

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

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

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

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

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

    tmask = control['tmask']
    for i in range (len(tmask)):
        tmask[i] = int(tmask[i])
    if (len(tmask) == 1):
        tmask.append(999)

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

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

    # nselfcal defaults to 2 but 0 is a valid input
    nselfcal = int(control.get('sciter', [2])[0])

    doplot = int(control.get('doplot', [True])[0])
    if doplot < 0:
        doplot = 0

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

    disk = int(control.get('disk', [1])[0])

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

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

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

    assert (os.path.exists(fitsdir)), 'FITS directory does not exist: ' + \
            fitsdir

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

    if (not solint and not phaseref_sources):
        raise 'You must set either solint or phaseref_sources'

    output_prefix = outdir + '/' + experiment 

def def_name(name):
    # shortcut to getting the default name for output plot files
    plotname = get_plotname(output_prefix, name)
    return plotname


def save_fits_image(data, name):
    fitsoutfile = outdir + '/' + experiment + '_' + name + '.FITS'

    # Rename any old plot files before writing the new one
    save_old_file(fitsoutfile)        
    runfittp(indata=data, outfile=fitsoutfile)




def print_start(text):
    print >>PYPELOG, '\n'
    print >>PYPELOG, '*' * 20
    print >>PYPELOG, 'Starting tmask ', text, ' at ', time.asctime()

def print_end(text):
    print >>PYPELOG, 'Ending tmask ', text, ' at ', time.asctime()

def print_vplot_err(bif, eif, stokes, uvdata):
    print >>PYPELOG, '\nPype warning, vplot failed for bif=, ', bif,  ', eif = ', \
            eif, ', stokes= ', stokes, ', indata= ',  uvdata


# the following are the main routines of the programme.

def evn_load(uvdata, msortdata, vbgludata, fxpoldata, uvcopdata, uvavgdata):
    # load the data then sort, vbglu and uvavg if required.

    # first zap any of the 'load' files from a previous run (in case the names
    # that are used have changed from the previous run due to e.g. different
    # number of heads)
    zap_old_data(uvdata)
    zap_old_data(msortdata)
    zap_old_data(vbgludata)
    zap_old_data(fxpoldata)
    zap_old_data(uvcopdata)
    zap_old_data(uvavgdata)

    global nfits
    uvdata, headdata = load_data(uvdata, msortdata, experiment, fitsdir,
            fits_file, heads, nfits)

    # and run VBGLU if necessary
    if len(heads) > 1:
        print >>PYPELOG, '\nNeed to VBGLU the data'

        uvdata = glue_data(headdata, vbgludata, fxpoldata)

    # select the freqid if necessary
    if freqid:
        runuvcop(uvdata, uvcopdata, freqid=freqid)
        runindxr(uvcopdata)
        uvdata = uvcopdata


    # now do averaging if requested
    if avg > 0.1e-4:
        runuvavg(indata=uvdata, outdata=uvavgdata, avg=avg)
        runindxr(uvavgdata)
        uvdata = uvavgdata


def evn_checkdata():
    # extract useful information from the dataset

    print >>PYPELOG, '\n', '*' *20
    print >>PYPELOG, 'Checking the data:'
    # Find the reference antennae numbers (this should be a prioritised list)
    refantlist = []
    for antenna in refantnames:
        refantlist.append(get_ant_num(uvdata, antenna))
    print >>PYPELOG, "refant number= ", refantlist

    # and the reference antennae for plotting
    plotref = []
    for antenna in plotrefnames:
        plotref.append(get_ant_num(uvdata, antenna))
    print >>PYPELOG, "plotref number= ", plotref
    #raw_input('hit return')

    uvdata.zap_table('AIPS PL', -1)
    assert(uvdata.table_highver('AIPS PL') == 0)

    #sources = uvdata.sources
    sutable = uvdata.table('SU', 1)
    sources = dict()
    for row in sutable:
        source = row.source.strip()
        if (not source_select) or (source in source_select):
            sources[source] = dict()
            sources[source]['id'] = row.id__no
            sources[source]['scanlen'] = list()

    assert(len(sources) > 0)

    # find out which sources are not targets. Note that selfcal_sources is a
    # superset of phaseref_sources
    for source in sources:
        if not source in (target_sources):
            selfcal_sources.append(source)

    assert(len(selfcal_sources) > 0), 'no fringe calibrators'

    for calibrator in (selfcal_sources):
        assert(calibrator in sources), ('fring calibrator ' + calibrator +  
            ' not in source list')
    for calibrator in bpass_calibrators:
        assert(calibrator in sources), 'bpass calibrator ' + calibrator + ' not in source list'
    for target in target_sources:
        assert(target in sources), 'target source ' + target + ' not in source list'

    print >>PYPELOG, "target sources= ", target_sources
    print >>PYPELOG, "phaseref sources= ", phaseref_sources
    print >>PYPELOG, "selfcal sources= ", selfcal_sources

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

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

    chanwidth = uvdata.header.cdelt[2]  # in Hz
    nchan = uvdata.header.naxis[2]
    nif = uvdata.header.naxis[3]
    ifwidth = chanwidth*nchan

    print >>PYPELOG, 'IF bandwidth= ', ifwidth*1.e-6, ' MHz'
    print >>PYPELOG, 'channel width= ', chanwidth*1.e-3, ' kHz'
    print >>PYPELOG, 'number channels per IF= ', nchan
    print >>PYPELOG, 'number IFs= ', nif
    print >>PYPELOG, 'antennas=', uvdata.antennas

    print >>PYPELOG, "fring SNR cutoff= ", fring_snr

    return (refantlist, sources, selfcal_sources, plotref, ifwidth,
            chanwidth, nchan, nif, stokes, plotstokes, polarizations)


def evn_flag(uvdata):
    # Delete any old flag tables.
    uvdata.zap_table('AIPS FG', -1)

    try:
        #uvflg.infile = uvflg_file
        print >>PYPELOG, "flagging using file:", uvflg_file
        runuvflg(indata=uvdata, infile=uvflg_file)
    except:
        raise "uvflg failed, check " + uvflg_file + " is o.k.!"


    # run vlog if necessary
    doglobal, vlba_missing = vlba_logs(uvdata, vlog_out, vlbacal_file,
            ifwidth)

    # run uvflg for any VLBA antennas.
    if doglobal:
        vlba_uvflg_file = vlog_out + '.FLAG'
        try:
            print >>PYPELOG, "flagging VLBAs using file:", vlba_uvflg_file
            runuvflg(indata=uvdata, infile=vlba_uvflg_file)
        except:
            raise "uvflg failed, check " + vlba_uvflg_file + " is o.k.!"


    if os.path.exists(chflg_file):
        runuvflg(indata=uvdata, infile=chflg_file)
    else:
        # flag the edge channels if no chflag file given
        print >>PYPELOG, "\n", chflg_file + " is missing"

        nflag = min(nchan//16, 8)
        print >>PYPELOG, 'Flagging outer ', nflag, ' channels instead'
        reason = 'subband edge'
        runuvflg(indata=uvdata, bchan=0, echan=nflag, reason=reason)
        runuvflg(indata=uvdata, bchan=nchan+1-nflag, echan=nchan, reason=reason)


def evn_plot1():
    # scan list
    outfile = output_prefix + '.SCAN'
    save_old_file(outfile)        
    runlistr(indata=uvdata, optype='SCAN', outprint=outfile)

    # list of stations and number of vis.
    outfile = output_prefix + '.DTSUM'
    save_old_file(outfile)        
    rundtsum(indata=uvdata, outprint=outfile)

    if doplot:
        # possm plot of the autocorrelations 
        aparm = set_default_aparms('possm')
        # aparm(7) gives freq. labelling, aparm(8) gives autocorrelation
        aparm[7] = 1
        aparm[8] = 1
        #print 'aparm=',aparm

        uvdata.zap_table('AIPS PL', -1)
        runpossm(indata=uvdata, aparm=aparm, stokes=plotstokes)
        plot(uvdata, def_name('POSSM_AUTOCORR'), dopng=dopng)


        # vplot of the raw data
        bparm = set_default_bparms('vplot')
        uvdata.zap_table('AIPS PL', -1)
        for poln in polarizations:
            for incif in range(1, nif+1):
                error = runvplot(indata=uvdata, bparm=bparm, antenna=plotref,
                        bif=incif, eif=incif, docalib=0, stokes=poln+poln,
                        solint=plotavg)
                if (error):
                    print_vplot_err(incif, incif, poln+poln, uvdata)
        plot(uvdata, def_name('VPLOT_UNCAL'), dopng=dopng)

        # possm plot of the cross-correlations 
        aparm = set_default_aparms('possm')
        uvdata.zap_table('AIPS PL', -1)
        runpossm(indata=uvdata, aparm=aparm, stokes=plotstokes,
            antennas=plotref)
        plot(uvdata, def_name('POSSM_UNCAL'), dopng=dopng)

        # possm plot of the cross-hand cross-correlations 
        if len(stokes) > 2:
            aparm = set_default_aparms('possm')
            uvdata.zap_table('AIPS PL', -1)
            runpossm(indata=uvdata, aparm=aparm, stokes='LR', antennas=plotref)
            plot(uvdata, def_name('POSSM_CPOL'), dopng=dopng)

def evn_ampcal(antab_file, uvdata):
    # run ANTAB and APCAL. Also VLOG for any VLBA antennas
    uvdata.zap_table('AIPS SN', -1)
    while uvdata.table_highver('AIPS CL') > 1:
        uvdata.zap_table('AIPS CL', 0)
    uvdata.zap_table('AIPS TY', -1)
    uvdata.zap_table('AIPS GC', -1)

    print >>PYPELOG, 'Antab file= ', antab_file
    runantab(antab_file, uvdata)


    doglobal, sparm = global_exper(uvdata)
    if doglobal:
        # if a global then we must run ANTAB again, but add missing VLBA
        # antennas to SPARM in ANTAB (uggh!).

        antab_file = vlog_out + '.TSYS'
        print >>PYPELOG, 'Antab file= ', antab_file
        runantab(antab_file, uvdata, sparm=sparm)


    uvdata.zap_table('AIPS PL', -1)
    runsnplt(uvdata, 'TY', 1, 'TSYS')
    plot(uvdata, def_name('TSYS'), dopng=dopng)

    runapcal(indata=uvdata)

    # use 'box', dobtween and doblank so no sources get left out
    runclcal(gainver=1, gainuse=2, indata=uvdata, snver=1,
            refant=refantlist[0], interpol='2PT', doblank=1, dobtween=1,
            samptype='BOX')

    uvdata.zap_table('AIPS PL', -1)
    runsnplt(uvdata, 'CL', 2, 'AMP')
    plot(uvdata, def_name('GAIN'), dopng=dopng)

    # Do the parallactic angle correction.
    runclcor(indata=uvdata, clcorprm=[0,1], opcode='PANG', gainver=2, gainuse=2)

def evn_fring(fringdata):
    # run fring (on averaged data set if necessary)
    #max_fring_chan = 128        # max channels for fringe fitting
    #if nchan > max_fring_chan:
    #    avspcdata = AIPSUVData (uvdata.name, 'AVSPC', uvdata.disk, 1)
    #    zap_old_data(avspcdata)
    #    runavspc(indata=uvdata, outdata=avspcdata, nchan=nchan, flagver=1,
    #            channel=int(nchan/max_fring_chan))
    #    fringdata = avspcdata

    ## first get the solution interval from COHER - doesn't work so well
    #for i in range(nif):
    #    runcoher (indata=fringdata, bif=i, sources= phaseref_sources +
    #            selfcal_sources)

    aparm = set_default_aparms('fring')
    dparm = set_default_dparms('fring')
    runfring (indata=fringdata, gainuse=2, refantlist=refantlist, snver=2, 
            solint=solint, calsour=selfcal_sources,
            aparm=aparm, dparm=dparm, snr=fring_snr)
    #runsnsmo(indata=uvdata, refant=refantlist[0], invers=2, smotype='VLBI',
    #        smotime=solint*6.)

    #if nchan > max_fring_chan:
    #    runtacop(indata=avspcdata, outdata=uvdata, inext='SN', inver=2,
    #            outver=2)
    #    avspcdata.zap()

    # apply the fring solutions CL2=>CL3, calibrating the target sources with
    # the appropriate phase calibrators. Zip together the target/phaseref pairs
    # and pass to runclcal. Also need to process the selfcal sources. Note also
    # that calsour and target must be passed to runclcal as lists.
    for (calibrator, target) in ( zip (phaseref_sources, target_sources) +
                    zip (selfcal_sources, selfcal_sources)):
        runclcal(indata=uvdata, opcode='CALI', interpol='AMBG', snver=2,
                gainver=2, gainuse=3, calsour=[calibrator], sources=[target],
                refant=refantlist[0])

def evn_plot2a(uvdata):
    uvdata.zap_table('AIPS PL', -1)

    # plot the bandpass table
    if nbp_table:
        aparm = set_default_aparms('possm')
        aparm[4] = 1.3
        aparm[5:7] = [0, 0]
        aparm[8] = 2
        #print >>PYPELOG, 'aparm=', aparm
        runpossm(indata=uvdata, aparm=aparm, stokes=plotstokes, solint=0)
        plot(uvdata, def_name('BANDPASS'), dopng=dopng)

    # plot the fringe solutions
    uvdata.zap_table('AIPS PL', -1)
    runsnplt(uvdata, 'CL', 3, 'PHAS')
    plot(uvdata, def_name('FRING_PHAS'), dopng=dopng)

    uvdata.zap_table('AIPS PL', -1)
    runsnplt(uvdata, 'SN', 2, 'DELA')
    plot(uvdata, def_name('FRING_DELAY'), dopng=dopng)

    uvdata.zap_table('AIPS PL', -1)
    runsnplt(uvdata, 'SN', 2, 'RATE')
    plot(uvdata, def_name('FRING_RATE'), dopng=dopng)


def evn_plot2b(uvdata):
    uvdata.zap_table('AIPS PL', -1)

    # plot the calibrated data versus frequency
    aparm = set_default_aparms('possm')
    aparm[1] = 1        # vector average 
    aparm[5] = 0
    aparm[6] = 0
    uvdata.zap_table('AIPS PL', -1)
    runpossm(indata=uvdata, aparm=aparm, docalib=1, doband=nbp_table,
            antennas=plotref, stokes=plotstokes)
    plot(uvdata, def_name('POSSM_CAL'), dopng=dopng)

    # plot the calibrated data versus time
    bparm = set_default_bparms('vplot')
    uvdata.zap_table('AIPS PL', -1)
    for poln in polarizations:
        for incif in range(1, nif+1):
            error = runvplot(indata=uvdata, bparm=bparm, antenna=plotref,
                    bif=incif, eif=incif, docalib=1, doband=nbp_table,
                    stokes=poln+poln, solint=plotavg)
            if (error):
                print_vplot_err(incif, incif, poln+poln, uvdata)
    plot(uvdata, def_name('VPLOT_CAL'), dopng=dopng)


def evn_multi(uvdata, sources, target_sources):
    # create multi files and make dirty maps and first clean maps
    for source in sources:
        print >>PYPELOG, '\nmapping source=', source
        splitdata = AIPSUVData(source, 'SPLIT', uvdata.disk, 1)
        try:
            splitdata.verify()
            assert (len(splitdata) > 0)
        except:
            print >>PYPELOG, '\nPype warning: There are no SPLIT data for ', source
            continue
        else:
            # Produce a multi file and map it
            multidata = AIPSUVData(splitdata.name, 'MULTI', uvdata.disk, 1)
            zap_old_data(multidata)
            runmulti(splitdata, multidata)
            runindxr(multidata)

            # run uvprm to fill in some useful keywords. Use the highest
            # available IF (longest u,v).
            for bif in range(nif, 0, -1):
                try:
                    runuvprm(multidata, nchan//2, bif=bif)
                    break
                except:
                    print >>PYPELOG, '\nPype warning: uvprm failed for IF= ', bif
                    continue
            #cellsize = getcell(multidata)

            try:
                if source in target_sources:
                    print >>PYPELOG, 'Uniformly weighted dirty map'
                    imgdata = AIPSImage(multidata.name, 'IIM001', uvdata.disk, 1)
                    plotname = (source + '_IMAPU' )
                    print >>PYPELOG, 'Creating image:', aipsuvname(imgdata)
                    create_image(multidata, imgdata, 0, 'U')
                    plot(imgdata, def_name(plotname), dopng=dopng)
                    save_fits_image(imgdata, plotname)

                    print >>PYPELOG, 'Naturally weighted dirty map'
                    imgdata = AIPSImage(multidata.name, 'IIM001', uvdata.disk, 2)
                    plotname = (source + '_IMAPN' )
                    print >>PYPELOG, 'Creating image:', aipsuvname(imgdata)
                    create_image(multidata, imgdata, 0, 'N')
                    plot(imgdata, def_name(plotname), dopng=dopng)
                    save_fits_image(imgdata, plotname)

                print >>PYPELOG, 'Uniformly weighted clean map'
                plotname = (source + '_ICLN_' + '1')
                imgdata = AIPSImage(multidata.name, 'ICL001', uvdata.disk, 1)
                print >>PYPELOG, 'Creating image:', aipsuvname(imgdata)
                create_image(multidata, imgdata, 100, 'U')
                plot(imgdata, def_name(plotname), dopng=dopng)
            except:
                print >>PYPELOG, '\nPype warning: First round of mapping has failed for ', source

def evn_map():
    for source in selfcal_sources:
        print >>PYPELOG, '\nmapping source=', source
        uvdata = AIPSUVData(source, 'MULTI', disk, 1)
        try:

            uvdata.verify()

            # delete old CL/SN tables
            uvdata.zap_table('AIPS SN', -1)
            while uvdata.table_highver('AIPS CL') > 1:
                uvdata.zap_table('AIPS CL', 0)

            nsc = nselfcal 

            ## make sure at least one source gets 'a&p' selfcal
            #if source == (selfcal_sources)[0] and nsc < 2:
            #    nsc = 2 

            # iterate on selfcal and mapping
            selfcal_map(uvdata, source, nsc, solint, refantlist[0],
                    prefix=output_prefix, dopng=dopng)
        except:
            print >>PYPELOG, '\nPype warning: Mapping has failed for ', source
            traceback.print_exc(file=PYPELOG)
            continue


def evn_plot3():
    for source in sources:
        print >>PYPELOG, '\nPlotting final calibrated data of source: ', source
        if source in target_sources:
            imageseq = 1
        else:
            imageseq = nselfcal+1

        uvdata = AIPSUVData(source, 'MULTI', disk, 1)
        imgdata = AIPSImage(source, 'ICL001', uvdata.disk, imageseq)
        try:
            uvdata.verify()
            imgdata.verify()
        except:
            print >>PYPELOG, '\nPype warning: There is no MULTI or ICLN file for ', source
            continue
        else:
            # uvplt the calibrated vis. versus u,v distance
            print >>PYPELOG, 'Plot amplitude versus u,v distance'
            uvdata.zap_table('AIPS PL', -1)
            bparm = set_default_bparms('uvplt')
            bparm[8] = 1000
            runuvplt(indata=uvdata, bparm=bparm, sources=[source])
            plot(uvdata, def_name(source + '_UVPLT'), dopng=dopng)

            # uvplt the u,v coverage
            print >>PYPELOG, 'Plot u,v coverage'
            uvdata.zap_table('AIPS PL', -1)
            bparm = set_default_bparms('uvplt')
            bparm[1:] = [6, 7, 2, 0]
            runuvplt(indata=uvdata, bparm=bparm, sources=[source])
            plot(uvdata, def_name(source + '_UVCOV'), dopng=dopng)

            if not (source in target_sources):
                # plot the model and the data together
                print >>PYPELOG, 'Plot model and data together'
                uvdata.zap_table('AIPS PL', -1)
                bparm = set_default_bparms('vplot')
                for poln in polarizations:
                    for incif in range(1, nif+1):
                        error = runvplot(indata=uvdata, bparm=bparm,
                                antenna=plotref, bif=incif, eif=incif,
                                docalib=1, in2data=imgdata, nmaps=1,
                                sources=[source], stokes=poln+poln,
                                solint=plotavg)
                        if (error):
                            print_vplot_err(incif, incif, poln+poln, uvdata)
                plot(uvdata, def_name(source + '_VPLOT_MODEL'), dopng=dopng)

                # plot the closure phase with the model subtracted (needs
                # single source file)
                try:
                    print >>PYPELOG, 'Plot closure phases of data-model'
                    uvdata.zap_table('AIPS PL', -1)
                    uvsubdata = plot_closure(uvdata, imgdata, source, nselfcal)
                    plot(uvsubdata, def_name(uvdata.name + '_CLPHS'), dopng=dopng)
                except:
                    print >>PYPELOG, '\nPype warning: closure phase plotting' \
                            ' failed for', source, '\n'
                    continue



def evn_sens(uvdata, selfcal_sources):
    # calculate the antenna sensitivities

    # copy tables from the pre-split multi source file
    tasavdata = AIPSUVData(uvdata.name, 'TMP', uvdata.disk, 1)
    zap_old_data(tasavdata)
    runtasav(indata=uvdata, outdata=tasavdata)
    # delete the old sn tables
    tasavdata.zap_table('AIPS SN', -1)
    assert (tasavdata.table_highver('AIPS CL') == 3)

    # copy the A&P solutions from each of the single source files.
    for source in selfcal_sources:
        multidata = AIPSUVData(source, 'MULTI', uvdata.disk, 1)
        try:
            multidata.verify()
        except:
            print >>PYPELOG, '\nPype warning: There is no MULTI file for ', source
            continue
        else:
            runtacop(indata=multidata, outdata=tasavdata, inext='SN', inver=2, 
                outver=0, ncount=1)

    # apply all SN tables from A&P selfcal to the CL table with the ANTAB
    # GAINS
    runclcal(indata=tasavdata, opcode='CALI', interpol='2PT', snver=0, 
            gainver=3, gainuse=4, refant=refantlist[0])

    # plot the resulting cl table
    tasavdata.zap_table('AIPS PL', -1)
    runsnplt(tasavdata, 'CL', 4, 'AMP')
    plot(tasavdata, def_name('SENS'), dopng=dopng)

def evn_save(uvdata):
    # save various useful products
    tasavdata = AIPSUVData(uvdata.name, 'TASAV', uvdata.disk, 1)
    zap_old_data(tasavdata)
    runtasav(indata=uvdata, outdata=tasavdata)

    # save a copy of the tables
    fitsoutfile = output_prefix + '.tasav.FITS'
    save_old_file(fitsoutfile)
    runfittp(indata=tasavdata, outfile=fitsoutfile)

    # print the history 
    hisfile = output_prefix + '.his'
    HISFILE = open (hisfile, 'w')
    i = 0
    for line in uvdata.history:
        i += 1
        print >>HISFILE, str(i), ' ', line

    # save fits files
    for source in selfcal_sources + target_sources:
        print >>PYPELOG, '\nsaving fits and ps images for source=', source
        # save ps and fits versions of the final images
        sequence = nselfcal+1
        if source in target_sources:
            sequence = 1
        imgdata = AIPSImage(source, 'ICL001', uvdata.disk, sequence)
        try:
            imgdata.verify()
        except:
            print >>PYPELOG, '\nPype warning: There is no image file, ', aipsuvname(imgdata), \
                    ', for ', source
            continue
        else:
            #imgdata.zap_table('AIPS PL', -1)
            #plotmap(imgdata)
            plot(imgdata, def_name(source + '_ICLN'), dopng=dopng)
            fitsoutfile = output_prefix + '_' + source + \
                        '_ICLN.FITS'
            save_old_file(fitsoutfile)        
            runfittp(indata=imgdata, outfile=fitsoutfile)

        # save the split data as fits
        splitdata = AIPSUVData(source, 'SPLIT', uvdata.disk, 1)
        try:
            splitdata.verify()
        except:
            print >>PYPELOG, '\nPype warning: There is no split file, ', \
                    aipsuvname(splitdata), ', for ', source
            continue
        else:
            fitsoutfile = output_prefix + '_' + source + \
                        '.UVDATA.FITS'
            save_old_file(fitsoutfile)        
            runfittp(indata=splitdata, outfile=fitsoutfile)



#AIPS.userno = 3601
#AIPS.disks = [ None, AIPSDisk(None, 1), AIPSDisk('http://jop31:8000', 1) ]
#AIPS.proxies = [ LocalProxy, ServerProxy('http://jop31:8000', 1) ]

######################
# Program entry point
#

# the current version of the code
# Note that $Revision: 270 $ is a SVN keyword
# Note that $Date: 2007-06-19 19:25:01 +0200 (Tue, 19 Jun 2007) $ is a SVN keyword
def mod_svntag(svntag, keyword):
    svntag = svntag.replace(keyword+':', '')
    svntag = svntag.replace(r'$', '')
    return svntag

svnurl = '$HeadURL: file:///export/jive/reynolds/svnroot/repos/Pypeline/tags/PYPELINE-4-3/Pypeline/EVN.py $'
svnurl = mod_svntag(svnurl, 'HeadURL')
svnrev = '$Revision: 270 $';
svnrev = mod_svntag(svnrev, 'Revision')
svndate = '$Date: 2007-06-19 19:25:01 +0200 (Tue, 19 Jun 2007) $';
svndate = mod_svntag(svndate, 'Date')
parenth = re.compile(r'\(.*\)')
svndate = parenth.sub(r'', svndate)

usage = 'usage: %prog [options] experiment.inp'
parser = optparse.OptionParser(usage=usage, version='%prog ' + svnrev + svndate)
(options, args) = parser.parse_args()
if len(args) != 1:
    parser.error("incorrect number of arguments")
today = time.asctime()




# some necessary global variables
#disk = 1
#niter = 300

# get the inputs and print them
#control = dict()
#control = parse_inp(args[0])
control = parse_inp(args[0])

# check the inputs and re-type where necessary
checkin(control)

vlog_out = output_prefix + '.vlba'
pypelog = output_prefix + '.pypelog'

# delete the log file if we are starting again. 
if tmask[0] <= 1 <= tmask[1]:
    if os.path.exists(pypelog):
        os.remove(pypelog)



AIPSTask.msgkill = msgkill

# redirect important messages to both pypelog and stdout
PYPELOG = crfuncs.Tee('a', pypelog, sys.stdout)
sys.stderr = PYPELOG
# messages from evn_funcs should also go to pypelog
evn_funcs.FuncLog = PYPELOG
evn_aips_tasks.FuncLog = PYPELOG

# redirect AIPS messages to AIPS log file, and filtered selection to pypelog
aips_startmsg = crfuncs.Msg_Filter(PYPELOG, 'begins', 'ended', 'died',
        'warning')
AIPS.log = crfuncs.Tee('w', aips_startmsg, outdir + '/' + 'AIPS.LOG')
#AIPS.log = open (os.getcwd() + '/' + 'AIPS.LOG', 'w')


# print some useful information to pypelog
print >>PYPELOG, '*' * 20
print >>PYPELOG, 'Pipeline started on ', today
print >>PYPELOG, 'Using svn revision: ', svnrev, ' ', svndate, ' ', svnurl 

print >>PYPELOG, "\ninputs:\n"
print >>PYPELOG, pprint.pformat(control)

print >>PYPELOG, 'doing tmask ', tmask[0], ' to ', tmask[1]

# get some necessary files assuming standard names
uvflg_file = getfile('uvflg')
assert(os.path.isfile(uvflg_file)), uvflg_file + ' does not exist!'
chflg_file = getfile('chflag')
antab_file = getfile('antab')
assert(os.path.isfile(antab_file)), antab_file + ' does not exist!'
vlbacal_file = getfile('vlbacal')

# Now we start the main part of the program. Most parts are wrapped in blocks
# so that they are only executed if tmask is appropriate. Need to be careful to
# set the program state correctly if parts are skipped (get right uvdata, sn
# tables, etc.).

# Any or all of the following datasets may be created by the loading procedure,
# depending on the contents of the dataset.
uvname = experiment.upper()
uvcopklass = 'FQ' + str(freqid)
uvdata = AIPSUVData(uvname, 'UVDATA', disk, 1)
msortdata = AIPSUVData(uvdata.name, 'MSORT', uvdata.disk, 1)
vbgludata = AIPSUVData(uvdata.name, 'VBGLU', uvdata.disk, 1)
fxpoldata = AIPSUVData(uvdata.name, 'FXPOL', uvdata.disk, 1)
uvcopdata = AIPSUVData(uvdata.name, uvcopklass, uvdata.disk, 1)
uvavgdata = AIPSUVData(uvdata.name, 'UVAVG', uvdata.disk, 1)

# 1, Load and sort the data - average too if requested
if tmask[0] <= 1 <= tmask[1]:
    print_start('1, load and sort the data')
    evn_load(uvdata, msortdata, vbgludata, fxpoldata, uvcopdata, uvavgdata)
    print_end('1')

# Always do this bit to make sure names are right at end of block. The order of
# the if statements here is significant.
# This is slightly awkward, but can't think of a better way...
if msortdata.exists():
    uvdata = msortdata
if vbgludata.exists():
    uvdata = vbgludata
if fxpoldata.exists():
    uvdata = fxpoldata
if uvcopdata.exists():
    uvdata = uvcopdata
if uvavgdata.exists():
    uvdata = uvavgdata

print 'uvdata=', uvdata
# check the data and extract useful information (freq, shape, etc.). Initialise
# some things as well
sources = []
selfcal_sources = []
plotref = []
nchan = None
nif = None
chanwidth = None
ifwidth = None
plotstokes = None


(refantlist, sources, selfcal_sources, plotref, ifwidth, chanwidth,
        nchan, nif, stokes, plotstokes, polarizations) = evn_checkdata()

# get the source list and scan details from the NX table
sources = read_nx(uvdata, sources)

# only print the source/scan info if reloading the data
if tmask[0] <= 1 <= tmask[1]:
    print >>PYPELOG, '\n', '*' * 20
    #print >>PYPELOG, "sources=\n", pprint.pformat(sources)

    # Make a copy of the sources dictionary so we can truncate the scanlen for
    # printing purposes.
    source_print = copy.deepcopy(sources)
    # truncate the 'scanlen' numbers to a sensible number of decimal places
    for source in source_print:
        for i in range(len(source_print[source]['scanlen'])):
            source_print[source]['scanlen'][i] = "%.3f" % source_print[source]['scanlen'][i] 

    print >>PYPELOG, "sources=\n", pprint.pformat(source_print)
    print >>PYPELOG, 'stokes=', stokes
    print >>PYPELOG, 'polarizations=', polarizations


if (not solint):
    # solint defaults to the typical scan length on the phase reference
    solint = getsolint(sources, phaseref_sources)

print >>PYPELOG, "solint=", solint, ' min'

# 2, A priori data flagging 
if tmask[0] <= 2 <= tmask[1]:
    print_start('2, flag using uvflg monitor data and flag band edges and run vlog')
    table_vers(uvdata=uvdata, cl=1, sn=0, fg=0, bp=0)

    evn_flag(uvdata)
    print_end('2')

# 3, plot the raw data
if tmask[0] <= 3 <= tmask[1]:
    print_start('3, plot the data - vs time and frequency')
    table_vers(uvdata=uvdata, cl=1, sn=0, fg=1, bp=0)

    evn_plot1()
    print_end('3')


# 4, Amplitude calibration and parallactic angle correction
if tmask[0] <= 4 <= tmask[1]:
    print_start('4, amplitude calibration and parallactic angle correction')
    table_vers(uvdata=uvdata, cl=1, sn=0, fg=1, bp=0)

    evn_ampcal(antab_file, uvdata)
    print_end('4')


# 5, Fringe fitting.
if tmask[0] <= 5 <= tmask[1]:
    print_start('5, fringe fitting')
    table_vers(uvdata=uvdata, cl=2, sn=1, fg=1, bp=0)

    evn_fring(uvdata)
    print_end('5')


# 6, Bandpass calibration
if tmask[0] <= 6 <= tmask[1] and nbp_table:
    print_start('6, bandpass calibration')
    table_vers(uvdata=uvdata, cl=3, sn=2, fg=1, bp=0)

    runbpass(refant=refantlist[0],indata=uvdata, calsour=bpass_calibrators)
    # Do the less time consuming plots now
    evn_plot2a(uvdata)
    print_end('6')


# 7, Plot the results after ampcal, fringe fitting and bandpass
if tmask[0] <= 7 <= tmask[1] and doplot:
    print_start('7, plot calibrated data vs time and frequency' )
    table_vers(uvdata=uvdata, cl=3, sn=2, fg=1, bp=nbp_table)

    evn_plot2b(uvdata)
    print_end('7')

# 8, Split
if tmask[0] <= 8 <= tmask[1]:
    print_start('8, split the calibrated data')
    table_vers(uvdata=uvdata, cl=3, sn=2, fg=1, bp=nbp_table)

    for source in sources:
        splitdata = AIPSUVData(source.upper(), 'SPLIT', uvdata.disk, 1)
        zap_old_data(splitdata)

    runsplit(sources=sources.keys(), indata=uvdata, gainuse=3,
            doband=nbp_table, bpver=nbp_table, outseq=1)

    print_end('8')

# 9, create multi files and make dirty maps and first clean maps
if tmask[0] <= 9 <= tmask[1]:
    print_start('9, create multi files and make first maps')
    evn_multi(uvdata, sources, target_sources)

    print_end('9')


# 10, continue mapping
if tmask[0] <= 10 <= tmask[1]:
    print_start('10, iterate self-cal and imaging')
    evn_map()
    print_end('10')

# 11, plot the final data
if tmask[0] <= 11 <= tmask[1] and doplot:
    print_start('11, plot the final self-calibrated data')
    evn_plot3()
    print_end('11')

# 12, calculate the antenna sensitivities
if tmask[0] <= 12 <= tmask[1]:
    print_start('''12, calculate the antenna sensitivities (using selfcal
            results)''')
    if nselfcal > 1:
        evn_sens(uvdata, selfcal_sources)
    else:
        print >>PYPELOG, '''No amplitude self-cal was done, so no sensitivity
            plots can be made (skipping tmask 12)'''

    print_end('12')

# 13, save useful data and plot final map
if tmask[0] <= 13 <= tmask[1]:
    print_start('13, save useful data and plot final map')
    table_vers(uvdata=uvdata, cl=3, sn=2, fg=1, bp=nbp_table)

    evn_save(uvdata)
    print_end('13')
    print >>PYPELOG, '*' * 20
    print >>PYPELOG, "Pipeline has completed and data are saved - now archive it"
    print >>PYPELOG, '*' * 20


# now wait for the pids of the child plot conversion processes to finish
for pid in evn_funcs.pid_list:
    print >>PYPELOG, '\n\nWaiting for the child process ( pid = ', pid, ') to'\
            ' finish...\n'
    os.waitpid(pid, 0)

#print 'sources =\n', pprint.pformat(sources)

print >>PYPELOG, '*' * 20
print >>PYPELOG, "The pipeline has ended normally at ", time.asctime()
print >>PYPELOG, '*' * 20
print >>PYPELOG, "\n\n"