# All files for project GL021 Lband were downloaded from the NRAO archive.
# Then renamed into a list of sequential filenames which AIPS FITLD could read
# easily. FITLD produces 6 files, presimably for different correlation passes.
# We use the first file, which has the continuum data for all stations,
# to make images of Arp220. 
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
from astropy.io import fits

#Ant   1 = EB       BX= -2716029.1182 BY=  4059907.1583 BZ=   973413.8253
#Ant   2 = HN       BX=  1741523.8390 BY=   716362.9947 BZ=   395289.2273
#Ant   5 = PT       BX=  1807322.6128 BY= -2421887.8920 BZ=  -351605.0083
#Ant   6 = Y        BX=  1840494.8333 BY= -2386975.7475 BZ=  -372140.9673
#Ant   7 = JV       BX= -2116984.1295 BY=  3748770.9591 BZ=  1159469.4243
#Ant   8 = NT       BX= -3395458.7510 BY=  5082351.0513 BZ=  -120532.4870
#Ant   9 = TR       BX= -3504782.1802 BY=  3787157.8325 BZ=  1150024.1625
#Ant  10 = WB       BX= -2705884.5374 BY=  3850285.8049 BZ=  1137906.2425
#Ant  11 = GO       BX=  1324673.4506 BY= -3065648.7048 BZ=  -249964.6369
#Ant  12 = KP       BX=  1772800.7696 BY= -2775645.3476 BZ=  -569688.7550
#Ant  13 = OV       BX=  1155119.2511 BY= -3094429.9497 BZ=   -88399.4985
#Ant  14 = BR       BX=   439371.1847 BY= -2677502.6691 BZ=   799796.9212
#Ant  15 = FD       BX=  2171294.5475 BY= -2159782.0902 BZ=  -695054.4020
#Ant  16 = MK       BX= -1291271.2903 BY= -5792820.8358 BZ= -1778720.1192
#Ant  17 = SC       BX=  2954027.8069 BY=  1696528.8254 BZ= -1994277.3759
#Ant  18 = MC       BX= -3074705.5538 BY=  4551017.6010 BZ=   522542.2383


vlba_antennas = [None, 2,3,4,5,12,13,14,15,16,17]
non_vlba_antennas = [None, 1,6,7,8,9,10,11,18]

masercal_bchan = 13
masercal_echan = 13
masercal_bif = 1
masercal_eif = 1
masercal_nchav = 512

whattodo = {'load_data': False, # ->CL1. Creates 6 files. Use the FIRST one for Lband continuum, skip others.
            'msort+indxr': False,
            'eops_corr': False, # ->CL2
            'preflagging': False,
            'par_ang_cor': False, # ->CL3
            'apriori_ampcal' : False, # SN1,2,3,CL4, # Doing the VLBA+VLA+EB
            'fring_cals' : True, #->SN4, CL5
            'bpass': True, # BP 1

            'image_ff_1': True, 
            'selfcal_ff' : True, # -> SN5, CL6 # Correct VLBA
            'image_ff_2': True, 
            'selfcal_ff_2' : True, # -> SN6, CL7 # Correct VLA+EB
            'image_ff_3': True, 

            'image_ampcal_1': True, 
            'selfcal_ampcal' : True, # -> SN7, CL8
            'image_ampcal_2': True, 
            
            'selfcal_target_maser' : True, # -> SN8, CL9. USING POINT MODEL. 
            'image_ARP220_maser_1': True,
            'selfcal_target_maser_2' : True, # -> SN9, CL10 Using image model, PHASE ONLY, 1min. Back from CL8.
            'image_ARP220_maser_2': True,
            'selfcal_target_maser_3' : True, # -> SN10 CL11. AMP + PHASE. Using model 2, long solint. Use CL10
            'image_ARP220_maser_3': True,
            'split': True,
            'fixwt': True,
            'multi': True,

            'flag_maser': True,
            'calculate_shifts': True, # Calculate shifts for Arp220 imaging. MUST be enabled to set variables used by imaging.
            'image_ARP220': True,
            'selfcal_arp220': True, # ->SN1, CL2
            'image_ARP220_V': True,
            ######
            'image_ARP220_final': True,
            'export_results': True,

            }

########## Initialize observation data ##########
AIPS.userno = 1000
clint = 1.0/60 # Seconds
NAME = 'GL021_L'
NAMERAW = NAME
CLASS = 'UVDATA'
MCLASS = 'L21L'
DISK = 1
AMPCALIM = NAME + '_AC'
FFCALIM = NAME + '_FC'
PCALIM = NAME + '_PC'
TCALIM = NAME + '_TC'
TCALIM2 = TCALIM + '2'
OUTPREFIX = NAME
target = 'ARP220'
ampcal = '1308+326'
ff = '1613+341'
uvrange = [None, 5000, 0]
sorted_class = 'MSORT'
logfile = 'log_'+OUTPREFIX+'_PTred.log'
refant = 2
#########^^^ 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(NAMERAW, CLASS, DISK, 1)
    fitld.datain = '/data1/eskil/RAWDATA/TEST_ARP220_ARCHIVE/GL021/file.uvfits1'
    #fitld.outdata = outdata
    fitld.outname = outdata.name
    #fitld.sour = [None, 'ARP220', '1613+341', '1308+326']
    fitld.ncount = 99
    fitld.digicor = 1
    fitld.clint = clint
    fitld.wtthresh = 0.7
    fitld.doconcat = 1
    fitld.douvcomp = -1
    if outdata.exists():
        print 'WARNING: Removing existing file ' + NAMERAW + '.' + CLASS
        outdata.zap()
        print '         File removed. Proceeding with FITLD...'
    fitld.go()

if whattodo['msort+indxr']:
    #for seq in [1,2,3,4,5,6 ]:
    for seq in [1,]: # We only use the first file
        msort = task('msort')
        data = UV(NAMERAW, CLASS, 1, seq)
        outdata = UV(NAMERAW, sorted_class, 1, seq)
        msort.default()
        msort.indata = data
        msort.outdata = outdata
        msort.sort = 'TB'
        if outdata.exists():
            print 'WARNING: Removing existing file ' + NAMERAW + '.' + sorted_class
            outdata.zap()
            print '         File removed. Proceeding with MSORT...'
        msort.go()
        indxr = task('indxr')
        indxr.default()
        data = UV(NAMERAW, sorted_class, 1, seq)
        indxr.indata = data
        indxr.cparm[3] = clint
        indxr.go()

