import os, sys
from AIPS import AIPS
from AIPSTask import AIPSTask as task
from AIPSData import AIPSUVData as UV
from AIPSData import AIPSImage as IM
from Wizardry.AIPSData import AIPSUVData as WUV
import time
import numpy as np

# Used for amp selfcal etc. We trust the VLBA.
vlba_antennas = [None, 7, 8, 9, 10, 11, 12, 13, 14, 15, 18]
evn_antennas = [None, 1, 2, 3, 4, 5]
# AR and YY need special attention

#Ant   1 = EF       BX= -3033407.4019 BY=  4058261.8171 BZ=  1115208.4638
#Ant   2 = ON       BX= -3367090.1054 BY=  3442412.9284 BZ=  1564608.3968
#Ant   3 = MC       BX= -3387647.9432 BY=  4552581.3919 BZ=   664336.8632
#Ant   4 = NT       BX= -3703604.6301 BY=  5086781.0174 BZ=    21262.1975
#Ant   5 = WB       BX= -3026932.3864 BY=  3848658.5180 BZ=  1279699.2271
#Ant   6 = GB       BX=  1768773.8310 BY=    38513.9638 BZ=   158506.9185
#Ant   7 = SC       BX=  2615142.7749 BY=  1643932.8542 BZ= -1852482.7995
#Ant   8 = HN       BX=  1393863.7199 BY=   674722.4385 BZ=   537083.6974
#Ant   9 = NL       BX=  1437460.9665 BY=  -932959.6585 BZ=   441628.5668
#Ant  10 = FD       BX=  1797723.5130 BY= -2205175.2118 BZ=  -553260.0247
#Ant  11 = LA       BX=  1424734.6372 BY= -2268865.2869 BZ=   -76098.5837
#Ant  12 = PT       BX=  1431406.6479 BY= -2463993.6381 BZ=  -209810.6327
#Ant  13 = KP       BX=  1393701.3387 BY= -2817425.9393 BZ=  -427894.3961
#Ant  14 = OV       BX=   773174.8314 BY= -3130636.7628 BZ=    53394.8842
#Ant  15 = BR       BX=    61209.3592 BY= -2707282.5699 BZ=   941591.2525
#Ant  16 = YY       BX=  1464891.5897 BY= -2429381.2214 BZ=  -230346.8089
#Ant  17 = AR       BX=  2654009.6505 BY=  1416749.0246 BZ= -1790501.9279
#Ant  18 = MK       BX= -1697410.3928 BY= -5806893.6653 BZ= -1636925.2948

whattodo = {'load_data': False, # ->CL1
            'load_tables': False,
            'msort': False,
            'indxr': False,
            'preflagging': False,
            'add_YY_to_table': False,
            'copy_JIVE_CL': False, #-> CL2. 
            # This Cl Contains A priori AMP, Parr. angle. 
            'fring_ampcal' :False, #->SN1, CL3
            'bpass_ampcal': False, # BP 1
            'image_ampcal_1': False, 
            # A priori ampcal is wrong for antennas 16,17. Adjust with selfcal.
            'selfcal_ampcal' : False, # -> SN2, 3. CL4
            'clip_ampcal': False, #USE PREFLAGGING AS WELL TO CLEAN FG
            'image_ampcal_2': False, 
            'fring_phasecal' : False, #->SN4, CL5
            'clip_phasecal': False, #USE PREFLAGGING AS WELL TO CLEAN FG
            # Also clip target at the same time
            'image_phasecal_1':False,
            'selfcal_phasecal' : False, # -> SN5, 6, CL6
            'image_phasecal_2': False, # 
            'fring_frcal' : False, #->SN7, CL7
            'image_frcal_1': False, 
            'selfcal_frcal' : False, # -> SN8, CL8
            'image_frcal_2': False, 
            'image_ARP220': False, # 
            'selfcal_arp220': False, # ->SN9, CL9
            'image_ARP220_final': True,
            'export_results': True,
            'save_calsols': True,
            }
 
########## Initialize observation data ##########
AIPS.userno = 1000
clint = 15.0/60 # Seconds, same as JIVE table
NAME = 'GC028_X'
TABNAME = NAME
CLASS = 'UVDATA'
TABCLASS = 'TAB'
DISK = 1
AMPCALIM = NAME + '_AC'
PCALIM = NAME + '_PC'
TCALIM = NAME + '_TC'
TCALBLANKIM = NAME + 'TCBL'
TCALIM2 = TCALIM + '2'
TCALBLANKIM2 = TCALBLANKIM + '2'
FRCALIM = NAME + '_FC'
OUTPREFIX = NAME
target = 'ARP220'
ampcal = 'J1516+1932'
phasecal = 'J1532+2344'
frcal = 'J1613+3412'
sorted_class = 'MSORT'
logfile = 'log_'+OUTPREFIX+'_PTred.log'
refant = 12
#########^^^ END OF CONFIG VARIABLES ^^^################
tic = time.time()
########## Initialize Log file ##########
try:
 os.system('cp ' + logfile + ' ' + logfile + '.old')
 os.system('rm -f ' + logfile)
except:
 pass
AIPS.log = open(logfile, 'a')
AIPS.log.write('whattodo = '+repr(whattodo)+'\n')

########## LOAD DATA ##########
if whattodo['load_data']:
    fitld = task('fitld')
    fitld.default()
    outdata = UV(NAME, CLASS, DISK, 1)
    fitld.datain = '/data1/eskil/RAWDATA/GC028A/gc028_1_1.IDI1'
    fitld.outdata = outdata
    fitld.ncount = 2
    fitld.digicor = -1 # Not for JIVE
    fitld.clint = clint
    fitld.wtthresh = 0.6 # From PI Letter
    fitld.doconcat = 1
    fitld.douvcomp = -1
    if outdata.exists():
        print 'WARNING: Removing existing file ' + NAME + '.' + CLASS
        outdata.zap()
        print '         File removed. Proceeding with FITLD...'
    fitld.go()

######### LOAD PIPELINE CAL TABLES ########
if whattodo['load_tables']:
    fitld = task('fitld')
    fitld.default()
    outdata = UV(TABNAME, TABCLASS, DISK, 1)
    fitld.datain = '/data1/eskil/RAWDATA/GC028A/gc028.tasav.FITS'
    fitld.outdata = outdata
    fitld.ncount = 1
    if outdata.exists():
        print 'WARNING: Removing existing file ' + TABNAME + '.' + TABCLASS
        outdata.zap()
        print '         File removed. Proceeding with FITLD...'
    fitld.go()

if whattodo['msort']:
    msort = task('msort')
    data = UV(NAME, CLASS, DISK, 1)
    outdata = UV(NAME, sorted_class, DISK, 1)
    msort.default()
    msort.indata = data
    msort.outdata = outdata
    msort.sort = 'TB'
    if outdata.exists():
        print 'WARNING: Removing existing file ' + NAME + '.' + sorted_class
        outdata.zap()
        print '         File removed. Proceeding with FITLD...'
    msort.go()
    if data.exists():
        print 'WARNING: Removing existing file ' + NAME + '.' + CLASS
        data.zap()
        print '         File removed. Proceeding with FITLD...'

if whattodo['indxr']:
    indxr = task('indxr')
    indxr.default()
    data = UV(NAME, sorted_class, DISK, 1)
    indxr.indata = data
    indxr.cparm[3] = clint
    indxr.go()

