# http://www.vlba.nrao.edu/astro/VOBS/astronomy/nov94/ 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 vlba_antennas = [None, 1,3,6,8,9,10,11,12,13,15] non_vlba_antennas = [None, 2, 4, 5, 7, 14, 16, 17] #Ant 1 = BR BX= 144371.6364 BY= -2172972.1533 BZ= 879022.0771 #Ant 2 = EB BX= -3945899.8424 BY= 4041440.2472 BZ= 1052638.9812 #Ant 3 = FD BX= 1783990.2669 BY= -1411891.4223 BZ= -615829.2461 #Ant 4 = GB BX= 1412793.8372 BY= 801429.9203 BZ= 96338.9585 #Ant 5 = GO BX= 1076239.6366 BY= -2429953.6259 BZ= -170739.4810 #Ant 6 = HN BX= 945622.4768 BY= 1372718.5791 BZ= 474514.3832 #Ant 7 = JB BX= -3308381.3724 BY= 3819560.0241 BZ= 1238694.5802 #Ant 8 = KP BX= 1478073.3453 BY= -2078599.5143 BZ= -490463.5991 #Ant 9 = LA BX= 1425078.8647 BY= -1531723.4364 BZ= -138667.7956 #Ant 10 = MK BX= -1120930.3314 BY= -5504539.5901 BZ= -1699494.9633 #Ant 11 = NL BX= 1233908.2384 BY= -209505.7379 BZ= 379059.3167 #Ant 12 = OV BX= 912576.4151 BY= -2482786.7268 BZ= -9174.3426 #Ant 13 = PT BX= 1461433.1257 BY= -1723551.3511 BZ= -272379.8524 #Ant 14 = RO BX= -3085381.9092 BY= 4842482.8782 BZ= 267317.3068 #Ant 15 = SC BX= 2004792.8623 BY= 2516855.3813 BZ= -1915052.2200 #Ant 16 = WB BX= -3907535.2268 BY= 3835271.5619 BZ= 1217131.3984 #Ant 17 = Y BX= 1489247.9766 BY= -1684235.0340 BZ= -292915.5016 # Note regarding maser calibration: # Earth rotation is 0.3km/s. # At Lband this means maximum relative difference of the # antennas of (0.3km/s / speed of light) * 1650MHz=1.65Khz. # this is nothing compared to the continuum channel width of 125kHz. # For finding rates, phases and amps, we do not need CVEL. # Rest frame freqs of maser lines are 1665 and 1667 MHz, masercal_bchan = 14 masercal_echan = 14 masercal_bif = 2 masercal_eif = 2 masercal_nchav = 512 whattodo = {'load_data': True, # ->CL1 'msort+indxr': True, 'eops_corr': True, # ->CL2, before vbglu since need CT tables etc. 'split_vbglu_dbcon': True, 'msort_dbcon': True, #'copy_ct_table': True, 'delete_extrafiles': True, # NO ION CORR, NO DATA FROM 1994 'preflagging': True, 'par_ang_cor': True, # ->CL2 'apriori_ampcal' : True, # SN1,CL3 'fring_cals' : True, #->SN2,3 CL4 'selfcal_ff_1' : True, # -> SN4, CL5 # Correct phases of FF, target 'selfcal_ff_2' : True, # -> SN5, CL6 # correct phases of FF, target 'selfcal_ff_3' : True, # -> SN6, CL7 # Correct phases of FF 'setjy': True, # RESET SU table, bpass complains on velocities. 'bpass': True, # BP 1 'image_ff_1': True, 'selfcal_ff_4' : True, # -> SN7, CL8 # Correct amps of ff, target 'selfcal_ff_5' : True, # -> SN8, CL9 # Correct amps of ff, target 'selfcal_ff_6' : True, # -> SN9, CL10 # Correct amps of ff 'image_ff_2': True, 'split': True, 'fixwt': True, 'multi': True, 'unflag_maser': True, 'selfcal_target_maser' : True, # -> SN1, CL2. USING POINT MODEL. 'image_ARP220_maser_1': True, 'selfcal_target_maser_2' : True, # -> SN2, CL3 Using image model, PHASE ONLY, 1min 'image_ARP220_maser_2': True, 'selfcal_target_maser_3' : True, # -> SN3, CL4. AMP + PHASE. Using model 2, 5min. 'image_ARP220_maser_3': True, 'flag_maser':True, # If excluding the maser from final images. Flag whole IF to get rid of all maser chans. 'calculate_shifts': True, # Calculate shifts for Arp220 imaging. MUST be enabled to set variables used by imaging. 'image_ARP220': True, 'selfcal_arp220': True, # ->SN4, CL5 'image_ARP220_2': True, 'selfcal_arp220_2': True, # ->SN5, CL6 'image_ARP220_3': True, 'image_ARP220_V': False, ###### 'image_ARP220_final': True, 'export_results': True, } ########## Initialize observation data ########## AIPS.userno = 1000 clint = 10.0/60 # Seconds NAME = 'GL015_L' NAMERAW = 'GL015_L_RAW' CLASS = 'UVDATA' DISK = 1 AMPCALIM = NAME + '_AC' FFCALIM = NAME + '_FC' PCALIM = NAME + '_PC' TCALIM = NAME + '_TC' TCALIM2 = TCALIM + '2' OUTPREFIX = NAME target = 'ARP220' ampcal = '1300+580' ff = 'OQ208' uvrange = [None, 5000, 0] sorted_class = 'MSORT' dbcon_class = 'DBCON' logfile = 'log_'+OUTPREFIX+'_PTred.log' refant = 4 #GB #########^^^ 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/GL015/file.uvfits1' #fitld.outdata = outdata fitld.outname = outdata.name fitld.sour = [None, 'ARP220', '1300+580', 'OQ208', '3C84'] fitld.ncount = 97 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 ]: 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 for seq in [1,2,3,4,5]: data = UV(NAMERAW, sorted_class, DISK, seq) # 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, seq) 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() if whattodo['split_vbglu_dbcon']: clinver = 2 sources = ['ARP220', '1300+580', 'OQ208', '3C84'] lfreqs = [2, 2, 2, 1] # Split out lowdata data = UV(NAMERAW, sorted_class, 1, 3) # lower IFs of continuum for i, sour in enumerate(sources): split = task('split') split.indata = data split.sour = [None, sour] split.outclass = 'GL15L' split.freqid = lfreqs[i] split.docal = 1 split.gainuse = clinver outdata = UV(sour, split.outclass, 1, 1) if outdata.exists(): print 'WARNING: Removing existing file ' + outdata.name + '.' + outdata.klass outdata.zap() print ' File removed. Proceeding with FITLD...' split.go() # Split out highdata data = UV(NAMERAW, sorted_class, 1, 5) # higher IFs of continuum hfreqs = [1, 1, 1, 2] for i, sour in enumerate(sources): split = task('split') split.indata = data split.sour = [None, sour] split.outclass = 'GL15H' split.freqid = hfreqs[i] split.docal = 1 split.gainuse = clinver outdata = UV(sour, split.outclass, 1, 1) if outdata.exists(): print 'WARNING: Removing existing file ' + outdata.name + '.' + outdata.klass outdata.zap() print ' File removed. Proceeding with FITLD...' split.go() # Combine all sources together in all-IF files for sour in sources: vbglu = task('vbglu') vbglu.indata = UV(sour, 'GL15L', 1, 1) vbglu.in2data = UV(sour, 'GL15H', 1, 1) outdata = UV(sour, 'GL15LH', 1, 1) vbglu.outdata = outdata if outdata.exists(): print 'WARNING: Removing existing file ' + outdata.name + '.' + outdata.klass outdata.zap() vbglu.go() multi = task('multi') multi.indata = UV(sour, 'GL15LH', 1, 1) outdata = UV(sour, 'GL15M', 1, 1) multi.outdata = outdata if outdata.exists(): print 'WARNING: Removing existing file ' + outdata.name + '.' + outdata.klass outdata.zap() multi.go() #dbcontmp1 = UV('3C84+OQ208', 'GL15M', 1, 1) #dbcon = task('dbcon') #dbcon.indata = UV('3C84', 'GL15M', 1, 1) #dbcon.in2data = UV('OQ208', 'GL15M', 1, 1) #dbcon.outdata = dbcontmp1 #if dbcontmp1.exists(): # print 'WARNING: Removing existing file ' + dbcontmp1.name + '.' + dbcontmp1.klass # dbcontmp1.zap() # print ' File removed. Proceeding with FITLD...' #dbcon.go() dbcontmp = UV('OQ13', 'GL15M', 1, 1) dbcon2 = task('dbcon') #dbcon2.indata = dbcontmp1 dbcon2.indata = UV('OQ208', 'GL15M', 1, 1) dbcon2.in2data = UV('1300+580', 'GL15M', 1, 1) dbcon2.outdata = dbcontmp if dbcontmp.exists(): print 'WARNING: Removing existing file ' + dbcontmp.name + '.' + dbcontmp.klass dbcontmp.zap() dbcon2.go() dbconfin = UV(NAMERAW, dbcon_class, 1, 1) dbcon3 = task('dbcon') dbcon3.indata = dbcontmp dbcon3.in2data = UV('ARP220', 'GL15M', 1, 1) dbcon3.outdata = dbconfin if dbconfin.exists(): print 'WARNING: Removing existing file ' + dbconfin.name + '.' + dbconfin.klass dbconfin.zap() dbcon3.go() if dbcontmp.exists(): print 'WARNING: Removing existing file ' + dbcontmp.name + '.' + dbcontmp.klass dbcontmp.zap() if whattodo['msort_dbcon']: msort = task('msort') data = UV(NAMERAW, dbcon_class, 1, 1) outdata = UV(NAME, sorted_class, 1, 1) msort.default() msort.indata = data msort.outdata = outdata msort.sort = 'TB' if outdata.exists(): print 'WARNING: Removing existing file ' + outdata.name + '.' + outdata.klass outdata.zap() print ' File removed. Proceeding with MSORT...' msort.go() indxr = task('indxr') indxr.default() data = UV(NAME, sorted_class, 1, 1) indxr.indata = data indxr.cparm[3] = clint indxr.go() #if whattodo['copy_ct_table']: # outdata = UV(NAME, sorted_class, DISK, 1) # inver = 0 # outver = 1 # # Remove all CT tables higher than infgver. # for i in range(outdata.table_highver('CT'), inver, -1): # outdata.zap_table('CT', i) # indata = UV(NAMERAW, sorted_class, DISK, 4)# use continuum file with CT table # # Copy CT table from raw to dbconned sorted file # tacop = task('tacop') # tacop.default() # tacop.indata = indata # tacop.inext = 'CT' # tacop.inver = inver # tacop.ncount = 1 # tacop.outdata = outdata # tacop.outvers = outver # tacop.go() if whattodo['delete_extrafiles']: sources = ['ARP220', '1300+580', 'OQ208', '3C84'] # Clean up files for sour in sources: outdata = UV(sour, 'GL15L', 1, 1) if outdata.exists(): print 'WARNING: Removing existing file ' + outdata.name + '.' + outdata.klass outdata.zap() outdata = UV(sour, 'GL15H', 1, 1) if outdata.exists(): print 'WARNING: Removing existing file ' + outdata.name + '.' + outdata.klass outdata.zap() outdata = UV(sour, 'GL15LH', 1, 1) if outdata.exists(): print 'WARNING: Removing existing file ' + outdata.name + '.' + outdata.klass outdata.zap() outdata = UV(sour, 'GL15M', 1, 1) if outdata.exists(): print 'WARNING: Removing existing file ' + outdata.name + '.' + outdata.klass outdata.zap() # Delete all raw files, now when we have sorted and indexed versions. for seq in [1,2,3,4,5]: data = UV(NAMERAW, CLASS, 1, seq) if data.exists(): print 'WARNING: Removing existing file ' + data.name + '.' +data.klass + '.' + str(data.seq) data.zap() print ' file removed.' # Delete sorted files now when we have dbconned versions. for seq in [1,2,3,4,5]: data = UV(NAMERAW, sorted_class, 1, seq) if data.exists(): print 'WARNING: Removing existing file ' + data.name + '.' +data.klass + '.' + str(data.seq) data.zap() print ' file removed.' # Delete DBCON file, now when we have sorted version outdata = UV(NAMERAW, 'DBCON', 1, 1) if outdata.exists(): print 'WARNING: Removing existing file ' + outdata.name + '.' + outdata.klass outdata.zap() # POSSM complains that CQ table is missing. But, according to the explain section of # of the FXVLB task, the corrections should be applied when splitting after loading the data, # since there is a CQ table generated by FITLD in the pre-VBGLU files. So this is # probably not needed. #if whattodo['add_cq_table']: # data = UV(NAME, sorted_class, DISK, 1) # fxvlb = task('fxvlb') # fxvlb.indata = data # fxvlb.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, 17, 50, 0, 0, 18, 4, 0] uvflg.antenna = [None, 13] uvflg.sour = [None, target] uvflg.bif = 7 uvflg.eif = 7 uvflg.reason = 'badamp' uvflg.opcode = 'FLAG' uvflg.go() uvflg = task('uvflg') uvflg.default() uvflg.indata = data uvflg.outfgver = outfgver uvflg.timeran = [None, 0, 14, 50, 0, 0, 14, 55, 0] uvflg.antenna = [None, 12] uvflg.sour = [None, target] uvflg.reason = 'badamp' uvflg.opcode = 'FLAG' uvflg.go() uvflg = task('uvflg') uvflg.default() uvflg.indata = data uvflg.outfgver = outfgver uvflg.timeran = [None, 0, 0, 0, 0, 0, 14, 0, 0] uvflg.antenna = [None, 12] uvflg.sour = [None, target] uvflg.reason = 'badamp' uvflg.opcode = 'FLAG' uvflg.go() uvflg = task('uvflg') uvflg.default() uvflg.indata = data uvflg.outfgver = outfgver uvflg.timeran = [None, 0, 13, 25, 0, 0, 13, 30, 0] uvflg.antenna = [None, 1,12] uvflg.sour = [None, target] uvflg.reason = 'badamp' uvflg.opcode = 'FLAG' uvflg.go() uvflg = task('uvflg') uvflg.default() uvflg.indata = data uvflg.outfgver = outfgver uvflg.timeran = [None, 0, 16, 46, 0, 0, 16, 48, 0] uvflg.antenna = [None, 9,3] uvflg.sour = [None, target] uvflg.reason = 'badamp' uvflg.opcode = 'FLAG' uvflg.go() uvflg = task('uvflg') uvflg.default() uvflg.indata = data uvflg.outfgver = outfgver uvflg.timeran = [None, 0, 13, 3, 0, 0, 13, 5, 0] uvflg.antenna = [None, 9,12] uvflg.sour = [None, target] uvflg.reason = 'badamp' uvflg.opcode = 'FLAG' uvflg.go() uvflg = task('uvflg') uvflg.default() uvflg.indata = data uvflg.outfgver = outfgver uvflg.timeran = [None, 0, 19, 27,0, 0, 19, 30, 0] uvflg.antenna = [None, 12] uvflg.sour = [None, target] uvflg.reason = 'badamp' uvflg.opcode = 'FLAG' uvflg.go() uvflg = task('uvflg') uvflg.default() uvflg.indata = data uvflg.outfgver = outfgver uvflg.timeran = [None, 0, 19, 27,0, 0, 19, 30, 0] uvflg.antenna = [None, 17] uvflg.bif = 10 uvflg.eif = 12 uvflg.sour = [None, target] uvflg.reason = 'badamp' uvflg.opcode = 'FLAG' uvflg.go() uvflg = task('uvflg') uvflg.default() uvflg.indata = data uvflg.outfgver = outfgver uvflg.timeran = [None, 0, 19, 20, 0, 0, 19, 50, 0] uvflg.antenna = [None, 15] uvflg.sour = [None, target] uvflg.reason = 'badamp' uvflg.opcode = 'FLAG' uvflg.bif = 4 uvflg.eif = 4 uvflg.go() uvflg = task('uvflg') uvflg.default() uvflg.indata = data uvflg.outfgver = outfgver uvflg.timeran = [None, 0, 14, 36, 0, 0, 14, 40, 0] uvflg.antenna = [None, 12] uvflg.sour = [None, target] uvflg.reason = 'badamp' uvflg.opcode = 'FLAG' uvflg.go() uvflg = task('uvflg') uvflg.default() uvflg.indata = data uvflg.outfgver = outfgver uvflg.timeran = [None, 0, 19, 28, 0, 0, 19, 40, 0] uvflg.antenna = [None, 1] uvflg.bif = 13 uvflg.eif = 13 uvflg.sour = [None, target] uvflg.reason = 'badamp' uvflg.opcode = 'FLAG' uvflg.go() uvflg = task('uvflg') uvflg.default() uvflg.indata = data uvflg.outfgver = outfgver uvflg.timeran = [None, 0, 13, 35, 0, 0, 13, 50, 0] uvflg.antenna = [None, 5] uvflg.bif = 11 uvflg.eif = 11 uvflg.sour = [None, target] uvflg.reason = 'badamp' uvflg.opcode = 'FLAG' uvflg.go() uvflg = task('uvflg') uvflg.default() uvflg.indata = data uvflg.outfgver = outfgver uvflg.timeran = [None, 0, 16, 20, 0, 0, 16, 22, 0] uvflg.antenna = [None, 1,3] uvflg.sour = [None, target] uvflg.reason = 'badamp' uvflg.opcode = 'FLAG' uvflg.go() uvflg = task('uvflg') uvflg.default() uvflg.indata = data uvflg.outfgver = outfgver uvflg.timeran = [None, 0, 13, 38, 0, 0, 13, 42, 0] uvflg.antenna = [None, 3] uvflg.sour = [None, target] uvflg.reason = 'badamp' uvflg.opcode = 'FLAG' uvflg.go() uvflg = task('uvflg') uvflg.default() uvflg.indata = data uvflg.outfgver = outfgver uvflg.timeran = [None, 0, 18, 0, 0, 1, 0, 0, 0] uvflg.antenna = [None, 14] uvflg.bif = 9 uvflg.eif = 9 uvflg.sour = [None, ff] uvflg.reason = 'badamp' 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, 11] uvflg.baselin = [None, 2] uvflg.bif = 9 uvflg.eif = 14 uvflg.sour = [None, ff] uvflg.reason = 'badamp' uvflg.opcode = 'FLAG' uvflg.go() uvflg = task('uvflg') uvflg.default() uvflg.indata = data uvflg.outfgver = outfgver uvflg.timeran = [None, 0, 18, 0, 0, 0, 18, 40, 14] uvflg.antenna = [None, 15] uvflg.sour = [None, ff] uvflg.reason = 'badamp' uvflg.opcode = 'FLAG' uvflg.go() uvflg = task('uvflg') uvflg.default() uvflg.indata = data uvflg.outfgver = outfgver uvflg.reason = 'EDGE' uvflg.opcode = 'FLAG' uvflg.bchan = 1 uvflg.echan = 2 uvflg.go() uvflg = task('uvflg') uvflg.default() uvflg.indata = data uvflg.outfgver = outfgver uvflg.reason = 'EDGE' uvflg.opcode = 'FLAG' uvflg.bchan = 15 uvflg.echan = 16 #uvflg.go() # Quack quack = task('quack') quack.default() quack.indata = data quack.opcode = 'BEG' quack.aparm[2] = 20.0/60 quack.go() quack = task('quack') quack.default() quack.indata = data quack.opcode = 'ENDB' quack.aparm[2] = 10.0/60 quack.go() quack = task('quack') quack.default() quack.indata = data quack.antenna = [None, 14] quack.bif = 9 quack.eif = 14 quack.opcode = 'BEG' quack.aparm[2] = 40.0/60 quack.go() uvflg.default() uvflg.indata = data uvflg.outfgver = outfgver uvflg.timeran = [None, 0, 18, 40, 0, 0, 18, 40, 10] uvflg.sour = [None, ff] uvflg.reason = 'badamp' uvflg.opcode = 'FLAG' uvflg.go() uvflg.default() uvflg.indata = data uvflg.outfgver = outfgver uvflg.timeran = [None, 0, 13, 0, 0, 0, 14, 0, 0] uvflg.sour = [None, ff] uvflg.antenna = [None, 17] uvflg.reason = 'badamp' uvflg.opcode = 'FLAG' uvflg.bif = 1 uvflg.eif = 8 uvflg.go() uvflg = task('uvflg') uvflg.default() uvflg.indata = data uvflg.outfgver = outfgver uvflg.timeran = [None, 0, 13, 4, 20, 0, 13, 4, 30] uvflg.sour = [None, target] uvflg.reason = 'badamp' uvflg.opcode = 'FLAG' uvflg.go() uvflg = task('uvflg') uvflg.default() uvflg.indata = data uvflg.outfgver = outfgver uvflg.timeran = [None, 0, 14, 37, 20, 0, 14, 37, 40] uvflg.sour = [None, target] uvflg.reason = 'badamp' uvflg.opcode = 'FLAG' uvflg.go() uvflg = task('uvflg') uvflg.default() uvflg.indata = data uvflg.outfgver = outfgver uvflg.timeran = [None, 0, 15, 51, 45, 0, 15, 52, 15] uvflg.sour = [None, target] uvflg.reason = 'badamp' uvflg.opcode = 'FLAG' uvflg.go() uvflg = task('uvflg') uvflg.default() uvflg.indata = data uvflg.outfgver = outfgver uvflg.timeran = [None, 0, 15, 4, 30, 0, 15, 5, 30] uvflg.sour = [None, target] uvflg.reason = 'badamp' uvflg.opcode = 'FLAG' uvflg.go() uvflg = task('uvflg') uvflg.default() uvflg.indata = data uvflg.outfgver = outfgver uvflg.timeran = [None, 0, 0, 0, 0, 0, 15, 0, 0] uvflg.sour = [None, ff] uvflg.reason = 'badamp' uvflg.opcode = 'FLAG' uvflg.antenna = [None, 10,17] uvflg.bif = 1 uvflg.eif = 8 uvflg.go() uvflg = task('uvflg') uvflg.default() uvflg.indata = data uvflg.outfgver = outfgver uvflg.timeran = [None,0, 18, 0, 0, 0, 20, 0, 0] uvflg.sour = [None, ff] uvflg.reason = 'badamp' uvflg.opcode = 'FLAG' uvflg.antenna = [None, 10,17] uvflg.bif = 9 uvflg.eif = 14 uvflg.go() uvflg = task('uvflg') uvflg.default() uvflg.indata = data uvflg.outfgver = outfgver uvflg.timeran = [None,0, 0, 0, 0, 0, 15, 0, 0] uvflg.sour = [None, ff] uvflg.reason = 'badamp' uvflg.opcode = 'FLAG' uvflg.antenna = [None, 5,15,16] uvflg.go() ########## PARALLACTIC ANGLE CORRECTION ########## if whattodo['par_ang_cor']: clinver = 1 cloutver = 2 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 = 2 cloutver = 3 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:gl015cal.vlba.edit.cont.gain' antab.tyver = tyoutver 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.invers = 1 # WX table version apcal.calin = 'PWD:gl015cal.vlba.weather' apcal.solint = 2 apcal.opcode = '' apcal.freqid = 1 apcal.sources = [None, ''] apcal.antenna = vlba_antennas apcal.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 clcal.gainver = clinver clcal.gainuse = cloutver clcal.refant = refant clcal.interpol = 'SELF' clcal.antenna = vlba_antennas clcal.go() if whattodo['fring_cals']: clinver = 3 cloutver = 4 snoutver = 2 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 = 5 fring.snver = snoutver fring.dparm[2] = 500 # ns #fring.dparm[3] = 50 #mHz #fring.aparm[7] = 3 f# SNR fring.dparm[8] = 5 # Zero rates and phases fring.dparm[9] = 1 # Do not fit rate fring.go() snsmo = task('snsmo') snsmo.indata = data snsmo.bparm = [None, 99,99,99,99,99,99,99,99,99,99] snsmo.samptype = 'BOX' snsmo.smotype = 'VLBI' snsmo.outvers = snoutver+1 snsmo.invers = snoutver 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, ff, ampcal] clcal.sour = [None, ''] # All sources clcal.gainver = clinver clcal.gainuse = cloutver clcal.refant = refant #clcal.interp = 'CALP' clcal.go() if whattodo['selfcal_ff_1']: 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) calib = task('calib') calib.default() calib.indata = data calib.calsour = [None, ff] calib.refant = refant calib.docalib = 1 calib.gainuse = clinver calib.flagver = 0 calib.snver = snoutver calib.solmode = 'P' calib.solint = 5 # calib.aparm[7] = 5 # SNR calib.aparm[1] = 4 # Min no antenna calib.aparm[6] = 1 # Print level calib.aparm[9] = 1 # Pass on failed to not flag antennas with no data in this scan calib.timeran = [None, 0, 0, 0, 0, 0, 15, 0, 0] calib.antenna = [None, 1,2,3,4,6,7,8,9,10,11,12,13,14,17] 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, ff] # Use sols from Ampcal clcal.sour = [None, ff, target] # Apply to all sources clcal.gainver = clinver clcal.gainuse = cloutver clcal.refant = refant clcal.antenna = calib.antenna clcal.go() if whattodo['selfcal_ff_2']: 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.calsour = [None, ff] calib.refant = refant calib.docalib = 1 calib.gainuse = clinver calib.flagver = 0 calib.snver = snoutver calib.solmode = 'P' calib.solint = 5 # calib.aparm[7] = 5 # SNR calib.aparm[1] = 3 # Min no antenna calib.aparm[6] = 1 # Print level calib.aparm[9] = 1 # Pass on failed to not flag antennas with no data in this scan calib.timeran = [None, 0, 18, 0, 0, 0, 20, 0, 0] calib.dofit= [None, 5,10,15,16,17] 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, ff] # Use sols from Ampcal clcal.sour = [None, ff, target] # Apply to all sources clcal.gainver = clinver clcal.gainuse = cloutver clcal.refant = refant clcal.antenna = calib.dofit clcal.go() if whattodo['selfcal_ff_3']: 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.calsour = [None, ff] calib.refant = refant calib.docalib = 1 calib.gainuse = clinver calib.flagver = 0 calib.snver = snoutver calib.solmode = '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.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] # Use sols from Ampcal clcal.sour = [None, ff] # Apply ONLY TO FF, not to cause failed interpolation on target clcal.gainver = clinver clcal.gainuse = cloutver clcal.refant = refant clcal.go() if whattodo['setjy']: # Needed to make bpass run setjy = task('setjy') setjy.default() data = UV(NAME, sorted_class, DISK, 1) setjy.indata = data setjy.sources = [None, ''] setjy.optype = 'RESE' # Reset setjy.go() setjy = task('setjy') setjy.default() data = UV(NAME, sorted_class, DISK, 1) setjy.indata = data setjy.sources = [None, ''] setjy.optype = 'VCAL' # Bpass complains about need restfreqs for all IFs. I just set them to the observed freqs, hope this means no shift. freqs = list(np.array([1.63562500, 1.63750000, 1.63962500, 1.64150000, 1.64762500, 1.64950000, 1.65562500, 1.65750000, 1.66362500, 1.66550000, 1.67162500, 1.67350000, 1.68782500, 1.68970000])*1e9) for i in range(1, 15): setjy.bif = i setjy.eif = i setjy.restfreq = [None, float(freqs[i-1])] setjy.go() if whattodo['bpass']: clinver = 7 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,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 = 7 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 = [None, 1,3,6,8,9,10,11,12,13] # Exclude antenna 15 SC with high amps imagr.antenna = vlba_antennas imagr.baselin = imagr.antenna imagr.bif = 1 imagr.eif = 8 # Exclude upper IFs with bad amps #imagr.uvrange = [None, 0, 43000] # Exclude bad longest baseline 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_4']: 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.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 to not flag antennas with no data in this scan calib.timeran = [None, 0, 0, 0, 0, 0, 15, 0, 0] calib.antenna = [None, 1,2,3,4,6,7,8,9,10,11,12,13,14,17] 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, ff] # Use sols from Ampcal clcal.sour = [None, ff, target] # Apply to all sources clcal.gainver = clinver clcal.gainuse = cloutver clcal.refant = refant clcal.antenna = calib.antenna clcal.go() if whattodo['selfcal_ff_5']: clinver = 8 cloutver = 9 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) 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 to not flag antennas with no data in this scan calib.timeran = [None, 0, 18, 0, 0, 0, 20, 0, 0] calib.dofit= [None, 5,10,15,16,17] 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, ff] # Use sols from Ampcal clcal.sour = [None, ff, target] # Apply to all sources clcal.gainver = clinver clcal.gainuse = cloutver clcal.refant = refant clcal.antenna = calib.dofit clcal.go() if whattodo['selfcal_ff_6']: clinver = 9 cloutver = 10 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) 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.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, ff] # Use sols from Ampcal clcal.sour = [None, ff] clcal.gainver = clinver clcal.gainuse = cloutver clcal.refant = refant clcal.go() if whattodo['image_ff_2']: clinver = 10 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 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['split']: clinver = 10 bpinver = 1 fginver = 2 data = UV(NAME, sorted_class, DISK, 1) split = task('split') split.default() split.sour = [None, target] split.indata = data split.docal = 1 split.gainuse =clinver split.doband = 1 split.bpver = 1 split.flagver = fginver split.outclass = 'SP' #for s in [target, ff, ampcal]: for s in [target, ]: tdata = UV(s, 'SP', 1, 1) if tdata.exists(): tdata.zap() split.go() if whattodo['fixwt']: for s in [target, ]: idata = UV(s, 'SP', 1, 1) odata = UV(s, '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']: for s in [target, ]: idata = UV(target, 'FIXWT', 1, 1) odata = UV(NAME, 'MULTI', 1, 1) multi = task('multi') multi.default() multi.indata = idata multi.outdata = odata if odata.exists(): odata.zap() multi.go() indxr = task('indxr') indxr.default() indxr.indata = odata indxr.cparm[3] = clint indxr.go() if whattodo['unflag_maser']: outfgver = 1 infgver = 0 data = UV(NAME, 'MULTI', DISK, 1) # Remove all FG tables higher than infgver. for i in range(data.table_highver('FG'), infgver, -1): data.zap_table('FG', i) if whattodo['selfcal_target_maser']: clinver = 1 cloutver = 2 snoutver = 1 data = UV(NAME, 'MULTI', 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, 'MULTI', DISK, 1) calib = task('calib') calib.default() calib.indata = data calib.calsour = [None, target] calib.refant = refant calib.docalib = 1 calib.gainuse = clinver calib.snver = snoutver calib.ichansel[1] = [None, masercal_bchan, masercal_echan, 1, masercal_bif] calib.solmode = 'P' calib.soltype = 'L1' calib.weightit = 1 calib.solint = 1 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.go() # replace contifs sols with maserif sncor = task('sncor') sncor.indata = data sncor.opcode = 'CPSN' sncor.sncorprm[1] = 2 # maser IF sncor.bif = 3 sncor.eif = 14 sncor.snver = snoutver sncor.go() # replace contifs sols with maserif sncor = task('sncor') sncor.indata = data sncor.opcode = 'CPSN' sncor.sncorprm[1] = 2 # maser IF sncor.bif = 1 sncor.eif = 1 sncor.snver = snoutver sncor.go() # Apply SN-table data = UV(NAME, 'MULTI', 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, 'MULTI', 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' clinver = 2 # 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, 'MULTI', 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.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, 'MULTI', DISK, 1) imagr = task('imagr') imagr.default() imagr.indata = data imagr.sources = [None, target] imagr.gainuse = clinver imagr.flagver = 0 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 = 2 cloutver = 3 snoutver = 2 data = UV(NAME, 'MULTI', 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, 'MULTI', DISK, 1) calib = task('calib') calib.default() calib.indata = data calib.calsour = [None, target] calib.refant = refant calib.docalib = 1 calib.gainuse = clinver 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 contifs sols with maserif sncor = task('sncor') sncor.indata = data sncor.opcode = 'CPSN' sncor.sncorprm[1] = 2 # maser IF sncor.bif = 3 sncor.eif = 14 sncor.snver = snoutver sncor.go() # replace contifs sols with maserif sncor = task('sncor') sncor.indata = data sncor.opcode = 'CPSN' sncor.sncorprm[1] = 2 # maser IF sncor.bif = 1 sncor.eif = 1 sncor.snver = snoutver sncor.go() # Apply SN-table data = UV(NAME, 'MULTI', 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, 'MULTI', 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' clinver = 3 # 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, 'MULTI', 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.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, 'MULTI', DISK, 1) imagr = task('imagr') imagr.default() imagr.indata = data imagr.sources = [None, target] imagr.gainuse = clinver imagr.flagver = 0 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 = 500 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_3']: clinver = 3 cloutver = 4 snoutver = 3 data = UV(NAME, 'MULTI', 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, 'MULTI', DISK, 1) calib = task('calib') calib.default() calib.indata = data calib.calsour = [None, target] calib.refant = refant calib.docalib = 1 calib.gainuse = clinver 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 = 15 calib.aparm[1] = 4 # Min no of antennas calib.aparm[7] = 3# 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] = 2 # maser IF sncor.bif = 3 sncor.eif = 14 sncor.snver = snoutver sncor.go() # replace contifs sols with maserif sncor = task('sncor') sncor.indata = data sncor.opcode = 'CPSN' sncor.sncorprm[1] = 2 # maser IF sncor.bif = 1 sncor.eif = 1 sncor.snver = snoutver sncor.go() # Apply SN-table data = UV(NAME, 'MULTI', 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, 'MULTI', 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' clinver = 4 # 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, 'MULTI', 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.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, 'MULTI', DISK, 1) imagr = task('imagr') imagr.default() imagr.indata = data imagr.sources = [None, target] imagr.gainuse = clinver imagr.flagver = 0 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['flag_maser']: outfgver = 1 infgver = 0 data = UV(NAME, 'MULTI', DISK, 1) # 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, 'MULTI', DISK, 1) # FLAG maser for final imaging uvflg = task('uvflg') uvflg.default() uvflg.indata = data uvflg.outfgver = outfgver uvflg.bif = masercal_bif uvflg.eif = masercal_eif uvflg.reason = 'Away with maser' uvflg.opcode = 'FLAG' #uvflg.bchan = 13 #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(NAME, 'MULTI', 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 = 4 # 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(NAME, 'MULTI', 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,2048] 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 = 100 imagr.gain = 0.1 imagr.uvrange = uvrange # remove large-scale disturbances imagr.robust = 0.5 imagr.bif = 3 imagr.eif = 3 imagr.go() # CLEAN WITH LARGE (NO) BOXES arpboxfile = '' data = UV(NAME, 'MULTI', 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,2048] 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 = 500 imagr.gain = 0.1 imagr.bcomp = [None, 100000, 100000] imagr.uvrange = uvrange # remove large-scale disturbances imagr.robust = 0.5 imagr.bif = 3 imagr.eif = 3 imagr.go() if whattodo['selfcal_arp220']: clinver = 4 cloutver = 5 snoutver = 4 data = UV(NAME, 'MULTI', 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, 'MULTI', 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 = refant calib.docalib = 1 calib.gainuse = clinver calib.flagver = 0 calib.snver = snoutver calib.solmode = 'P' calib.solint = 7 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.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(NAME, 'MULTI', 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, 'MULTI', 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+'2' arpfiledisk = 1 arpfileseq = 1 arpfileclass = 'ICL001' arpfileclass2 = 'ICL002' if whattodo['image_ARP220_2']: clinver = 5 # 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(NAME, 'MULTI', 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,2048] 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 = 100 imagr.gain = 0.1 imagr.uvrange = uvrange # remove large-scale disturbances imagr.robust = 0.5 #imagr.bif = 2 #imagr.eif = 2 imagr.go() # CLEAN WITH LARGE (NO) BOXES arpboxfile = '' data = UV(NAME, 'MULTI', 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,2048] 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 = 200 imagr.gain = 0.1 imagr.bcomp = [None, 100000, 100000] imagr.uvrange = uvrange # remove large-scale disturbances imagr.robust = 0.5 #imagr.bif = 2 #imagr.eif = 2 imagr.go() if whattodo['selfcal_arp220_2']: clinver = 4 cloutver = 6 snoutver = 5 data = UV(NAME, 'MULTI', 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, 'MULTI', DISK, 1) calib = task('calib') calib.default() calib.indata = data calib.in2data = IM(TCALIM+'2', '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 = 7 calib.aparm[1] = 3 # min No. of antennas to include first 45 min. calib.aparm[7] = 1 # 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(NAME, 'MULTI', 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, 'MULTI', 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+'3' arpfiledisk = 1 arpfileseq = 1 arpfileclass = 'ICL001' arpfileclass2 = 'ICL002' if whattodo['image_ARP220_3']: clinver = 6 # 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(NAME, 'MULTI', 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,2048] 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 = 100 imagr.gain = 0.1 imagr.uvrange = uvrange # remove large-scale disturbances imagr.robust = 0.5 #imagr.bif = 2 #imagr.eif = 2 imagr.go() # CLEAN WITH LARGE (NO) BOXES arpboxfile = '' data = UV(NAME, 'MULTI', 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,2048] 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 = 500 imagr.gain = 0.1 imagr.bcomp = [None, 100000, 100000] imagr.uvrange = uvrange # remove large-scale disturbances imagr.robust = 0.5 #imagr.bif = 2 #imagr.eif = 2 imagr.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 = 6 # 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(NAME, 'MULTI', 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,2048] imagr.cellsize = [None, 0.5e-3, 0.5e-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 = 2 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 = 6 # 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, 'MULTI', 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.bif = 1 #imagr.eif = 8 imagr.go() # CLEAN WITH LARGE (NO) BOXES arpboxfile = '' data = UV(NAME, 'MULTI', 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 = 200 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'