########## EARTH ORIENTATION PARAMETERS CORRECTION ##########
# => CL2
if whattodo['eops_corr']:
    clinver = 1
    cloutver = 2
    data = UV(NAMERAW, sorted_class, DISK, 1)
    # Remove all CL tables higher than cloutver-1.
    for i in range(data.table_highver('CL'), cloutver-1, -1):
        data.zap_table('CL', i)
    data = UV(NAMERAW, sorted_class, DISK, 1)
    clcor = task('clcor')
    clcor.default()
    clcor.indata = data
    clcor.gainver = clinver
    clcor.gainuse = cloutver
    clcor.opcode = 'EOPS'
    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()


########## PREFLAGGING ##########
#Flag data known to be bad (whole antennas etc.)
if whattodo['preflagging']:
    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()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.timeran = [None, 0, 0, 0, 0, 0, 17, 0, 0]
    uvflg.antenna = [None, 5, 12, 14, 15]
    uvflg.sour = [None, ff]
    uvflg.reason = 'badscan'
    uvflg.opcode = 'FLAG'
    uvflg.go()
    
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.timeran = [None, 0, 0, 0, 0, 0, 18, 0, 0]
    uvflg.antenna = [None, 1, 13]
    uvflg.sour = [None, ff]
    uvflg.reason = 'badscan'
    uvflg.opcode = 'FLAG'
    uvflg.go()
    
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.timeran = [None, 1, 4, 0, 0, 2, 0, 0, 0]
    uvflg.antenna = [None, 6, 17]
    uvflg.sour = [None, ff]
    uvflg.reason = 'badscan'
    uvflg.opcode = 'FLAG'
    uvflg.go()
    
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.timeran = [None, 1, 1, 0, 0, 1, 2, 0, 0]
    uvflg.antenna = [None, 6]
    uvflg.sour = [None, ff]
    uvflg.reason = 'badscan'
    uvflg.opcode = 'FLAG'
    uvflg.go()
    
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.timeran = [None, 0, 0, 0, 0, 0, 20, 0, 0]
    uvflg.antenna = [None, 16]
    uvflg.sour = [None, ff]
    uvflg.reason = 'badscan'
    uvflg.opcode = 'FLAG'
    uvflg.go()
    
    # Quack
    quack = task('quack')
    quack.default()
    quack.indata = data
    quack.opcode = 'BEG'
    quack.aparm[2] = 30.0/60
    quack.go()
    quack = task('quack')
    quack.default()
    quack.indata = data
    quack.opcode = 'ENDB'
    quack.aparm[2] = 10.0/60
    #quack.go()
    
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.antenna = non_vlba_antennas
    uvflg.reason = 'NOCAL'
    uvflg.opcode = 'FLAG'
    #uvflg.go()
    
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.antenna = [None, 17]
    uvflg.sour = [None, ampcal]
    uvflg.reason = 'nofring'
    uvflg.opcode = 'FLAG'
    uvflg.go()
    
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.antenna = [None, 1]
    uvflg.sour = [None, ff]
    uvflg.reason = 'nofring'
    uvflg.timeran = [None, 0, 17, 0, 0, 0, 17, 14, 15]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.reason = 'badamp'
    uvflg.timeran = [None, 0, 16, 0, 0, 0, 16, 15, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    
    # Quack EB
    quack = task('quack')
    quack.default()
    quack.indata = data
    quack.opcode = 'ENDB'
    quack.aparm[2] = 20.0/60
    quack.antenna = [None, 1]
    quack.go()
    
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.antenna = [None,4,5]
    #uvflg.sour = [None, ampcal]
    uvflg.bchan = 1
    uvflg.echan = 2
    uvflg.reason = 'EDGE'
    uvflg.opcode = 'FLAG'
    uvflg.go()
    
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.timeran = [None, 1, 2, 52, 0, 1, 2, 56, 0]
    uvflg.bif = 2
    uvflg.eif = 2
    uvflg.bchan = 15
    uvflg.echan = 15
    uvflg.opcode = 'FLAG'
    #uvflg.go()
    
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.timeran = [None, 1, 2, 52, 0, 1, 2, 56, 0]
    uvflg.bif = 3
    uvflg.eif = 3
    uvflg.bchan = 2
    uvflg.echan = 2
    uvflg.opcode = 'FLAG'
    #uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 14]
    uvflg.timeran = [None, 0, 21, 35, 0, 0, 21, 35, 30]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 12]
    uvflg.timeran = [None, 0, 19, 29,0,0,19,30,30]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 13]
    uvflg.timeran = [None, 0, 23, 18, 30, 0, 23, 20, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 12]
    uvflg.timeran = [None, 0, 22, 0, 0, 1, 0, 0, 0]
    uvflg.opcode = 'FLAG'
    uvflg.bif = 4
    uvflg.eif = 4
    uvflg.bchan = 1
    uvflg.echan = 4
    uvflg.go()
    
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 12]
    uvflg.timeran = [None, 0, 18, 20, 30, 0, 18, 21, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 17]
    uvflg.timeran = [None, 0, 18, 20, 0, 0, 18, 21, 0]
    uvflg.opcode = 'FLAG'
    uvflg.bif = 2
    uvflg.eif = 2
    uvflg.bchan = 9
    uvflg.echan = 9
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 17]
    uvflg.timeran = [None, 1, 2, 15, 0, 1, 2, 22, 0]
    uvflg.opcode = 'FLAG'
    uvflg.bif = 1
    uvflg.eif = 1
    uvflg.bchan = 9
    uvflg.echan = 9
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 17]
    uvflg.timeran = [None, 1, 0, 30, 0, 1, 0, 31, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 5]
    uvflg.timeran = [None, 1, 2, 2, 0, 1, 2, 4, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 17]
    uvflg.timeran = [None, 1, 2, 44, 0, 1, 2, 56, 0]
    uvflg.opcode = 'FLAG'
    uvflg.bif = 2
    uvflg.eif = 2
    uvflg.bchan = 15
    uvflg.echan = 15
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 17]
    uvflg.timeran = [None, 1, 2, 44, 0, 1, 2, 56, 0]
    uvflg.opcode = 'FLAG'
    uvflg.bif = 3
    uvflg.eif = 3
    uvflg.bchan = 1
    uvflg.echan = 6
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 15]
    uvflg.timeran = [None, 1, 2, 43, 30, 1, 2, 44, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 6]
    uvflg.timeran = [None, 1, 2, 52, 30, 1, 2, 55, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 6]
    uvflg.timeran = [None, 1, 2, 24, 0, 1, 2, 30, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 6]
    uvflg.timeran = [None, 1, 2, 34, 0, 1, 2, 36, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 3]
    uvflg.timeran = [None, 1, 2, 52, 0, 1, 2, 54, 0]
    uvflg.opcode = 'FLAG'
    uvflg.bif = 3
    uvflg.eif = 3
    uvflg.bchan = 1
    uvflg.echan = 2
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 5]
    uvflg.timeran = [None, 1, 2, 15, 0, 1, 2, 16, 30]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 4]
    uvflg.timeran = [None, 1, 2, 19, 0, 1, 2, 21, 0]
    uvflg.opcode = 'FLAG'
    uvflg.bif = 4
    uvflg.eif = 4
    uvflg.bchan = 13
    uvflg.echan = 13
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 3]
    uvflg.timeran = [None, 1, 3, 16, 0, 1, 3, 22, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 2]
    uvflg.timeran = [None, 1, 3, 15, 0, 1, 3, 16, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 4]
    uvflg.timeran = [None, 1, 3, 28, 0, 1, 3, 29, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 2]
    # Flag all end since very high noise in this ant here
    uvflg.timeran = [None, 1, 3, 35, 0, 1, 5,0,0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 12]
    uvflg.timeran = [None, 0, 21, 35, 0, 0, 21, 35, 30]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 2]
    uvflg.timeran = [None, 0, 18, 42, 30, 0, 18, 44, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 17]
    uvflg.timeran = [None, 0, 18, 37, 0, 0, 18, 38, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 5]
    uvflg.timeran = [None, 0, 18, 9, 0, 0, 18, 10, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 13]
    uvflg.timeran = [None, 0,18, 7, 45, 0, 18, 8, 15]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 5]
    uvflg.timeran = [None, 0, 18, 15, 0, 0, 18, 16, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 15]
    uvflg.timeran = [None, 0, 18, 13, 15, 0, 18, 13, 45]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 12]
    uvflg.timeran = [None, 1, 4, 16, 0, 1, 4, 17, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 6]
    uvflg.timeran = [None, 1, 4, 18, 0, 1, 4, 18, 30]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 5]
    uvflg.timeran = [None, 1, 4, 19, 30, 1, 4, 20, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 17]
    uvflg.timeran = [None, 0, 16, 0, 0, 0, 16, 30, 0]
    uvflg.opcode = 'FLAG'
    uvflg.bif = 3
    uvflg.eif = 3
    uvflg.bchan = 5
    uvflg.echan = 5
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 17]
    uvflg.timeran = [None, 0, 16, 18, 0, 0, 16, 26, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 2]
    uvflg.timeran = [None, 0, 16, 43, 0, 0, 16, 44, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 13]
    uvflg.timeran = [None, 0, 16, 41, 0, 0, 16, 42, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 17]
    uvflg.timeran = [None, 0, 16, 55, 0, 0, 16, 56, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 16]
    uvflg.timeran = [None, 0, 17, 17, 0, 0, 17, 18, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 14]
    uvflg.timeran = [None, 0, 17, 18, 30, 0, 17, 19, 30]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 3]
    uvflg.timeran = [None, 0, 17, 27, 30, 0, 17, 28, 30]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 12]
    uvflg.timeran = [None, 0, 17, 30, 30, 0, 17, 31, 30]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 4]
    uvflg.timeran = [None, 0, 17, 45, 0, 0, 17, 46, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 3]
    uvflg.timeran = [None, 0, 17, 48, 0, 0, 17, 49, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 6]
    uvflg.timeran = [None, 0, 17, 47, 0, 0, 18, 5, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 6]
    uvflg.timeran = [None, 0, 21, 30, 0, 0, 22, 0, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 6]
    uvflg.timeran = [None, 1, 2, 35, 0, 1, 3, 15, 0]
    uvflg.opcode = 'FLAG'
    uvflg.go()