########## PREFLAGGING ##########
#Flag data known to be bad (whole antennas etc.)
if whattodo['preflagging']:
    # Copy flag table provided from JIVE
    # This will flag when antennas off source etc.
    indata = UV(TABNAME, TABCLASS, DISK, 1)
    outdata = UV(NAME, sorted_class, DISK, 1)
    # Remove all FG tables
    outdata.zap_table('FG', -1)
    tacop = task('tacop')
    tacop.default()
    tacop.indata = indata
    tacop.outdata = outdata
    tacop.inext = 'FG'
    tacop.inver = 1
    tacop.outver = 1
    tacop.go()
    data = UV(NAME, sorted_class, DISK, 1)
    infgver = 1
    outfgver = 2
    # Remove all FG tables higher than infgver.
    for i in range(data.table_highver('FG'), infgver, -1):
     data.zap_table('FG', i)
    data = UV(NAME, sorted_class, DISK, 1)
    # Copy flag table 1 to flagtable 2
    tacop = task('tacop')
    tacop.default()
    tacop.indata = data
    tacop.inext = 'FG'
    tacop.inver = infgver
    tacop.ncount = 1
    tacop.outdata = data
    tacop.outvers = outfgver
    tacop.go()
    # REMARKS ON BAD ANTENNAS FROM PI LETTER:
#  Ef: station reported shadowing by surrounding hills at low elevations
#    near beginning and end of experiment; subreflector control was hung
#    between UT 07:25 and 08:03, pointing may be affected in a few earlier
#    scans, but OK after 08:03; lost UT 11:30-11:38 due to break in
#    communication with the Mark5 unit.
#  
#    On: OK. Has RCP only (RCP recorded in all channels). Data in "LCP"
#    channels were flagged before producing the FITS files. Station
#    reported heavy fog. (NB standard plots of amplitude/phase vs time
#    show LCP data only, so for On and Nt these show cross-pols, plotted before
#    flagging)
#  
#    Nt: OK. Has RCP only (same as Onsala).
#  
   #Flag antenan 4 NT and 2 Onsala LCP since not present (PiLetter)
    data = UV(NAME, sorted_class, DISK, 1)
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    # FLAG BOTH POLS SINCE THESE ANTENNAS ONLY ADD NOISE AND RIPPLES
    #uvflg.stokes = 'LL'
    uvflg.antennas = [None, 2,4]
    uvflg.reason = 'NO LCP, PIletter'
    uvflg.opcode = 'FLAG'
    uvflg.go()
   #Flag timeran at antenna 2 Onsala which seems bad (microseconds spikes in delay etc)
    data = UV(NAME, sorted_class, DISK, 1)
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.antennas = [None, 2]
    #uvflg.timeran = [None, 0, 16, 0, 0, 0, 18, 0, 0] # bad delay
    uvflg.timeran = [None, 0, 0, 0, 0, 0, 18, 0, 0]
    uvflg.reason = 'BadD+A'
    uvflg.opcode = 'FLAG'
    uvflg.go()
    
#  
#    Mc: no good LCP data in subband 3 (AIPS IF 4) before UT 15:28, due to
#    VC #8 being wrongly connected on the patch panel. Problem with phase
#    stability at X-band (noted also in the Network Monitoring Experiment
#    from the same session).
    data = UV(NAME, sorted_class, DISK, 1)
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.antennas = [None, 3]
    # Do not use the timerange, flag all times for these IFs 
    # LL pol since lrge rates at the small end time remaining.
    #uvflg.timeran = [None, 0, 0, 0, 0, 0, 15, 28, 0]
    #uvflg.stokes = 'LL'
    #uvflg.bif = 3
    #uvflg.eif = 4
    uvflg.reason = 'PIletter'
    uvflg.opcode = 'FLAG'
    uvflg.go()
#    Wb: used in phased array mode. Station reported poor phase stability
#    on long baselines.
#  
#    Gb: no fringes found despite fairly exhaustive searching as outlined
#    in email from Bob Campbell. Correlated with nominal clocks.
   #Flag GB, antenna 6, since no fringes found (also mentioned in PI-letter)
    data = UV(NAME, sorted_class, DISK, 1)
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.antennas = [None, 6]
    uvflg.reason = 'NO fringes'
    uvflg.opcode = 'FLAG'
    uvflg.go()
#  
#    Fd: seems to have a possible problem with low sensitivity and high
#    polarization leakage. In cross-correlations, RCP amplitude in subband
#    1 (AIPS IF 2) is low relative to LCP, while in subbands 2 and 3 (AIPS
#    IF 3 and 4), LCP amplitude is lower than RCP.
    data = UV(NAME, sorted_class, DISK, 1)
   #Flag IF 1 at antenna 10 FD since it seems bad (microseconds spikes in delay etc)
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.antennas = [None, 10]
    uvflg.bif = 1
    uvflg.eif = 1 
    uvflg.reason = 'PIletter, no fringsols'
    uvflg.opcode = 'FLAG'
    uvflg.go()
   # Flag timeran, IF 4 antenna 10 FD since almost no FRING solutions found 
   # Also see PI letter notes
    data = UV(NAME, sorted_class, DISK, 1)
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.antennas = [None, 10]
    uvflg.timeran = [None, 0, 16, 0, 0, 1, 2, 0, 0]
    uvflg.bif = 4
    uvflg.eif = 4 
    uvflg.reason = 'PIletter'
    uvflg.opcode = 'FLAG'
    uvflg.go()

    # FD HAS MANY PROBLEMS (large variations in amp sols etc;. FLAG IT ALL!
    data = UV(NAME, sorted_class, DISK, 1)
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.antennas = [None, 10]
    uvflg.reason = 'Problems'
    uvflg.opcode = 'FLAG'
    uvflg.go()