# FLAG antennas with NO AMPCAL. Can be solved and add these?
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = outfgver
    uvflg.sour = [None, target]
    uvflg.antenn = [None, 7,8,9,10,11,18] # Flag ants with no ampcal
    uvflg.opcode = 'FLAG'
    uvflg.go()
    
 
########## PARALLACTIC ANGLE CORRECTION ##########
if whattodo['par_ang_cor']:
    clinver = 2
    cloutver = 3
    data = UV(NAME, sorted_class, DISK, 1)
    # Remove all CL tables higher than 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)
    clcor = task('clcor')
    clcor.default()
    clcor.indata = data
    clcor.gainver = clinver
    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()

############# AMPLITUDE calibration ##############
if whattodo['apriori_ampcal']:
    # Define initial versionnumbers, increased further down
    clinver = 3
    cloutver = 4

    data = UV(NAME, sorted_class, DISK, 1)
    tacop = task('tacop')
    tacop.default()
    tacop.indata = data
    tacop.inext = 'CL'
    tacop.inver = clinver
    tacop.ncount = 1
    tacop.outdata = data
    tacop.outvers = cloutver
    
    snoutver = 1
    tyinver = 0
    tyoutver = 1
    gcinver = 0
    gcoutver = 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)
    # Remove all TY tables higher than version tyoutver-1.
    for i in range(data.table_highver('TY'), tyoutver-1, -1):
        data.zap_table('TY', i)
    # Remove all GC tables higher than version gcoutver-1.
    for i in range(data.table_highver('GC'), gcoutver-1, -1):
        data.zap_table('GC', i)
    
    antab = task('antab')
    data = UV(NAME, sorted_class, DISK, 1)
    antab.indata = data
    antab.calin = 'PWD:gl021acal.vlba.tsys'
    antab.tyver = tyoutver
    antab.blver = -1
    antab.offset = 0.5
    antab.go()
    
    antab = task('antab')
    data = UV(NAME, sorted_class, DISK, 1)
    antab.indata = data
    antab.calin = 'PWD:gl021acal.vlba.gain'
    antab.gcver = gcoutver
    antab.blver = -1
    antab.offset = 0.5
    antab.go()
    

    # Run APCAL 
    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)
    apcal = task('apcal')
    apcal.indata = data
    apcal.tyver = tyoutver
    apcal.gcver = gcoutver
    apcal.snver = snoutver
    apcal.calin = 'PWD:gl021cal.vlba.weather'
    apcal.solint = 2 
    apcal.opcode = ''
    apcal.sources = [None, '']
    apcal.antenna = vlba_antennas
    apcal.go()

    # Needed to calibrate VLA
    setjy = task('setjy')
    setjy.default()
    data = UV(NAME, sorted_class, DISK, 1)
    setjy.indata = data
    setjy.sources = [None, '']
    setjy.optype = 'REJY' # Reset
    setjy.go()
    
    setjy = task('setjy')
    setjy.default()
    data = UV(NAME, sorted_class, DISK, 1)
    setjy.indata = data
    setjy.sources = [None, ff]
    setjy.zerosp = [None, 3.8, 0, 0, 0] # From VLBA only clean image
    setjy.go()
    
    setjy = task('setjy')
    setjy.default()
    data = UV(NAME, sorted_class, DISK, 1)
    setjy.indata = data
    setjy.sources = [None, target, ampcal]
    setjy.zerosp = [None, 1.0, 0, 0, 0] # N/A for now, set as 1.0 according to APCAL help file
    setjy.go()
    
    antab = task('antab')
    data = UV(NAME, sorted_class, DISK, 1)
    antab.indata = data
    antab.calin = 'PWD:gl021acal.y'
    antab.gcver = gcoutver
    antab.tyver = tyoutver
    antab.blver = -1
    antab.offset = 0.5
    antab.go()
    
    # Run APCAL 
    data = UV(NAME, sorted_class, DISK, 1)
    apcal = task('apcal')
    apcal.indata = data
    apcal.tyver = tyoutver
    apcal.gcver = gcoutver
    apcal.snver = snoutver+1
    apcal.solint = 2 
    apcal.opcode = ''
    apcal.sources = [None, '']
    apcal.antenna = [None, 6]
    apcal.go()
    
    
    # Calibrate EB
    antab = task('antab')
    data = UV(NAME, sorted_class, DISK, 1)
    antab.indata = data
    antab.calin = 'PWD:gl021acal.eb.gain+tsys'
    antab.gcver = gcoutver
    antab.tyver = tyoutver
    antab.blver = -1
    antab.offset = 0.5
    antab.go()
    
    # Run APCAL 
    data = UV(NAME, sorted_class, DISK, 1)
    apcal = task('apcal')
    apcal.indata = data
    apcal.tyver = tyoutver
    apcal.gcver = gcoutver
    apcal.snver = snoutver+2
    apcal.solint = 2 
    apcal.opcode = ''
    apcal.sources = [None, '']
    apcal.antenna = [None, 1]
    apcal.go()
    
    setjy = task('setjy')
    setjy.default()
    data = UV(NAME, sorted_class, DISK, 1)
    setjy.indata = data
    setjy.sources = [None, '']
    setjy.optype = 'REJY' # Reset
    setjy.go()

    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)
    # Make CL table
    clcal = task('clcal') 
    clcal.sources = [None, '']
    clcal.indata = data
    clcal.timeran = [None, 0]
    clcal.samptype = 'BOX'
    clcal.bparm = [None, 0.00001]
    clcal.doblank = 1
    clcal.snver = snoutver
    clcal.invers = snoutver+2
    clcal.gainver = clinver
    clcal.gainuse = cloutver
    clcal.refant = refant
    clcal.interpol = 'SELF'
    clcal.antenna = vlba_antennas + [1, 6,]
    clcal.go()

if whattodo['fring_cals']:
    clinver = 4
    cloutver = 5
    snoutver = 4
    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, ff, ampcal]
    fring.refant = refant
    fring.solint = 0.5
    fring.snver = snoutver
    fring.dparm[2] = 500 # ns
    fring.dparm[3] = 50 #mHz
    #fring.aparm[7] = 3 # SNR
    fring.dparm[9] = 1 # Do not fit rate
    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, ff, ampcal]
    clcal.sour = [None, ''] # All sources
    clcal.gainver = clinver
    clcal.gainuse = cloutver
    clcal.refant = refant
    clcal.go()

if whattodo['bpass']:
    
    clinver = 5
    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.antennas = vlba_antennas + [6,1]
    bpass.calsour = [None,ff]
    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_ff_1']:
    clinver = 5
    imagename = FFCALIM + '1'
    imagedisk = 1
    imageseq = 1
    data = UV(NAME, sorted_class, DISK, 1)
    imagr = task('imagr')
    imagr.default()
    imagr.indata = data
    imagr.sources = [None, ff]
    imagr.gainuse = clinver
    imagr.flagver = 0
    imagr.nchav = 512
    imagr.stokes = 'I'
    imagr.cellsize = [None, 0.4e-3, 0.4e-3]
    imagr.imsize = [None, 1024, 1024]
    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.antenna = vlba_antennas
    imagr.baselin = imagr.antenna
    imagr.uvrange = uvrange
    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_ff']:
    clinver = 5
    cloutver = 6
    snoutver = 5
    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.nmaps = 1
    calib.in2data = IM(FFCALIM + '1', 'ICL001', 1, 1)
    calib.calsour = [None, ff]
    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[1] = 4 # Min no antenna
    calib.aparm[6] = 1 # Print level
    #calib.aparm[9] = 1 # Pass on failed sols, for scans without these ants
    calib.antenna = vlba_antennas
    calib.doband = 1
    calib.bpver = 1
    calib.uvrange = uvrange
    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.antenna = vlba_antennas
    clcal.calsour = [None, ff] # Use sols from Ampcal 
    clcal.sour = [None, ''] # Apply to all sources
    clcal.gainver = clinver
    clcal.gainuse = cloutver
    clcal.refant = refant
    clcal.go()

if whattodo['image_ff_2']:
    clinver = 6
    imagename = FFCALIM+'2'
    imagedisk = 1
    imageseq = 1
    data = UV(NAME, sorted_class, DISK, 1)
    imagr = task('imagr')
    imagr.default()
    imagr.indata = data
    imagr.sources = [None, ff]
    imagr.gainuse = clinver
    imagr.flagver = 0
    imagr.nchav = 512
    imagr.stokes = 'I'
    imagr.cellsize = [None, 0.4e-3, 0.4e-3]
    imagr.imsize = [None, 1024, 1024]
    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.antenna = vlba_antennas
    imagr.baselin = imagr.antenna
    imagr.uvrange = uvrange
    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_ff_2']:
    clinver = 6
    cloutver = 7
    snoutver = 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):
        data.zap_table('SN', i)
    data = UV(NAME, sorted_class, DISK, 1)
    calib = task('calib')
    calib.default()
    calib.indata = data
    calib.in2data = IM(FFCALIM + '2', 'ICL001', 1, 1)
    calib.nmaps = 1
    calib.calsour = [None, ff]
    calib.refant = refant
    calib.docalib = 1
    calib.gainuse = clinver
    calib.flagver = 0
    calib.snver = snoutver
    calib.solmode = 'A&P'
    calib.solint = 5 # 
    calib.aparm[1] = 4 # Min no antennas
    calib.aparm[7] = 5 # SNR
    calib.cmethod = 'DFT'
    calib.aparm[6] = 1 # Print level
    calib.antennas = vlba_antennas + [6,1]
    calib.dofit = [None, 6,1]
    calib.doband = 1
    calib.bpver = 1
    calib.timeran = [None, 0, 18, 0, 0, 0, 20, 0, 0]  # use good scan
    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.antenna = [None, 6,1]
    clcal.calsour = [None, ff] # Use sols from Ampcal 
    clcal.sour = [None, ''] # Apply to all sources
    clcal.gainver = clinver
    clcal.gainuse = cloutver
    clcal.refant = refant
    clcal.go()