#  
#    Ov: OK (below horizon for first scan shown in standardplots)
#  
#    Y: very low amplitude for LCP in subband 1 (AIPS IF 2, corresponds to
#    BBC #4). Also low amplitude for RCP in subband 3 (AIPS IF 4,
#    corresponding to BBC #7) - raw autocorrelations look bad (close to
#    0). After the 2-bit van Vleck correction, there is a large spike in
#    the subband 3 RCP autocorrelations.
#  
#    Ar: used an incorrect schedule, ~40s out of synch with other
#    stations. Hence the last 41-43s of every scan are missed.
    data = UV(NAME, sorted_class, DISK, 1)
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.antennas = [None, 5]
    uvflg.bif = 2
    uvflg.eif = 2
    uvflg.stokes = 'LL'
    uvflg.timeran = [None, 0, 7, 56, 0, 0, 7, 58, 0]
    uvflg.reason = 'bad delay'
    uvflg.opcode = 'FLAG'
    uvflg.go()

   
   #Flag all data for elevations less than 10 deg
    data = UV(NAME, sorted_class, DISK, 1)
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.aparm = [None, -90, 10]
    uvflg.reason = 'el.lt.10deg'
    uvflg.opcode = 'FLAG'
    #uvflg.go()
   
    data = UV(NAME, sorted_class, DISK, 1)
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.antenna = [None, 16]
    uvflg.bif = 2
    uvflg.eif = 2
    uvflg.stokes = 'LL'
    uvflg.bchan = 1
    uvflg.echan = 8
    uvflg.opcode = 'FLAG'
    uvflg.go()
    
    data = UV(NAME, sorted_class, DISK, 1)
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.antenna = [None, 1]
    uvflg.timeran = [None, 0, 7, 56, 0, 0, 7, 58, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    
    data = UV(NAME, sorted_class, DISK, 1)
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.antenna = [None, 1]
    uvflg.timeran = [None, 0, 16, 50, 0, 0, 17, 0, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    
    data = UV(NAME, sorted_class, DISK, 1)
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    #uvflg.antenna = [None, 1]
    #uvflg.baselin = [None, 14, 16]
    uvflg.timeran = [None, 0, 12, 0, 0, 0, 12, 15, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    
    data = UV(NAME, sorted_class, DISK, 1)
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.antenna = [None, 4]
    uvflg.timeran = [None, 0, 10, 23, 0, 0, 10, 27, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    
    data = UV(NAME, sorted_class, DISK, 1)
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.antenna = [None, 16]
    uvflg.timeran = [None, 0, 20, 30, 0, 0, 21, 30, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    
    data = UV(NAME, sorted_class, DISK, 1)
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.antenna = [None, 17]
    uvflg.timeran = [None, 0, 14, 12, 0, 0, 14, 16, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
   
    
    data = UV(NAME, sorted_class, DISK, 1)
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.antenna = [None, 10]
    uvflg.timeran = [None, 0, 12, 42, 0, 0, 12, 45, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    
    data = UV(NAME, sorted_class, DISK, 1)
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.timeran = [None, 0, 11, 31, 0, 0, 11, 32, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    
    data = UV(NAME, sorted_class, DISK, 1)
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.antenna = [None, 16]
    uvflg.timeran = [None, 0, 23, 7, 0, 0, 23, 9, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
   

    data = UV(NAME, sorted_class, DISK, 1)
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.antenna = [None, 10]
    uvflg.bif = 4
    uvflg.eif = 4
    uvflg.opcode = 'FLAG'
    uvflg.go()
    
    # Quack
    quack = task('quack')
    quack.default()
    quack.indata = data
    quack.opcode = 'ENDB'
    quack.antenna = [None, 1]
    quack.baselin = [None, 5]
    quack.aparm[2] = 10.0/60
    quack.go()
    
     # noisy baseline. On is noisy, but this is the worst
    data = UV(NAME, sorted_class, DISK, 1)
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.antenna = [None, 2]
    uvflg.baselin = [None, 4]
    uvflg.opcode = 'FLAG'
    uvflg.go()

    data = UV(NAME, sorted_class, DISK, 1)
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.antenna = [None, 9]
    uvflg.timeran = [None, 1, 0, 0, 0, 1, 30, 0, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    
    data = UV(NAME, sorted_class, DISK, 1)
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.antenna = [None, 7]
    uvflg.timeran = [None, 0, 21, 15, 0, 0, 21, 20, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    data = UV(NAME, sorted_class, DISK, 1)
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.antenna = [None, 7]
    uvflg.timeran = [None, 0, 21, 35, 0, 0, 21, 40, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    data = UV(NAME, sorted_class, DISK, 1)
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.antenna = [None, 8]
    uvflg.timeran = [None, 0, 22, 45, 0, 0, 22, 52, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    data = UV(NAME, sorted_class, DISK, 1)
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.antenna = [None, 8]
    uvflg.timeran = [None, 0, 22, 25, 0, 0, 22, 35, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    
    data = UV(NAME, sorted_class, DISK, 1)
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.antennas = [None, 8]
    uvflg.timeran = [None, 0, 22, 0, 0, 1, 0, 0, 0]
    uvflg.reason = 'Badamp'
    uvflg.opcode = 'FLAG'
    uvflg.go()
    
    data = UV(NAME, sorted_class, DISK, 1)
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.antennas = [None, 1]
    uvflg.timeran = [None, 0, 17, 0, 0, 0, 18, 0, 0]
    uvflg.reason = 'Badamp'
    uvflg.opcode = 'FLAG'
    uvflg.go()
    
    data = UV(NAME, sorted_class, DISK, 1)
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, ampcal]
    uvflg.timeran = [None, 1, 0, 30, 0, 1, 2, 0, 0]
    uvflg.reason = 'Badrate'
    uvflg.opcode = 'FLAG'
    uvflg.go()
    
    data = UV(NAME, sorted_class, DISK, 1)
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.antenna = [None, 17]
    uvflg.timeran = [None, 0, 16, 50, 0, 0, 16, 55, 0]
    uvflg.reason = 'dropout'
    uvflg.opcode = 'FLAG'
    uvflg.go()
    
    data = UV(NAME, sorted_class, DISK, 1)
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.antenna = [None, 17]
    uvflg.baselin = [None, 7]
    uvflg.timeran = [None, 0, 16, 50, 0, 0, 16, 55, 0]
    uvflg.reason = 'BADBASEL'
    uvflg.opcode = 'FLAG'
    uvflg.go()
    
    data = UV(NAME, sorted_class, DISK, 1)
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.antenna = [None, 14]
    uvflg.timeran = [None, 0, 12, 30, 0, 0, 12, 41, 0]
    uvflg.reason = 'BADAMP'
    uvflg.opcode = 'FLAG'
    uvflg.go()
    
    # Quack
    quack = task('quack')
    quack.default()
    quack.indata = data
    quack.opcode = 'ENDB'
    quack.antenna = [None, 17]
    quack.aparm[2] = 15.0/60
    quack.go()
    

######### Y is missing from pipeline ampcal, no ampcal available. So set it to ones and correct with selfcal. ####
if whattodo['add_YY_to_table']:
    cloutver = 4 # 2 is the one we correct, 3 contains pipeline FRING.
    data = UV(TABNAME, TABCLASS, DISK, 1)
    # Remove all CL tables higher than version cloutver-1.
    for i in range(data.table_highver('CL'), cloutver-1, -1):
        data.zap_table('CL', i)
    indata = UV(TABNAME, TABCLASS, DISK, 1)
    outdata = UV(TABNAME, TABCLASS, DISK, 1)
    tacop = task('tacop')
    tacop.default()
    tacop.indata = indata
    tacop.outdata = outdata
    tacop.inext = 'CL'
    tacop.inver = 2
    tacop.outver = 4
    tacop.ncount = 1
    tacop.go()
    # Now add ones for antenna Y (16) by copy-paste
    # from CL1 which is naive CL from pipeline with no 
    # corrections.
    wizard_data = WUV(TABNAME, TABCLASS, DISK, 1)
    oldcl = wizard_data.table('CL', 1)
    newcl = wizard_data.table('CL', 4)
    
    #Copy entries of ones for AR to the corrected table.
    for row in oldcl:
        if row.antenna_no == 16: # Y is antenna 16 in this dataset.
            newcl.append(row)
    oldcl.close()
    newcl.close()
   
    # Add parrallactic corrections for this antenna
    data = UV(TABNAME, TABCLASS, DISK, 1)
    clcor = task('clcor')
    clcor.default()
    clcor.indata = data
    clcor.antenna = [None, 16]
    clcor.gainver = cloutver
    clcor.gainuse = cloutver
    clcor.opcode = 'PANG'
    clcor.clcorprm = [None, 1, 0] 
    if not os.path.exists('./usno_finals.erp'):
     os.system('wget http://gemini.gsfc.nasa.gov/solve_save/usno_finals.erp')
    clcor.infile = 'PWD:usno_finals.erp'
    clcor.go()

########## Use a priori Amp Cal an Parallactic angle correction from Jive ##########
if whattodo['copy_JIVE_CL']:
    indata = UV(TABNAME, TABCLASS, DISK, 1)
    outdata = UV(NAME, sorted_class, DISK, 1)
    cloutver = 2
    # Remove all CL tables higher than version cloutver-1.
    for i in range(outdata.table_highver('CL'), cloutver-1, -1):
        outdata.zap_table('CL', i)
    tacop = task('tacop')
    tacop.default()
    tacop.indata = indata
    tacop.outdata = outdata
    tacop.inext = 'CL'
    tacop.inver = 4 # With Y
    tacop.ncount = 1
    tacop.outver = cloutver
    tacop.go()

#   FRING on AMPCAL
if whattodo['fring_ampcal']:
    clinver = 2
    cloutver = 3
    snoutver = 1
    data = UV(NAME, sorted_class, DISK, 1)
    # Remove all SN tables higher than version snoutver-1.
    for i in range(data.table_highver('SN'), snoutver-1, -1):
        data.zap_table('SN', i)
    data = UV(NAME, sorted_class, DISK, 1)
    fring=task('fring')
    fring.default()
    fring.indata = data
    fring.docalib = 1
    fring.gainuse = clinver
    fring.calsour = [None, ampcal]
    fring.refant = refant
    fring.solint = 5 # Longer than scan
    fring.snver = snoutver
    #fring.dparm[2] = 200 # ns
    fring.go()

    # Apply SN-table
    data = UV(NAME, sorted_class, DISK, 1)
    # Remove all CL tables higher than version cloutver-1.
    for i in range(data.table_highver('CL'), cloutver-1, -1):
        data.zap_table('CL', i)
    data = UV(NAME, sorted_class, DISK, 1)
    # Make CL table
    clcal = task('clcal') 
    clcal.default()
    clcal.indata = data
    clcal.snver = snoutver
    clcal.invers = snoutver
    clcal.calsour = [None, ampcal]
    clcal.sour = [None, ''] # All sources
    clcal.gainver = clinver
    clcal.gainuse = cloutver
    clcal.refant = refant
    clcal.go()

if whattodo['bpass_ampcal']:
    clinver = 3
    bpoutver = 1
    data = UV(NAME, sorted_class, DISK, 1)
    # Remove all BP tables higher than version bpoutver-1.
    for i in range(data.table_highver('BP'), bpoutver-1, -1):
        data.zap_table('BP', i)
    data = UV(NAME, sorted_class, DISK, 1)
    bpass = task('bpass')
    bpass.default()
    bpass.indata = data
    bpass.calsour = [None,ampcal]
    bpass.timeran = [None, 0] # Use all scans to get all ants
    bpass.gainuse = clinver
    bpass.docal = 1
    bpass.refant = 1
    bpass.solint = -1 # Whole timerange
    bpass.flagver = 0
    bpass.outvers = bpoutver
    bpass.doband = -1
    bpass.bpassprm[10]=1 # Normalize amps using all channels
    bpass.go()

if whattodo['image_ampcal_1']:
    clinver = 3
    imagename = AMPCALIM + '1'
    imagedisk = 1
    imageseq = 1
    data = UV(NAME, sorted_class, DISK, 1)
    imagr = task('imagr')
    imagr.default()
    imagr.indata = data
    imagr.sources = [None, ampcal]
    imagr.gainuse = clinver
    imagr.flagver = 0
    imagr.nchav = 512
    imagr.stokes = 'I'
    imagr.cellsize = [None, 0.1e-3, 0.1e-3]
    imagr.imsize = [None, 512,512]
    imagr.docalib = 1
    imagr.dotv = -1
    imagr.outname = imagename
    imagr.outdisk = imagedisk
    imagr.outseq = imageseq
    imagr.niter = 500
    imagr.doband = 1
    imagr.bpver = 1
    #imagr.uvrange = [None, 7000,0]
    imagr.antenna = vlba_antennas
    imagr.baselin = imagr.antenna
    if imagr.niter == 0:
        imageclass = 'IIM001'
    else:
        imageclass = 'ICL001'
    # Remove previous old image
    image = IM(imagename, imageclass, imagedisk, imageseq)
    if image.exists():
     image.zap()
    # Also remove old beam image
    image = IM(imagename, 'IBM001', imagedisk, imageseq)
    if image.exists():
     image.zap()
    imagr.robust = 0.5
    imagr.go()


if whattodo['selfcal_ampcal']:
    clinver = 3
    cloutver = 4 # 
    snoutver = 2 # and 3
    data = UV(NAME, sorted_class, DISK, 1)
    # Remove all SN tables higher than version snoutver-1.
    for i in range(data.table_highver('SN'), snoutver-1, -1):
        data.zap_table('SN', i)
    data = UV(NAME, sorted_class, DISK, 1)
    calib = task('calib')
    calib.default()
    calib.indata = data
    calib.in2data = IM(AMPCALIM + '1', 'ICL001', 1, 1)
    calib.calsour = [None, ampcal]
    calib.refant = refant
    calib.docalib = 1
    calib.gainuse = clinver
    calib.flagver = 0
    calib.snver = snoutver
    calib.solmode = 'A&P'
    calib.solint = 5 # 
    calib.aparm[7] = 5 # SNR
    #calib.cmethod = 'DFT'
    calib.aparm[6] = 1 # Print level
    calib.weightit = 1
    calib.soltype = 'L1R'
    calib.doband = 1
    calib.bpver = 1
    calib.refant = 12
    calib.aparm[1] = 3 # Min no antennas
    calib.dofit = [None, 16, 17] #
    #calib.aparm[9] = 1 # Keep failed, has to since else flag many ants
    calib.go()
    # CLIP solutions
    data = UV(NAME, sorted_class, DISK, 1)
    snsmo = task('snsmo')
    snsmo.default()
    snsmo.indata = data
    snsmo.samptype = 'MWF'
    snsmo.cparm[1] = 2 # hours, clipping time ampl
    snsmo.cparm[2] = 2 # hours, clipping time phas
    snsmo.cparm[6] = 0.2 # Max deviation, ampl
    snsmo.cparm[7] = 400 # degrees Max deviation, phas
    snsmo.antenna = vlba_antennas # Do not smooth 16,17
    snsmo.inver = snoutver
    snsmo.outver = snoutver + 1
    snsmo.smotype = 'BOTH' # 
    snsmo.refant = 12
    snsmo.doblank = -1
    snsmo.go()
    # Apply SN-table
    data = UV(NAME, sorted_class, DISK, 1)
    # Remove all CL tables higher than version cloutver-1.
    for i in range(data.table_highver('CL'), cloutver-1, -1):
        data.zap_table('CL', i)
    data = UV(NAME, sorted_class, DISK, 1)
    # Make CL table
    clcal = task('clcal') 
    clcal.default()
    clcal.indata = data
    clcal.snver = snoutver+1
    clcal.invers = snoutver+1
    clcal.calsour = [None, ampcal] # Use sols from Ampcal 
    clcal.sour = [None, ''] # Apply to all sources
    clcal.gainver = clinver
    clcal.gainuse = cloutver
    clcal.refant = refant
    clcal.go()
    ## Fit EVN, no EB
    #clinver = 4
    #cloutver = 5
    #snoutver = 4
    #data = UV(NAME, sorted_class, DISK, 1)
    #calib = task('calib')
    #calib.default()
    #calib.indata = data
    #calib.in2data = IM(AMPCALIM + '1', 'ICL001', 1, 1)
    #calib.calsour = [None, ampcal]
    #calib.refant = refant
    #calib.docalib = 1
    #calib.gainuse = clinver
    #calib.flagver = 0
    #calib.snver = snoutver
    #calib.solmode = 'A&P'
    #calib.solint = 5 # 
    #calib.aparm[7] = 5 # SNR
    ##calib.cmethod = 'DFT'
    #calib.aparm[6] = 1 # Print level
    #calib.weightit = 1
    #calib.soltype = 'L1R'
    #calib.doband = 1
    #calib.bpver = 1
    #calib.refant = 12
    #calib.aparm[1] = 4 # Min no antennas
    #calib.dofit = evn_antennas
    #calib.aparm[3] = 1 # need to average pols, since few ants here with LL
    #calib.aparm[9] = 1 # Keep failed, has to since else flag many ants
    #calib.go()
    ## Apply SN-table
    #data = UV(NAME, sorted_class, DISK, 1)
    ## Remove all CL tables higher than version cloutver-1.
    #for i in range(data.table_highver('CL'), cloutver-1, -1):
    #    data.zap_table('CL', i)
    #data = UV(NAME, sorted_class, DISK, 1)
    ## Make CL table
    #clcal = task('clcal') 
    #clcal.default()
    #clcal.indata = data
    #clcal.snver = snoutver
    #clcal.invers = snoutver
    #clcal.calsour = [None, ampcal] # Use sols from Ampcal 
    #clcal.sour = [None, ''] # Apply to all sources
    #clcal.gainver = clinver
    #clcal.gainuse = cloutver
    #clcal.refant = refant
    #clcal.go()
    ## Fit EB
    #clinver = 5
    #cloutver = 6
    #snoutver = 4
    #data = UV(NAME, sorted_class, DISK, 1)
    #calib = task('calib')
    #calib.default()
    #calib.indata = data
    #calib.in2data = IM(AMPCALIM + '1', 'ICL001', 1, 1)
    #calib.calsour = [None, ampcal]
    #calib.refant = refant
    #calib.docalib = 1
    #calib.gainuse = clinver
    #calib.flagver = 0
    #calib.snver = snoutver
    #calib.solmode = 'A&P'
    #calib.solint = 5 # 
    #calib.aparm[7] = 5 # SNR
    ##calib.cmethod = 'DFT'
    #calib.aparm[6] = 1 # Print level
    #calib.weightit = 1
    #calib.soltype = 'L1R'
    #calib.doband = 1
    #calib.bpver = 1
    #calib.refant = 12
    #calib.aparm[1] = 4 # Min no antennas
    #calib.dofit = [None, 1]
    #calib.aparm[3] = 1 # need to average pols, since few ants here with LL
    #calib.aparm[9] = 1 # Keep failed, has to since else flag many ants
    #calib.go()
    ## Apply SN-table
    #data = UV(NAME, sorted_class, DISK, 1)
    ## Remove all CL tables higher than version cloutver-1.
    #for i in range(data.table_highver('CL'), cloutver-1, -1):
    #    data.zap_table('CL', i)
    #data = UV(NAME, sorted_class, DISK, 1)
    ## Make CL table
    #clcal = task('clcal') 
    #clcal.default()
    #clcal.indata = data
    #clcal.snver = snoutver
    #clcal.invers = snoutver
    #clcal.calsour = [None, ampcal] # Use sols from Ampcal 
    #clcal.sour = [None, ''] # Apply to all sources
    #clcal.gainver = clinver
    #clcal.gainuse = cloutver
    #clcal.refant = refant
    #clcal.go()

if whattodo['clip_ampcal']:
    outfgver = 2
    clinver = 4
    data = UV(NAME, sorted_class, DISK, 1)
    clip = task('clip')
    clip.default()
    clip.sour = [None, ampcal]
    clip.aparm[6] = 1.5e-3 # Flag very sensitive baselines
    clip.docal = 1
    clip.gainuse = clinver
    clip.stokes = 'I'
    clip.flagver = outfgver
    clip.doband = 1
    clip.bpver = 1
    clip.indata = data
    clip.outfgver = outfgver
    clip.go()

if whattodo['image_ampcal_2']:
    clinver = 4
    imagename = AMPCALIM+'2'
    imagedisk = 1
    imageseq = 1
    data = UV(NAME, sorted_class, DISK, 1)
    imagr = task('imagr')
    imagr.default()
    imagr.indata = data
    imagr.sources = [None, ampcal]
    imagr.gainuse = clinver
    imagr.flagver = 0
    imagr.nchav = 512
    imagr.stokes = 'I'
    imagr.cellsize = [None, 0.1e-3, 0.1e-3]
    imagr.imsize = [None, 512,512]
    imagr.docalib = 1
    imagr.dotv = -1
    imagr.doband = 1
    imagr.bpver = 1
    imagr.outname = imagename
    imagr.outdisk = imagedisk
    imagr.outseq = imageseq
    imagr.niter = 300
    imagr.uvrange = [None, 7000, 0]
    #imagr.antenna = vlba_antennas
    imagr.baselin = imagr.antenna
    if imagr.niter == 0:
        imageclass = 'IIM001'
    else:
        imageclass = 'ICL001'
    # Remove previous old image
    image = IM(imagename, imageclass, imagedisk, imageseq)
    if image.exists():
     image.zap()
    # Also remove old beam image
    image = IM(imagename, 'IBM001', imagedisk, imageseq)
    if image.exists():
     image.zap()
    imagr.robust = 0.5
    imagr.go()

#   FRING on PHASECAL
if whattodo['fring_phasecal']:
    clinver = 4 # With ampcal corrections
    cloutver = 5
    snoutver = 4
    refant = 1
    data = UV(NAME, sorted_class, DISK, 1)
    # Remove all SN tables higher than version snoutver-1.
    for i in range(data.table_highver('SN'), snoutver-1, -1):
        data.zap_table('SN', i)
    data = UV(NAME, sorted_class, DISK, 1)
    fring=task('fring')
    fring.default()
    fring.indata = data
    fring.docalib = 1
    fring.gainuse = clinver
    fring.calsour = [None, phasecal]
    fring.refant = refant
    fring.solint = 3 # Longer than scan
    fring.dparm[2] = 200#ns, delay win, already fitted bulk on ampcal
    fring.dparm[3] = 50#mHz, rate win
    #fring.aparm[3] = 1 # AVG RR,LL
    #fring.aparm[5] = 1 # AVG IFs
    fring.snver = snoutver
    fring.doband = 1
    fring.bpver = 1
    fring.go()

    # Apply SN-table
    data = UV(NAME, sorted_class, DISK, 1)
    # Remove all CL tables higher than version cloutver-1.
    for i in range(data.table_highver('CL'), cloutver-1, -1):
        print "WARNING: Deleting CL" + str(i)
        data.zap_table('CL', i)
    data = UV(NAME, sorted_class, DISK, 1)
    # Make CL table
    clcal = task('clcal') 
    clcal.default()
    clcal.indata = data
    clcal.snver = snoutver
    clcal.invers = snoutver
    clcal.calsour = [None, phasecal]
    clcal.sour = [None, phasecal, target]
    clcal.gainver = clinver
    clcal.gainuse = cloutver
    clcal.refant = refant
    clcal.go()

if whattodo['clip_phasecal']:
    outfgver = 2
    data = UV(NAME, sorted_class, DISK, 1)
    clip = task('clip')
    clip.default()
    clip.sour = [None, phasecal, target]
    clip.aparm[1] = 1.5 # Max Jy
    clip.stokes = 'I'
    clip.docal = 1
    clip.flagver = outfgver
    clip.doband = 1
    clip.bpver = 1
    clip.indata = data
    clip.outfgver = outfgver
    clip.go()

if whattodo['image_phasecal_1']:
    clinver = 5
    imagename = PCALIM + '1'
    imagedisk = 1
    imageseq = 1
    data = UV(NAME, sorted_class, DISK, 1)
    imagr = task('imagr')
    imagr.default()
    imagr.indata = data
    imagr.sources = [None, phasecal]
    imagr.gainuse = clinver
    imagr.flagver = 0
    imagr.nchav = 512
    imagr.stokes = 'I'
    imagr.cellsize = [None, 0.1e-3, 0.1e-3]
    imagr.imsize = [None, 512,512]
    imagr.docalib = 1
    imagr.dotv = -1
    imagr.doband = 1
    imagr.bpver = 1
    imagr.outname = imagename
    imagr.outdisk = imagedisk
    imagr.outseq = imageseq
    imagr.niter = 500
    imagr.uvrang = [None, 7000, 0]
    if imagr.niter == 0:
        imageclass = 'IIM001'
    else:
        imageclass = 'ICL001'
    # Remove previous old image
    image = IM(imagename, imageclass, imagedisk, imageseq)
    if image.exists():
     image.zap()
    # Also remove old beam image
    image = IM(imagename, 'IBM001', imagedisk, imageseq)
    if image.exists():
     image.zap()
    imagr.robust = 0.5
    imagr.go()


if whattodo['selfcal_phasecal']:
    clinver = 5
    cloutver = 6
    snoutver = 5 # and 6
    data = UV(NAME, sorted_class, DISK, 1)
    # Remove all SN tables higher than version snoutver-1.
    for i in range(data.table_highver('SN'), snoutver-1, -1):
        print "WARNING: Deleting SN" + str(i)
        data.zap_table('SN', i)
    data = UV(NAME, sorted_class, DISK, 1)
    calib = task('calib')
    calib.default()
    calib.indata = data
    calib.in2data = IM(PCALIM + '1', 'ICL001', 1, 1)
    calib.calsour = [None, phasecal]
    calib.refant = refant
    calib.docalib = 1
    calib.gainuse = clinver
    calib.flagver = 0
    calib.snver = snoutver
    calib.solmode = 'A&P'
    calib.solint = 3 # Longer than scan
    calib.aparm[1] = 4 # Minant
    calib.aparm[7] = 5 # SNR
    #calib.aparm[3] = 1 # AVG RR,LL, 
    calib.aparm[6] = 1 # Print level
    #calib.aparm[5] = 1 # AVG ifs
    #calib.dofit = vlba_antennas # Cannot calibrate EVN easily,will lose LL
    calib.weightit = 1
    calib.soltype = 'L1R'
    calib.doband = 1
    calib.bpver = 1
    calib.uvran = [None, 7000,0]
    calib.go()

    # CLIP solutions
    data = UV(NAME, sorted_class, DISK, 1)
    snsmo = task('snsmo')
    snsmo.default()
    snsmo.indata = data
    snsmo.samptype = 'MWF'
    snsmo.cparm[1] = 2 # hours, clipping time ampl
    snsmo.cparm[6] = 0.2 # Max deviation, ampl
    snsmo.inver = snoutver
    snsmo.outver = snoutver + 1
    snsmo.smotype = 'AMPL' # 
    snsmo.refant = refant
    snsmo.doblank = -1
    snsmo.go()
    
    # Apply SN-table
    data = UV(NAME, sorted_class, DISK, 1)
    # Remove all CL tables higher than version cloutver-1.
    for i in range(data.table_highver('CL'), cloutver-1, -1):
        print "WARNING: Deleting CL" + str(i)
        data.zap_table('CL', i)
    data = UV(NAME, sorted_class, DISK, 1)
    # Make CL table
    clcal = task('clcal') 
    clcal.default()
    clcal.indata = data
    clcal.snver = snoutver + 1
    clcal.invers = snoutver + 1
    clcal.calsour = [None, phasecal]
    clcal.sour = [None, phasecal, target] # Apply to all sources
    clcal.gainver = clinver
    clcal.gainuse = cloutver
    clcal.refant = refant
    clcal.go()

if whattodo['image_phasecal_2']:
    clinver = 6
    imagename = PCALIM + '2'
    imagedisk = 1
    imageseq = 1
    data = UV(NAME, sorted_class, DISK, 1)
    imagr = task('imagr')
    imagr.default()
    imagr.indata = data
    imagr.sources = [None, phasecal]
    imagr.gainuse = clinver
    imagr.flagver = 0
    imagr.nchav = 512
    imagr.stokes = 'I'
    imagr.cellsize = [None, 0.1e-3, 0.1e-3]
    imagr.imsize = [None, 512,512]
    imagr.docalib = 1
    imagr.dotv = -1
    imagr.doband = 1
    imagr.bpver = 1
    imagr.outname = imagename
    imagr.outdisk = imagedisk
    imagr.outseq = imageseq
    imagr.niter = 300
    imagr.uvran = [None, 7000,0]
    if imagr.niter == 0:
        imageclass = 'IIM001'
    else:
        imageclass = 'ICL001'
    # Remove previous old image
    image = IM(imagename, imageclass, imagedisk, imageseq)
    if image.exists():
     image.zap()
    # Also remove old beam image
    image = IM(imagename, 'IBM001', imagedisk, imageseq)
    if image.exists():
     image.zap()
    imagr.robust = 0.5
    imagr.go()

#   FRING on FRCAL
if whattodo['fring_frcal']:
    clinver = 6
    cloutver = 7
    snoutver = 7
    data = UV(NAME, sorted_class, DISK, 1)
    # Remove all SN tables higher than version snoutver-1.
    for i in range(data.table_highver('SN'), snoutver-1, -1):
        data.zap_table('SN', i)
    data = UV(NAME, sorted_class, DISK, 1)
    fring=task('fring')
    fring.default()
    fring.indata = data
    fring.docalib = 1
    fring.gainuse = clinver
    fring.calsour = [None, frcal]
    fring.refant = refant
    fring.solint = 5 # Longer than scan
    fring.snver = snoutver
    fring.go()

    # Apply SN-table
    data = UV(NAME, sorted_class, DISK, 1)
    # Remove all CL tables higher than version cloutver-1.
    for i in range(data.table_highver('CL'), cloutver-1, -1):
        data.zap_table('CL', i)
    data = UV(NAME, sorted_class, DISK, 1)
    # Make CL table
    clcal = task('clcal') 
    clcal.default()
    clcal.indata = data
    clcal.snver = snoutver
    clcal.invers = snoutver
    clcal.calsour = [None, frcal]
    clcal.sour = [None, frcal]
    clcal.gainver = clinver
    clcal.gainuse = cloutver
    clcal.refant = refant
    clcal.go()

if whattodo['image_frcal_1']:
    clinver = 7
    imagename = FRCALIM + '1'
    imagedisk = 1
    imageseq = 1
    data = UV(NAME, sorted_class, DISK, 1)
    imagr = task('imagr')
    imagr.default()
    imagr.indata = data
    imagr.sources = [None, frcal]
    imagr.gainuse = clinver
    imagr.flagver = 0
    imagr.nchav = 512
    imagr.stokes = 'I'
    imagr.cellsize = [None, 0.1e-3, 0.1e-3]
    imagr.imsize = [None, 512,512]
    imagr.docalib = 1
    imagr.dotv = -1
    imagr.outname = imagename
    imagr.outdisk = imagedisk
    imagr.outseq = imageseq
    imagr.niter = 300
    imagr.doband = 1
    imagr.bpver = 1
    imagr.uvrange = [None, 7000,0]
    imagr.robust = 0.5
    imagr.clbox = [None, [None, 204,189,302, 313]]
    imagr.nbox = 1
    if imagr.niter == 0:
        imageclass = 'IIM001'
    else:
        imageclass = 'ICL001'
    # Remove previous old image
    image = IM(imagename, imageclass, imagedisk, imageseq)
    if image.exists():
     image.zap()
    # Also remove old beam image
    image = IM(imagename, 'IBM001', imagedisk, imageseq)
    if image.exists():
     image.zap()
    imagr.go()

if whattodo['selfcal_frcal']:
    clinver = 7
    cloutver = 8
    snoutver = 8
    data = UV(NAME, sorted_class, DISK, 1)
    # Remove all SN tables higher than version snoutver-1.
    for i in range(data.table_highver('SN'), snoutver-1, -1):
        data.zap_table('SN', i)
    # FIT ANTENNA 3, with proper timeran not to flag things
    data = UV(NAME, sorted_class, DISK, 1)
    calib = task('calib')
    calib.default()
    calib.indata = data
    calib.in2data = IM(FRCALIM + '1', 'ICL001', 1, 1)
    calib.calsour = [None, frcal]
    calib.refant = refant
    calib.docalib = 1
    calib.gainuse = clinver
    calib.flagver = 0
    calib.snver = snoutver
    calib.solmode = 'A&P'
    calib.solint = 0.5 # 
    calib.aparm[7] = 5 # SNR
    calib.cmethod = 'DFT'
    calib.aparm[6] = 1 # Print level
    calib.aparm[9] = 1 # Pass on failed sols, for scans without these ants
    calib.doband = 1
    calib.bpver = 1
    calib.go()
    # Apply SN-table
    data = UV(NAME, sorted_class, DISK, 1)
    # Remove all CL tables higher than version cloutver-1.
    for i in range(data.table_highver('CL'), cloutver-1, -1):
        data.zap_table('CL', i)
    data = UV(NAME, sorted_class, DISK, 1)
    # Make CL table
    clcal = task('clcal') 
    clcal.default()
    clcal.indata = data
    clcal.snver = snoutver
    clcal.invers = snoutver
    clcal.calsour = [None, frcal]
    clcal.sour = [None, frcal]
    clcal.gainver = clinver
    clcal.gainuse = cloutver
    clcal.refant = refant
    clcal.go()

if whattodo['image_frcal_2']:
    clinver = 8
    imagename = FRCALIM+'2'
    imagedisk = 1
    imageseq = 1
    data = UV(NAME, sorted_class, DISK, 1)
    imagr = task('imagr')
    imagr.default()
    imagr.indata = data
    imagr.sources = [None, frcal]
    imagr.gainuse = clinver
    imagr.flagver = 0
    imagr.nchav = 512
    imagr.stokes = 'I'
    imagr.cellsize = [None, 0.1e-3, 0.1e-3]
    imagr.imsize = [None, 512, 512]
    imagr.docalib = 1
    imagr.dotv = -1
    imagr.doband = 1
    imagr.bpver = 1
    imagr.outname = imagename
    imagr.outdisk = imagedisk
    imagr.outseq = imageseq
    imagr.niter = 300
    imagr.uvrange = [None, 7000,0]
    imagr.robust = 0.5
    imagr.clbox = [None, [None, 204,189,302, 313]]
    imagr.nbox = 1
    if imagr.niter == 0:
        imageclass = 'IIM001'
    else:
        imageclass = 'ICL001'
    # Remove previous old image
    image = IM(imagename, imageclass, imagedisk, imageseq)
    if image.exists():
     image.zap()
    # Also remove old beam image
    image = IM(imagename, 'IBM001', imagedisk, imageseq)
    if image.exists():
     image.zap()
    imagr.go()


############### IMAGE ARP220 ######################
# Load these outside if to be available also for export later
arpfilename = TCALIM
arpfiledisk = 1
arpfileseq = 1
arpfileclass = 'ICL001'
arpfileclass2 = 'ICL002'
if whattodo['image_ARP220']:
    clinver = 8
    # FIRST CLEAN WITH SMALL BOXES
    # Remove previous old images
    image = IM(arpfilename, arpfileclass, arpfiledisk, arpfileseq)
    if image.exists():
        image.zap()
    image = IM(arpfilename, arpfileclass2, arpfiledisk, arpfileseq)
    if image.exists():
        image.zap()
    # Also remove old beam images
    image = IM(arpfilename, 'IBM001', arpfiledisk, arpfileseq)
    if image.exists():
        image.zap()
    image = IM(arpfilename, 'IBM002', arpfiledisk, arpfileseq)
    if image.exists():
        image.zap()
    # Clean with latest gainuse
    data = UV(NAME, sorted_class, DISK, 1)
    imagr = task('imagr')
    imagr.default()
    imagr.indata = data
    imagr.sources = [None, target]
    imagr.im2parm = [None, 1, 5.0, 10.0, 0] # Require island to have SNR 5 and peak SNR10. To clean (for later removal) of the strongest sources
    imagr.gainuse = clinver
    imagr.flagver = 0
    imagr.nchav = 512
    imagr.doband = 1
    imagr.bpver = 1
    imagr.imsize = [None, 8192,4096]
    imagr.cellsize = [None, 0.1e-3, 0.1e-3]
    imagr.stokes = 'I'
    imagr.nfield = 2
    imagr.overlap = 1
    imagr.do3dimag = 1
    imagr.rashift = [None, 0.6, -0.35]
    imagr.decshift = [None, 0.02, 0.18]
    imagr.docalib = 1
    imagr.dotv = -1
    imagr.outname = arpfilename
    imagr.outdisk = arpfiledisk
    imagr.outseq = arpfileseq
    imagr.niter = 500
    imagr.gain = 0.1
    imagr.uvrange = [None, 7000,0] # remove large-scale disturbances
    imagr.robust = 0.5
    imagr.go()
    # CLEAN WITH LARGE (NO) BOXES
    arpboxfile = ''
    data = UV(NAME, sorted_class, DISK, 1)
    imagr = task('imagr')
    imagr.default()
    imagr.indata = data
    imagr.sources = [None, target]
    imagr.gainuse = clinver
    imagr.flagver = 0
    imagr.nchav = 512
    imagr.doband = 1
    imagr.bpver = 1
    imagr.imsize = [None, 8192,4096]
    imagr.cellsize = [None, 0.1e-3, 0.1e-3]
    imagr.stokes = 'I'
    imagr.nfield = 2
    imagr.overlap = 1
    imagr.do3dimag = 1
    imagr.rashift = [None, 0.6, -0.35]
    imagr.decshift = [None, 0.02, 0.18]
    imagr.docalib = 1
    imagr.dotv = -1
    imagr.outname = arpfilename
    imagr.outdisk = arpfiledisk
    imagr.outseq = arpfileseq
    imagr.outver = 1
    imagr.niter = 1000
    imagr.gain = 0.1
    imagr.bcomp = [None, 100000, 100000]
    imagr.uvrange = [None, 7000,0] # remove large-scale disturbances
    imagr.robust = 0.5
    imagr.go()


if whattodo['selfcal_arp220']:
    clinver = 8
    cloutver = 9
    snoutver = 9
    data = UV(NAME, sorted_class, DISK, 1)
    # Remove all SN tables higher than version snoutver-1.
    for i in range(data.table_highver('SN'), snoutver-1, -1):
        print "WARNING: Deleting SN" + str(i)
        data.zap_table('SN', i)
    data = UV(NAME, sorted_class, DISK, 1)
    calib = task('calib')
    calib.default()
    calib.indata = data
    calib.in2data = IM(TCALIM, 'ICL001', 1, 1)
    calib.calsour = [None, target]
    calib.nmaps = 2
    calib.refant = refant
    calib.docalib = 1
    calib.gainuse = clinver
    calib.flagver = 0
    calib.snver = snoutver
    calib.solmode = 'P'
    calib.solint = 10 # Longer than scan
    calib.aparm[1] = 3 # min No. of antennas 
    calib.aparm[7] = 2 # SNR
    calib.aparm[3] = 1 # AVG RR,LL, 
    calib.aparm[6] = 1 # Print level
    calib.aparm[5] = 1 # AVG ifs
    calib.cmethod = 'DFT'
    calib.weightit = 1
    calib.soltype = 'L1R'
    calib.doband = 1
    calib.bpver = 1
    calib.uvrange = [None, 7000,0]
    calib.go()
    
    # Apply SN-table
    data = UV(NAME, sorted_class, DISK, 1)
    # Remove all CL tables higher than version cloutver-1.
    for i in range(data.table_highver('CL'), cloutver-1, -1):
        print "WARNING: Deleting CL" + str(i)
        data.zap_table('CL', i)
    data = UV(NAME, sorted_class, DISK, 1)
    # Make CL table
    clcal = task('clcal') 
    clcal.default()
    clcal.indata = data
    clcal.snver = snoutver
    clcal.invers = snoutver
    clcal.calsour = [None, target]
    clcal.sour = [None, target] # Apply to target
    clcal.gainver = clinver
    clcal.gainuse = cloutver
    clcal.refant = refant
    clcal.go()

############### IMAGE ARP220 ######################
# Load these outside if to be available also for export later
arpfilename = OUTPREFIX + '_FIN'
arpfiledisk = 1
arpfileseq = 1
arpfileclass = 'ICL001'
arpfileclass2 = 'ICL002'
if whattodo['image_ARP220_final']:
    clinver = 9
    # FIRST CLEAN WITH SMALL BOXES
    # Remove previous old images
    image = IM(arpfilename, arpfileclass, arpfiledisk, arpfileseq)
    if image.exists():
        image.zap()
    image = IM(arpfilename, arpfileclass2, arpfiledisk, arpfileseq)
    if image.exists():
        image.zap()
    # Also remove old beam images
    image = IM(arpfilename, 'IBM001', arpfiledisk, arpfileseq)
    if image.exists():
        image.zap()
    image = IM(arpfilename, 'IBM002', arpfiledisk, arpfileseq)
    if image.exists():
        image.zap()
    # Clean with latest gainuse
    data = UV(NAME, sorted_class, DISK, 1)
    imagr = task('imagr')
    imagr.default()
    imagr.indata = data
    imagr.sources = [None, target]
    imagr.im2parm = [None, 1, 5.0, 10.0, 0] # Require island to have SNR 5 and peak SNR10. To clean (for later removal) of the strongest sources
    imagr.gainuse = clinver
    imagr.flagver = 0
    imagr.nchav = 512
    imagr.doband = 1
    imagr.bpver = 1
    imagr.imsize = [None, 8192, 8192]
    imagr.cellsize = [None, 0.1e-3, 0.1e-3]
    imagr.stokes = 'I'
    imagr.nfield = 2
    imagr.overlap = 1
    imagr.do3dimag = 1
    imagr.rashift = [None, 0.6, -0.35]
    imagr.decshift = [None, 0.02, 0.18]
    imagr.docalib = 1
    imagr.dotv = -1
    imagr.outname = arpfilename
    imagr.outdisk = arpfiledisk
    imagr.outseq = arpfileseq
    imagr.niter = 500
    imagr.gain = 0.1
    imagr.uvrange = [None, 7000,0] # remove large-scale disturbances
    imagr.robust = 0.5
    imagr.go()
    # CLEAN WITH LARGE (NO) BOXES
    arpboxfile = ''
    data = UV(NAME, sorted_class, DISK, 1)
    imagr = task('imagr')
    imagr.default()
    imagr.indata = data
    imagr.sources = [None, target]
    imagr.gainuse = clinver
    imagr.flagver = 0
    imagr.nchav = 512
    imagr.doband = 1
    imagr.bpver = 1
    imagr.imsize = [None, 8192, 8192]
    imagr.cellsize = [None, 0.1e-3, 0.1e-3]
    imagr.stokes = 'I'
    imagr.nfield = 2
    imagr.overlap = 1
    imagr.do3dimag = 1
    imagr.rashift = [None, 0.6, -0.35]
    imagr.decshift = [None, 0.02, 0.18]
    imagr.docalib = 1
    imagr.dotv = -1
    imagr.outname = arpfilename
    imagr.outdisk = arpfiledisk
    imagr.outseq = arpfileseq
    imagr.outver = 1
    imagr.niter = 1500
    imagr.gain = 0.1
    imagr.bcomp = [None, 100000, 100000]
    imagr.uvrange = [None, 7000,0] # remove large-scale disturbances
    imagr.robust = 0.5
    imagr.go()

############### EXPORT RESULTS ######################
if whattodo['export_results']:
    east = IM(arpfilename, arpfileclass, arpfiledisk, arpfileseq)
    west = IM(arpfilename, arpfileclass2, arpfiledisk, arpfileseq)
   
    # Export EAST
    fittp = task('fittp')
    fittp.default()
    fittp.indata = east
    fittp.dataout = 'PWD:'+OUTPREFIX+'_EAST_IMAGE.FITS'
    if os.path.exists(fittp.dataout[4:]):
     os.system('rm ' + fittp.dataout[4:])
    fittp.go()
    # Export WEST
    fittp = task('fittp')
    fittp.default()
    fittp.indata = west
    fittp.dataout = 'PWD:'+OUTPREFIX+'_WEST_IMAGE.FITS'
    if os.path.exists(fittp.dataout[4:]):
     os.system('rm ' + fittp.dataout[4:])
    fittp.go()
    # Export FRCAL
    im = IM(NAME+'_FC2', 'ICL001', 1, 1)
    fittp = task('fittp')
    fittp.default()
    fittp.indata = im
    fittp.dataout = 'PWD:'+NAME+'_J1613+3412_IMAGE.FITS'
    if os.path.exists(fittp.dataout[4:]):
     os.system('rm ' + fittp.dataout[4:])
    fittp.go()
    # Export ACAL
    im = IM(NAME+'_AC2', 'ICL001', 1, 1)
    fittp = task('fittp')
    fittp.default()
    fittp.indata = im
    fittp.dataout = 'PWD:'+NAME+'_J1516+1932_IMAGE.FITS'
    if os.path.exists(fittp.dataout[4:]):
     os.system('rm ' + fittp.dataout[4:])
    fittp.go()
    # Export PCAL
    im = IM(NAME+'_PC2', 'ICL001', 1, 1)
    fittp = task('fittp')
    fittp.default()
    fittp.indata = im
    fittp.dataout = 'PWD:'+NAME+'_J1532+2344_IMAGE.FITS'
    if os.path.exists(fittp.dataout[4:]):
     os.system('rm ' + fittp.dataout[4:])
    fittp.go()
    
    fittp = task('fittp')
    fittp.default()
    fittp.indata = IM(NAME+'_FIN','IBM001', 1, 1)
    fittp.dataout = 'PWD:'+NAME+'_EAST_BEAM.FITS'
    if os.path.exists(fittp.dataout[4:]):
     os.system('rm ' + fittp.dataout[4:])
    fittp.go()
    fittp = task('fittp')
    fittp.default()
    fittp.indata = IM(NAME+'_FIN','IBM002', 1, 1)
    fittp.dataout = 'PWD:'+NAME+'_WEST_BEAM.FITS'
    if os.path.exists(fittp.dataout[4:]):
     os.system('rm ' + fittp.dataout[4:])
    fittp.go()

if whattodo['save_calsols']:
    snver = 9
    # Remove all PL tables
    data = UV(NAME, sorted_class, DISK, 1)

    data.zap_table('PL', -1)
    snplt = task('snplt')
    snplt.indata = data
    snplt.bif = 1
    snplt.eif = 1
    snplt.stokes = 'RR'
    snplt.optype = 'PHAS'
    snplt.sour = [None, target]
    snplt.inext = 'SN'
    snplt.inver = snver
    snplt.nplots = 6
    #snplt.pixrange = [None, -180,180]
    snplt.dotv = -1
    snplt.go()
    outname = NAME+'_ARP220_SELFCAL_SOLS.EPS'
    if os.path.exists(outname):
        os.system('rm ' + outname)
    lwpla = task('lwpla')
    lwpla.indata = data
    lwpla.plver = 1
    lwpla.inver = 99
    lwpla.outfile = 'PWD:'+outname
    #lwpla.docolor = 1
    lwpla.go()

tac = time.time()
processtime = tac-tic
print 'OVERALL PROCESS TIME: ',processtime,' SECONDS. or ', processtime/60.0, 'minutes'