if whattodo['image_ff_3']:
    clinver = 7
    imagename = FFCALIM+'3'
    imagedisk = 1
    imageseq = 1
    data = UV(NAME, sorted_class, DISK, 1)
    imagr = task('imagr')
    imagr.default()
    imagr.indata = data
    imagr.sources = [None, ff]
    imagr.gainuse = clinver
    imagr.flagver = 0
    imagr.nchav = 512
    imagr.stokes = 'I'
    imagr.cellsize = [None, 0.4e-3, 0.4e-3]
    imagr.imsize = [None, 1024, 1024]
    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.timeran = [None, 0, 0, 0, 0, 0, 18,0, 0] 
    #imagr.antenna = [None, -11]
    imagr.uvrange = uvrange
    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['image_ampcal_1']:
    clinver = 6
    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.4e-3, 0.4e-3]
    imagr.imsize = [None, 1024, 1024]
    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 = uvrange
    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 = 7
    cloutver = 8
    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)
    calib = task('calib')
    calib.default()
    calib.indata = data
    calib.in2data = IM(AMPCALIM + '1', 'ICL001', 1, 1)
    calib.nmaps = 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[1] = 4 # Min no antennas
    calib.aparm[7] = 5 # SNR
    calib.cmethod = 'DFT'
    calib.aparm[6] = 1 # Print level
    calib.doband = 1
    calib.bpver = 1
    calib.uvrange = uvrange
    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, ampcal]
    clcal.gainver = clinver
    clcal.gainuse = cloutver
    clcal.refant = refant
    clcal.go()

if whattodo['image_ampcal_2']:
    clinver = 8
    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.4e-3, 0.4e-3]
    imagr.imsize = [None, 1024, 1024]
    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 = uvrange
    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_target_maser']:
    clinver = 8
    cloutver = 9
    snoutver = 8
    bpinver = 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)
        print "Deleting SN "+str(i)
    data = UV(NAME, sorted_class, DISK, 1)
    calib = task('calib')
    calib.default()
    calib.indata = data
    calib.calsour = [None, target]
    calib.refant = refant
    calib.docalib = 1
    calib.gainuse = clinver
    calib.doband = 1
    calib.bpver = bpinver
    calib.snver = snoutver
    calib.ichansel[1] = [None, masercal_bchan, masercal_echan, 1, masercal_bif]
    calib.solmode = 'P'
    calib.soltype = 'L1R'
    calib.weightit = 1
    calib.solint = 2
    calib.aparm[1] = 3 # Min no of antennas
    calib.aparm[7] = 1 # SNR
    calib.aparm[3] = 1 # AVG RR,LL, 
    calib.aparm[6] = 1 # Print level
    calib.uvrange = uvrange
    calib.go()

    # replace contifs sols with maserif
    sncor = task('sncor')
    sncor.indata = data
    sncor.opcode = 'CPSN'
    sncor.sncorprm[1] = 1 # maser IF
    sncor.bif = 2
    sncor.eif = 4
    sncor.snver = snoutver
    sncor.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, target]
    clcal.sour = [None, target]
    clcal.gainver = clinver
    clcal.gainuse = cloutver
    clcal.refant = refant
    clcal.opcode = 'CALP'
    clcal.go()

if whattodo['image_ARP220_maser_1']:
    arpfilename = NAME + 'MAS1'
    arpfiledisk = 1
    arpfileseq = 1
    arpfileclass = 'ICL001'
    bpinver = 1
    clinver = 9
    # FIRST CLEAN WITH SMALL BOXES
    # Remove previous old images
    image = IM(arpfilename, arpfileclass, arpfiledisk, arpfileseq)
    if image.exists():
        image.zap()
    # Also remove old beam images
    image = IM(arpfilename, 'IBM001', arpfiledisk, arpfileseq)
    if image.exists():
        image.zap()
    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.
    imagr.gainuse = clinver
    imagr.flagver = 0
    imagr.doband = 1
    imagr.bpver = bpinver
    imagr.nchav = masercal_nchav
    imagr.stokes = 'I'
    imagr.nfield = 1
    imagr.docalib = 1
    imagr.dotv = -1
    imagr.uvrange = uvrange
    imagr.outname = arpfilename
    imagr.outdisk = arpfiledisk
    imagr.outseq = arpfileseq
    imagr.niter = 100
    imagr.cellsize = [None, 0.40e-3, 0.40e-3]
    imagr.imsize =[None, 1024,1024]  
    imagr.gain = 0.1
    imagr.bchan = masercal_bchan
    imagr.echan = masercal_echan
    imagr.bif = masercal_bif
    imagr.eif = masercal_eif
    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.doband = 1
    imagr.bpver = bpinver
    imagr.nchav = masercal_nchav
    imagr.stokes = 'I'
    imagr.nfield = 1
    imagr.docalib = 1
    imagr.dotv = -1
    imagr.uvrange = uvrange
    imagr.outname = arpfilename
    imagr.outdisk = arpfiledisk
    imagr.outseq = arpfileseq
    imagr.outver = 1
    imagr.niter = 300
    imagr.cellsize = [None, 0.40e-3, 0.40e-3]
    imagr.imsize =[None, 1024,1024]  
    imagr.gain = 0.1
    imagr.bcomp = [None, 100000, 100000]
    imagr.bchan = masercal_bchan
    imagr.echan = masercal_echan
    imagr.bif = masercal_bif
    imagr.eif = masercal_eif
    imagr.go()

if whattodo['selfcal_target_maser_2']:
    clinver = 8
    cloutver = 10
    bpinver = 1
    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):
        data.zap_table('SN', i)
        print "Deleting SN "+str(i)
    data = UV(NAME, sorted_class, DISK, 1)
    calib = task('calib')
    calib.default()
    calib.indata = data
    calib.calsour = [None, target]
    calib.refant = refant
    calib.docalib = 1
    calib.gainuse = clinver
    calib.doband = 1
    calib.bpver = bpinver
    calib.snver = snoutver
    calib.ichansel[1] = [None, masercal_bchan, masercal_echan, 1, masercal_bif]
    calib.solmode = 'P'
    calib.soltype = 'L1R'
    calib.weightit = 1
    calib.in2data = IM(NAME + 'MAS1', 'ICL001', 1, 1)
    calib.nmaps = 1
    calib.solint = 1
    calib.aparm[1] = 3 # Min no of antennas
    calib.aparm[7] = 1# SNR
    calib.aparm[3] = 1 # AVG RR,LL, 
    calib.aparm[6] = 1 # Print level
    calib.uvrange = uvrange
    calib.go()
    
    # replace if 2-4 sols with if1
    sncor = task('sncor')
    sncor.indata = data
    sncor.opcode = 'CPSN'
    sncor.sncorprm[1] = 1 # maser IF
    sncor.bif = 2
    sncor.eif = 4
    sncor.snver = snoutver
    sncor.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, target]
    clcal.sour = [None, target]
    clcal.gainver = clinver
    clcal.gainuse = cloutver
    clcal.refant = refant
    clcal.opcode = 'CALP'
    clcal.go()

if whattodo['image_ARP220_maser_2']:
    arpfilename = NAME + 'MAS2'
    arpfiledisk = 1
    arpfileseq = 1
    arpfileclass = 'ICL001'
    bpinver = 1
    clinver = 10
    # FIRST CLEAN WITH SMALL BOXES
    # Remove previous old images
    image = IM(arpfilename, arpfileclass, arpfiledisk, arpfileseq)
    if image.exists():
        image.zap()
    # Also remove old beam images
    image = IM(arpfilename, 'IBM001', arpfiledisk, arpfileseq)
    if image.exists():
        image.zap()
    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.
    imagr.gainuse = clinver
    imagr.flagver = 0
    imagr.doband = 1
    imagr.bpver = bpinver
    imagr.nchav = masercal_nchav
    imagr.stokes = 'I'
    imagr.nfield = 1
    imagr.docalib = 1
    imagr.dotv = -1
    imagr.uvrange = uvrange
    imagr.outname = arpfilename
    imagr.outdisk = arpfiledisk
    imagr.outseq = arpfileseq
    imagr.niter = 100
    imagr.cellsize = [None, 0.40e-3, 0.40e-3]
    imagr.imsize =[None, 1024,1024]  
    imagr.gain = 0.1
    imagr.bchan = masercal_bchan
    imagr.echan = masercal_echan
    imagr.bif = masercal_bif
    imagr.eif = masercal_eif
    imagr.antenna = vlba_antennas
    imagr.baselin = imagr.antenna
    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.doband = 1
    imagr.bpver = bpinver
    imagr.nchav = masercal_nchav
    imagr.stokes = 'I'
    imagr.nfield = 1
    imagr.docalib = 1
    imagr.dotv = -1
    imagr.uvrange = uvrange
    imagr.outname = arpfilename
    imagr.outdisk = arpfiledisk
    imagr.outseq = arpfileseq
    imagr.outver = 1
    imagr.niter = 300
    imagr.cellsize = [None, 0.40e-3, 0.40e-3]
    imagr.imsize =[None, 1024,1024]  
    imagr.gain = 0.1
    imagr.bcomp = [None, 100000, 100000]
    imagr.bchan = masercal_bchan
    imagr.echan = masercal_echan
    imagr.bif = masercal_bif
    imagr.eif = masercal_eif
    imagr.antenna = vlba_antennas
    imagr.baselin = imagr.antenna
    imagr.go()

if whattodo['selfcal_target_maser_3']:
    clinver = 10
    cloutver = 11
    bpinver = 1
    snoutver = 10

    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)
        print "Deleting SN "+str(i)
    data = UV(NAME, sorted_class, DISK, 1)
    calib = task('calib')
    calib.default()
    calib.indata = data
    calib.calsour = [None, target]
    calib.refant = refant
    calib.docalib = 1
    calib.gainuse = clinver
    calib.doband = 1
    calib.bpver = bpinver
    calib.snver = snoutver
    calib.ichansel[1] = [None, masercal_bchan, masercal_echan, 1, masercal_bif]
    calib.solmode = 'A&P'
    calib.soltype = 'L1R'
    calib.weightit = 1
    calib.in2data = IM(NAME + 'MAS2', 'ICL001', 1, 1)
    calib.nmaps = 1
    calib.solint = 10
    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.uvrange = uvrange
    calib.dofit = [None, 1, 6,13]
    calib.go()
    
    # replace if 2-4 sols with if1
    sncor = task('sncor')
    sncor.indata = data
    sncor.opcode = 'CPSN'
    sncor.sncorprm[1] = 1 # maser IF
    sncor.bif = 2
    sncor.eif = 4
    sncor.snver = snoutver
    sncor.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, target]
    clcal.sour = [None, target]
    clcal.gainver = clinver
    clcal.gainuse = cloutver
    clcal.refant = refant
    clcal.opcode = 'CALP'
    clcal.go()

if whattodo['image_ARP220_maser_3']:
    arpfilename = NAME + 'MAS3'
    arpfiledisk = 1
    arpfileseq = 1
    arpfileclass = 'ICL001'
    bpinver = 1
    clinver = 11
    # FIRST CLEAN WITH SMALL BOXES
    # Remove previous old images
    image = IM(arpfilename, arpfileclass, arpfiledisk, arpfileseq)
    if image.exists():
        image.zap()
    # Also remove old beam images
    image = IM(arpfilename, 'IBM001', arpfiledisk, arpfileseq)
    if image.exists():
        image.zap()
    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.
    imagr.gainuse = clinver
    imagr.flagver = 0
    imagr.doband = 1
    imagr.bpver = bpinver
    imagr.nchav = masercal_nchav
    imagr.stokes = 'I'
    imagr.nfield = 1
    imagr.docalib = 1
    imagr.dotv = -1
    imagr.uvrange = uvrange
    imagr.outname = arpfilename
    imagr.outdisk = arpfiledisk
    imagr.outseq = arpfileseq
    imagr.niter = 100
    imagr.cellsize = [None, 0.40e-3, 0.40e-3]
    imagr.imsize =[None, 1024,1024]  
    imagr.gain = 0.1
    imagr.bchan = masercal_bchan
    imagr.echan = masercal_echan
    imagr.bif = masercal_bif
    imagr.eif = masercal_eif
    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.doband = 1
    imagr.bpver = bpinver
    imagr.nchav = masercal_nchav
    imagr.stokes = 'I'
    imagr.nfield = 1
    imagr.docalib = 1
    imagr.dotv = -1
    imagr.uvrange = uvrange
    imagr.outname = arpfilename
    imagr.outdisk = arpfiledisk
    imagr.outseq = arpfileseq
    imagr.outver = 1
    imagr.niter = 400
    imagr.cellsize = [None, 0.40e-3, 0.40e-3]
    imagr.imsize =[None, 1024,1024]  
    imagr.gain = 0.1
    imagr.bcomp = [None, 100000, 100000]
    imagr.bchan = masercal_bchan
    imagr.echan = masercal_echan
    imagr.bif = masercal_bif
    imagr.eif = masercal_eif
    imagr.go()

if whattodo['split']:
    clinver = 11
    bpinver = 1
    fginver = 0
    
    data = UV(NAME, sorted_class, DISK, 1)
    split = task('split')
    split.default()
    split.indata = data
    split.docal = 1
    split.gainuse =clinver
    split.doband = 1
    split.bpver = 1
    split.flagver = fginver
    split.outclass = MCLASS
    split.sour = [None, target]
    tdata = UV(target, MCLASS, 1, 1)
    if tdata.exists():
        tdata.zap()
    split.go()

if whattodo['fixwt']:
    idata = UV(target, MCLASS, 1, 1)
    odata = UV(target, 'FIXWT', 1, 1)
    fixwt = task('fixwt')
    fixwt.default()
    fixwt.indata = idata
    fixwt.outdata = odata
    fixwt.solint = 5
    if odata.exists():
        odata.zap()
    fixwt.go()

if whattodo['multi']:
    idata = UV(target, 'FIXWT', 1, 1)
    odata = UV(target, MCLASS, 1, 1)
    multi = task('multi')
    multi.default()
    multi.indata = idata
    multi.outdata = odata
    if odata.exists():
        odata.zap()
    multi.go()

if whattodo['flag_maser']:
    # FLAG maser for final imaging
    data = UV(target, MCLASS, DISK, 1)
    uvflg = task('uvflg')
    uvflg.default()
    uvflg.indata = data
    uvflg.outfgver = 1
    uvflg.bif = masercal_bif
    uvflg.eif = masercal_eif
    uvflg.reason = 'Away with maser'
    uvflg.opcode = 'FLAG'
    uvflg.bchan = 10
    uvflg.echan = 16
    uvflg.go()

if whattodo['calculate_shifts']:
    # The maser self-calibration align the position to the wrong place, since
    # the maser was supposed to be in the wrong place. We fix this by changing
    # the coordinates of the output FITS images when saving to disk.  But, we
    # need to calculate the shifts from the "wrong" maser position here, to get
    # the images aligned pixelwise.
    # From the phase-referenced epoch BB297A we get the maser peak position as
    # R.A. 15h34m57s.22435
    # Dec. 23d30'11''.6644
    # which in decimal degrees, for AIPS, means
    maser_ra_bb297 = (15.0+34.0/60.0+57.22435/3600.0)*360.0/24 # degrees
    maser_dec_bb297 = 23+30.0/60.0+11.6644/3600.0 # degrees
    # Self-cal on maser will set the maser position to the source position in
    # the UVFITS header.
    # From AIPS memo 117, the UVFITS coordinates are defined as
    # RAAPP Source apparent equatorial position RA coordinate
    # DECAPP Source apparent equatorial position DEC coordinate
    # RAEPO Source J2000 equatorial position RA coordinate
    # DECPO Source J2000 equatorial position DEC coordinate
    # RAOBS Pointing right ascension of equinox
    # RADEC Pointing declination of equinox

    # We want to compare the RAEPO, DECEPO pair with the maser coordinates.
    # For the phased-referenced epochs, Arp 220 pointing is set between the
    # two nuclei at 15:34:57.2500  23:30:11.330, or in decimal deg
    # from BP129_C header:
    arp220_ra_bb297 = 233.73854166666666287711
    arp220_dec_bb297 = 23.50314722222221419656
    # Relative to these coordinates, we use the shifts 
    arp220_rashifts_bb297 = np.array([0.6, -0.35]) # in arcsec
    arp220_decshifts_bb297 = np.array([0.02, 0.18]) # in arcsec
    # to image the two Arp220 nuclei.
    # What are the corresponding shift for this epoch, i.e.
    # what would be the relative shifts with respect to the maser?

    # First, get the current maser coordinates
    wdata = WUV(target, MCLASS, 1, 1)
    sutab = wdata.table('SU',1)
    for row in sutab:
        if 'ARP220' in row.source:
            this_ra = row.raepo
            this_dec = row.decepo
    sutab.close()
    # this_ra, and this_dec are now the coordinates of the maser in this 
    # epoch. 
    # Then, calculate the difference between the actual maser position,
    # not the one in this epoch, and the point between the two nuclei
    #dmas_ra = maser_ra_bb297-this_ra
    #dmas_dec = maser_dec_bb297-this_dec
    dra = arp220_ra_bb297-maser_ra_bb297
    ddec = arp220_dec_bb297-maser_dec_bb297
    # now convert this to arcseconds, to have same unit as IMAGR wants
    # i.e. units of shifts given above
    dra_arcsec = dra*np.cos((np.pi/180.0)*arp220_dec_bb297)*3600 # scale with cos(DEC)
    ddec_arcsec = ddec*3600

    masercal_rashift = arp220_rashifts_bb297 + dra_arcsec
    masercal_decshift = arp220_decshifts_bb297 + ddec_arcsec
    # Convert to right type for imagr
    masercal_rashift = [float(i) for i in masercal_rashift]
    masercal_decshift = [float(i) for i in masercal_decshift]


############### 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 = 1
    # 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()
    data = UV(target, MCLASS, 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, 2048,1024]
    imagr.cellsize = [None, 0.5e-3, 0.5e-3]
    imagr.stokes = 'I'
    imagr.nfield = 2
    imagr.overlap = 1
    imagr.do3dimag = 1
    imagr.rashift = [None,] + masercal_rashift
    imagr.decshift = [None,] + masercal_decshift
    imagr.docalib = 1
    imagr.dotv = -1
    imagr.outname = arpfilename
    imagr.outdisk = arpfiledisk
    imagr.outseq = arpfileseq
    imagr.niter = 500
    imagr.gain = 0.1
    imagr.uvrange = uvrange # remove large-scale disturbances
    imagr.robust = 0.5
    imagr.bif = 1
    imagr.eif = 1
    imagr.go()
    # CLEAN WITH LARGE (NO) BOXES
    arpboxfile = ''
    data = UV(target, MCLASS, 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, 2048,1024]
    imagr.cellsize = [None, 0.5e-3, 0.5e-3]
    imagr.stokes = 'I'
    imagr.nfield = 2
    imagr.overlap = 1
    imagr.do3dimag = 1
    imagr.rashift = [None,] + masercal_rashift
    imagr.decshift = [None,] + masercal_decshift
    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 = uvrange # remove large-scale disturbances
    imagr.robust = 0.5
    imagr.bif = 1
    imagr.eif = 1
    imagr.go()

if whattodo['selfcal_arp220']:
    clinver = 1
    cloutver = 2
    snoutver = 1
    data = UV(target, MCLASS, 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(target, MCLASS, DISK, 1)
    calib = task('calib')
    calib.default()
    calib.indata = data
    # Do not use any model, just a point works fine for phasecal
    calib.in2data = IM(TCALIM, 'ICL001', 1, 1)
    calib.calsour = [None, target]
    calib.nmaps = 2
    calib.refant = 5
    calib.docalib = 1
    calib.gainuse = clinver
    calib.flagver = 0
    calib.snver = snoutver
    calib.solmode = 'P'
    calib.solint = 5
    #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 = 'L1'
    #calib.doband = 1
    #calib.bpver = 1
    calib.uvrange = uvrange
    calib.go()
    
    # Apply SN-table
    data = UV(target, MCLASS, 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(target, MCLASS, 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 = TCALIM+'V'
arpfiledisk = 1
arpfileseq = 1
arpfileclass = 'VCL001'
arpfileclass2 = 'VCL002'
if whattodo['image_ARP220_V']:
    clinver = 2
    # 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, 'VBM001', arpfiledisk, arpfileseq)
    if image.exists():
        image.zap()
    image = IM(arpfilename, 'VBM002', arpfiledisk, arpfileseq)
    if image.exists():
        image.zap()
    # CLEAN WITH LARGE (NO) BOXES
    arpboxfile = ''
    data = UV(target, MCLASS, 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, 2048,1024]
    imagr.cellsize = [None, 0.4e-3, 0.4e-3]
    imagr.stokes = 'V'
    imagr.nfield = 2
    imagr.overlap = 1
    imagr.do3dimag = 1
    imagr.rashift = [None,] + masercal_rashift
    imagr.decshift = [None,] + masercal_decshift
    imagr.docalib = 1
    imagr.dotv = -1
    imagr.outname = arpfilename
    imagr.outdisk = arpfiledisk
    imagr.outseq = arpfileseq
    imagr.outver = 1
    imagr.niter = 1
    imagr.gain = 0.1
    imagr.uvrange = uvrange # remove large-scale disturbances
    imagr.robust = 0.5
    #imagr.antenna = [None, -14] # Skip AR
    imagr.go()


############### IMAGE ARP220 ######################
# Load these outside if to be available also for export later
arpfilename = NAME + '_FIN'
arpfiledisk = 1
arpfileseq = 1
arpfileclass = 'ICL001'
arpfileclass2 = 'ICL002'
if whattodo['image_ARP220_final']:
    clinver = 2
    # 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(target, MCLASS, 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,] + masercal_rashift
    imagr.decshift = [None,] + masercal_decshift
    imagr.docalib = 1
    imagr.dotv = -1
    imagr.outname = arpfilename
    imagr.outdisk = arpfiledisk
    imagr.outseq = arpfileseq
    imagr.niter = 100
    imagr.gain = 0.1
    imagr.uvrange = uvrange # remove large-scale disturbances
    imagr.robust = 0.5
    imagr.go()
    # CLEAN WITH LARGE (NO) BOXES
    arpboxfile = ''
    data = UV(target, MCLASS, 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,] + masercal_rashift
    imagr.decshift = [None,] + masercal_decshift
    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 = uvrange # 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:'+NAME+'_EAST_IMAGE.FITS'
    if os.path.exists(fittp.dataout[4:]):
     os.system('rm ' + fittp.dataout[4:])
    fittp.go()
    # Change coordinates of file based on BB297A maser reference position
    f = fits.open(NAME+'_EAST_IMAGE.FITS', mode='update')
    hdr = f[0].header
    hdr['CRVAL1'] = arp220_ra_bb297+arp220_rashifts_bb297[0]/(3600.0*np.cos(arp220_dec_bb297*np.pi/180))
    hdr['CRVAL2'] = arp220_dec_bb297+arp220_decshifts_bb297[0]/(3600.0)
    f.flush()
    f.close()

    # Export WEST
    fittp = task('fittp')
    fittp.default()
    fittp.indata = west
    fittp.dataout = 'PWD:'+NAME+'_WEST_IMAGE.FITS'
    if os.path.exists(fittp.dataout[4:]):
     os.system('rm ' + fittp.dataout[4:])
    fittp.go()
    # Change coordinates of file based on BB297A maser reference position
    f = fits.open(NAME+'_WEST_IMAGE.FITS', mode='update')
    hdr = f[0].header
    hdr['CRVAL1'] = arp220_ra_bb297+arp220_rashifts_bb297[1]/(3600.0*np.cos(arp220_dec_bb297*np.pi/180))
    hdr['CRVAL2'] = arp220_dec_bb297+arp220_decshifts_bb297[1]/(3600.0)
    f.flush()
    f.close()
    
    # Export MASER
    masim = IM(NAME+'MAS3', 'ICL001', 1, 1)
    fittp = task('fittp')
    fittp.default()
    fittp.indata = masim
    fittp.dataout = 'PWD:'+NAME+'_MASER_IMAGE.FITS'
    if os.path.exists(fittp.dataout[4:]):
     os.system('rm ' + fittp.dataout[4:])
    fittp.go()
    # Change coordinates of file based on BB297A maser reference position
    f = fits.open(NAME+'_MASER_IMAGE.FITS', mode='update')
    hdr = f[0].header
    hdr['CRVAL1'] = maser_ra_bb297
    hdr['CRVAL2'] = maser_dec_bb297
    f.flush()
    f.close()
    
tac = time.time()
processtime = tac-tic
print 'OVERALL PROCESS TIME: ',processtime,' SECONDS. or ', processtime/60.0, 'minutes'