# n05l2pt.py # script and notes for reducing N05L2 using ParselTongue # 2005 Nov 25 James M Anderson --JIVE start # 2005 Dec 23 JMA --change name to n05l2pt because python can't handle .'s # well # This script is intended to produce the basic calibration necessary for # looking at N05L2. ################################################################################ # USER DEFINED STUFF HERE ################################################################################ ################################################################################ # set some runlevel information, for debugging purposes. # Set to 15 to rerun Phase Calibration Tests # Set to 9 to rerun from Ionospheric Calibration point. restart_runlevel = 15 # Should the program pause after plots to make sure the user has seen them? pause_after_plots = 1 # Should the program pause before overwriting datafiles? pause_before_overwrite = 0 ################################################################################ # some variables to set for Cormac's pipeline stuff # Also set a bunch of Global variables for my own use. experiment = 'n05l2' passwd = '' userno = 52 date = '050602' # observing date fits = '/jop30_0/anderson/astro/ionosphere/observations/2005May/N05L2/' nfits = 6 # number of IDI files to combine processing_option = 'IO_COD' #None # The ionospheric processing option name # set to None for no ionosphering processing # Working options # None # 'IO_COD' # 'IO_COR' # 'IO_ESA' # 'IO_ESR' # 'IO_IGR' # 'IO_IGS' # 'IO_JPL' # 'IO_JPR' # 'IO_UPC' # 'IO_UPR' refantname = 'EF' solint = float(30.) # main SOLINT in seconds !!! t_int = float(0.5) # correlator integration time in seconds !!! cl_table_interval = 2. # desired CL table interval in seconds !!! # By default, the AIPS Version defaults to 'NEW'. # Here at JIVE, that means 31DEC04, which is not very new at all. # In order to ensure that this script will work, I have hardcoded # date version. You may feel free to try replasing these with # 'TST' or 'NEW', etc. as you desire, but check the results. # Where a specific date has been left in an AIPSTask call, then # there is a good chance it will only work with that version. aips_version_try = '31DEC05' aips_version_std = '31DEC05' aips_version_old = '31DEC04' #AIPSTask.version = 'NEW' imsize = 512 #How big should the images be made? kntr_imsize = 256 # How big should the kntr plots be? # this stuff is used to find the instrumental phase in the absence of PCAL info calibrator_instrument = '3C273B' # # 0/16:47:40 - 0/16:55:39 calibrator_instrument_timer = [0, 16, 47, 40, 0, 16, 55, 39] # put all FRING antennas here, with the desired refant first calibrator_instrument_search = ['EF', 'WB', 'JB', 'MC', 'TR', 'NT'] calibrator_instrument_solint = float(60) # in seconds !!! # for n05l2 does not work below 20 # This is where the main calibrator names should be stored. If it is None, # then all objects are assumed to be calibrators, and FRING will be run on # all of them fring_calibrators = None # after user flagging, what is the best FG version. Should be set to None or 0 # to get the last version or when you are just starting out with processing. final_flagver = 8 # How many self-cal iterations Num_SelfCal_Iterations = 9 # Should image plots be dumped to disk in color or greyscale? Do_Image_PLots_In_Color = 1 ################################################################################ # The user should not need to modify things below here. ################################################################################ ################################################################################ # some import commands The User should not need to change this. ################################################################################ from AIPS import AIPS, AIPSDisk from AIPSTask import AIPSTask, AIPSList from AIPSData import AIPSUVData, AIPSImage from AIPSTV import AIPSTV import Wizardry.AIPSData import LocalProxy from xmlrpclib import ServerProxy import copy, optparse, os, sys import re, string import numarray.ma import numarray.ieeespecial import numarray.nd_image import inspect import warnings warnings.defaultaction = "always" import math ################################################################################ # Set up the reminaing global variable information ################################################################################ AIPS.userno = userno AIPSTask.isbatch = 0 # AIPS behaves differently in batch mode, so # even though this is a "batch" mode script, don't # tell that to AIPS. If you really want batch mode, # set this to 32000 (stupid AIPS magic numbers). ################################################################################ # set up some path information dirname = os.path.dirname(fits) if not dirname: dirname = os.getcwd() ################################################################################ # open the log file AIPS.log = open(dirname + '/' + experiment + '.log', "a") ################################################################################ # replace the warnings module's showwarning function with my own to print # to the AIPS.log file def showwarning(message, category, filename, lineno, file=None): """Hook to write a warning to AIPS.log file""" if file is None: file = AIPS.log try: file.write(warnings.formatwarning(message, category, filename, lineno)) file.flush() if(file != sys.stderr): sys.stderr.write(warnings.formatwarning(message, category, filename, lineno)) sys.stderr.flush() except IOError: pass # the file (probably stderr) is invalid - this warning gets lost. # Set the warnings function to use the AIPS log warnings.showwarning = showwarning fits_file = experiment + '_1_1' + '.IDI' # declare a global uvdataset to work on uvdata = AIPSUVData(experiment.upper(), 'UVDATA', 1, 1) # what value of the CL gain table should be globally used Gainuse_Global = 1 # What is the CL version of the non-FRINGed stuff? Gainuse_Before_Fring = 1 # what value of the bandpass table should be globally used Bandpass_Global = 0 # set a stupid refant refant = 1 # Now check that we have a valid flag table version if(type(final_flagver) != type(0)): final_flagver = 0 ################################################################################ # I want to have the possibility of using the TV tv_local = AIPSTV() ################################################################################ # Things which will be filled in properly in the calibrate_data function ################################################################################ uvflg_file = 'uvflg' chflg_file = 'chflag' antab_file = 'antab' ################################################################################ # Define some useful functions ################################################################################ ################################################################################ # 6&*!%#@^ Python doesn't have the equivalent of __LINE__ def current_line(): return inspect.getouterframes(inspect.currentframe())[1][2] ################################################################################ # Wait for the user to say that it is ok to move on def pause_prog(text): print '#' * 78 junk = raw_input(text) ################################################################################ # Find the antenna number def get_antenna_number(ant_name): ant = 0 for antenna in uvdata.antennas: ant = ant + 1 if antenna == ant_name: return ant raise RuntimeError, "antenna not found!" return None ################################################################################ # Find the antenna dictionary def get_antenna_dictionary(): ant = 0 a_dict = {} for antenna in uvdata.antennas: ant = ant + 1 a_dict[antenna] = ant return a_dict ################################################################################ # dump a dictionary to a file in text mode def print_dict_to_file(d, middle_name, overwrite=1, copy_to_stdout=0): """prints a dict to a file based on a middle_name""" # generate the full file name filename = dirname + '/' + experiment + '.' + middle_name + '.dict' if(os.path.isfile(filename)): if(overwrite): warnings.warn("File '%s' already exists. Deleting."%filename) os.remove(filename) try: fp = open(filename, "w") dk = d.keys() dk.sort() for k in dk: print >> fp, "%30s"%k, d[k] if(copy_to_stdout): print "%30s"%k, d[k] finally: fp.close() ################################################################################ # write out the last plot to disk def lwpla_last_plot(data_file, middle_name, overwrite=1): """plots out the last plot to disk""" # generate the full file name filename = dirname + '/' + experiment + '.' + middle_name + '.ps' if(os.path.isfile(filename)): if(overwrite): warnings.warn("File '%s' already exists. Deleting."%filename) os.remove(filename) else: warnings.warn("File '%s' already exists. AIPS will probably append to it."%filename) # now run lwpla lwpla = AIPSTask('lwpla') lwpla.indata = data_file lwpla.outfile = filename lwpla.dparm[6] = 4 lwpla.plver = 0 lwpla.inver = data_file.table_highver('AIPS PL') if(Do_Image_PLots_In_Color): lwpla.ofmfile = 'RAINBOW' lwpla.dparm[9] = 1 lwpla.dodark = 1 lwpla.docolor=1 lwpla.plcolors = [None, # 0 # Bright Lines [None, 1, 1, 0], # 1 Border lines, tick marks, internal character labels [None, 1, 1, 1], # 2 Contours, model lines [None, 1, 0.5, 1], # 3 Polarization vectors [None, 0, 1, 1], # 4 Stars (incl labels), symbols # Dark Lines [None, 0,0,0], # 5 internal character labels [None, 0,0,0], # 6 Contours, model lines [None, 0,0,0], # 7 Polarization vectors [None, 0,0,0], # 8 Stars (incl labels), symbols # Other Stuff [None, 0,0,0], # 9 Character labels outside the plot area [None, 0.8, 0.8, 1]] # 10 Overall background lwpla() ################################################################################ # a tool to help with snplt def snplot(inext, invers, optype, opcode='', stokes_use='', sources_use='', timerange=[0,0,0,0], snuvdata = uvdata, nplots=10): snplt = AIPSTask('snplt', version = aips_version_std) snplt.indata = snuvdata snplt.inext = inext snplt.invers = invers if(type(sources_use) == type('string')): snplt.sources[1] = sources_use else: snplt.sources[1:] = sources_use if(timerange != None): if(timerange[0] == None): snplt.timerang = timerange else: snplt.timerang[1:] = timerange snplt.optype = optype snplt.opcode = opcode snplt.nplots = nplots #snplt.factor = 0.5 snplt.cutoff = 1e-6 snplt.do3col = 1 snplt.dotv = 1 if(snplt.dotv > 0): tv_local.clear() snplt() if(pause_after_plots): pause_prog('Done with plotting. Continue?') ################################################################################ # something to help run possm def possm_1(sources_use, stokes_use, gainuse, solint): possm = AIPSTask('possm', version = aips_version_std) possm.indata = uvdata possm.solint = solint / 60.0 # convert from s to min if len(uvdata.stokes) == 1: possm.stokes = stokes else: possm.stokes = stokes_use possm.freqid = 1 if(gainuse >= 0): possm.docalib = 2 possm.gainuse = gainuse if(type(sources_use) == type('string')): possm.sources[1] = sources_use else: possm.sources[1:] = sources_use #possm.baseline[1:3] = [refant,0] #possm.timerang[1:] = [0, 0, 0, 1] possm.aparm[2] = 1 possm.aparm[5] = -180 possm.aparm[6] = 180 #possm.aparm[7] = 1 possm.aparm[8] = 1 possm.aparm[9] = 1 possm.codetype = 'A&P' possm.nplots = 9 possm.dotv = 1 if(possm.dotv > 0): tv_local.clear() possm() if(pause_after_plots): pause_prog('Done with plotting. Continue?') ################################################################################ def possm_2(sources_use, stokes_use, gainuse, solint, refant_use, flagver=0, bpver = 0): possm = AIPSTask('possm', version = aips_version_std) possm.indata = uvdata possm.solint = solint / 60.0 # convert from s to min possm.flagver = flagver if len(uvdata.stokes) == 1: possm.stokes = stokes else: possm.stokes = stokes_use possm.freqid = 1 if(gainuse >= 0): possm.docalib = 2 possm.gainuse = gainuse if(type(sources_use) == type('string')): possm.sources[1] = sources_use else: possm.sources[1:] = sources_use if(bpver > 0): possm.doband = 1 possm.bpver = bpver possm.baseline[1] = refant_use #possm.timerang[1:] = [0, 0, 0, 1] possm.aparm[2] = 1 possm.aparm[5] = -180 possm.aparm[6] = 180 #possm.aparm[7] = 1 #possm.aparm[8] = 1 possm.aparm[9] = 1 possm.codetype = 'A&P' possm.nplots = 9 possm.dotv = 1 if(possm.dotv > 0): tv_local.clear() possm() if(pause_after_plots): pause_prog('Done with plotting. Continue?') ################################################################################ # run the EDITR task to flag stuff def editr1(sources_use='',gainuse=0,flagver=0,timerange=[0,0,0,0],uvdata_use=uvdata): # print out the telescopes for antenna in uvdata.antennas: print antenna_dict[antenna], antenna # create and editr object editr = AIPSTask('editr', version = aips_version_std) editr.indata = uvdata_use if(type(sources_use) == type('string')): editr.sources[1] = sources_use editr.reason = 'Manual EDITR ' + sources_use else: editr.sources[1:] = sources_use editr.reason = 'Manual EDITR many' + sources_use if(timerange != None): if(timerange[0] == None): editr.timerang = timerange else: editr.timerang[1:] = timerange editr.docalib = 2 editr.gainuse = gainuse editr.flagver = flagver editr.stokes = '' editr.dohist = 1 editr.crowded = 1 editr.antuse = [None, 1] editr.antuse[2:] = range(1,len(uvdata.antennas)+1) #editr.inputs() editr() ################################################################################ # my kntr runner def run_kntr(imgdata, rms, blc=[0,0], trc=[0,0],dotv=1,pixrange=None): kntr = AIPSTask('kntr', version = aips_version_std) kntr.indata = imgdata kntr.docont = 1 kntr.dovect = 0 kntr.dogrey = 1 kntr.blc[1:] = blc kntr.trc[1:] = trc kntr.ltype = 4 kntr.clev = rms for i in range(1,17): kntr.levs[i+4] = 2**i for i in range(1,5): kntr.levs[5-i] = -(2**i) kntr.cbplot = 1 if(pixrange): if(pixrange[0] == None): kntr.pixrange = pixrange else: kntr.pixrange[1:] = pixrange if(Do_Image_PLots_In_Color): kntr.functype = 'LN' kntr.ofmfile = 'RAINBOW' kntr.darkline = 0 kntr.dowedge = 3 if(pixrange): kntr.dowedge = 1 kntr.darkline = 0 kntr.dotv = dotv if(kntr.dotv > 0): tv_local.clear() kntr() ################################################################################ # my kntr runner which takes in a C position def run_kntr2(imgdata, rms, C_center = None, width=None, dotv=1, pixrange=None): if(C_center == None): C_center = numarray.array([imgdata.header.naxis[1],imgdata.header.naxis[0]],type=numarray.Float32) C_center = C_Center * 0.5 - 1.0 if(width == None): width = max(imgdata.header.naxis[1],imgdata.header.naxis[0]) # get the C values blc = C_center - width/2 trc = C_center + width/2 if(blc[0] < 0): blc[0] = 0 if(blc[1] < 0): blc[1] = 0 if(trc[0] >= imgdata.header.naxis[0]): trc[0] = imgdata.header.naxis[0] -1 if(trc[1] >= imgdata.header.naxis[1]): trc[1] = imgdata.header.naxis[1] -1 # convert to FORTRAN blc = blc[-1:None:-1] + 1 trc = trc[-1:None:-1] + 1 # convert to simple lists blc = [blc[0],blc[1]] trc = [trc[0],trc[1]] return run_kntr(imgdata, rms, blc, trc, dotv, pixrange) ################################################################################ # get image statistics, assuming that the source is at the center def get_imstats(imgdata, blc=[0,0], trc=[0,0]): imean = AIPSTask('imean', version = aips_version_std) imean.indata = imgdata imean.blc[1:] = blc imean.trc[1:] = trc imean() return (imean.pixavg, imean.pixstd) ################################################################################ # bounds check aips_image def bounds_check_aips(aips_image, blc, trc): # This function assumes that blc is actually the bottom left corner # blc and trc are assumed to be FORTRAN lists bad_message = 'fully out of bounds' if(blc[1] < 1): if(trc[1] < 1): raise Exception, bad_message elif(trc[1] > aips_image.header['naxis'][0]): trc[1] = aips_image.header['naxis'][0] blc[1] = 1 elif(blc[1] > aips_image.header['naxis'][0]): raise Exception, bad_message if(blc[2] < 1): if(trc[2] < 1): raise Exception, bad_message elif(trc[2] > aips_image.header['naxis'][1]): trc[2] = header['naxis'][1] blc[2] = 1 elif(blc[2] > aips_image.header['naxis'][1]): raise Exception, bad_message ################################################################################ # run jmfit for a single component def jmfit_1(imgdata, outfile, center=None, half_width = 20): """run JMFIT for a single Gaussian component center is a FORTRAN array/list """ # JMFIT can only handle so many pixels if(half_width > 49): half_width = 49 jmfit = AIPSTask('jmfit', version = aips_version_std) jmfit.indata = imgdata if(center == None): center = [imgdata.header['naxis'][0]/2,imgdata.header['naxis'][1]/2] else: jmfit.gpos[1] = [None, center[0], center[1]] jmfit.gwidth[1] = [None, 1.5*imgdata.header.bmaj/imgdata.header.cdelt[0], 1.5*imgdata.header.bmin/imgdata.header.cdelt[0], imgdata.header.bpa] blc=[None,center[0]-half_width,center[1]-half_width] trc=[None,center[0]+half_width,center[1]+half_width] bounds_check_aips(imgdata,blc,trc) jmfit.blc = blc jmfit.trc = trc jmfit.ngauss = 1 jmfit.niter = 120 jmfit.domax[1] = 1 jmfit.dopos[1] = [None, 1, 1] jmfit.dowidth[1] = [None, 1, 1, 1] jmfit.docrt = -1 jmfit.outprint = outfile jmfit.prtlev=1 #jmfit.inp() jmfit() return (jmfit.fmax[1], jmfit.fpos[1][1:], jmfit.fwidth[1][1:]) ################################################################################ # use IMAGR to make an image def run_imagr(uvdata, cellsize, imsize, sources_use='', out_data=None, gainuse=0, niter=1000,timerange=[0,0,0,0],flagver=0,bpver=0,stokes='I',flux=0,dotv=1,clbox=None,robust_use=+8): # First, check what is going on with out_data if(out_data == None): out_data = AIPSImage(uvdata.name, 'ICL001', uvdata.disk, 1) # Next, make the image. Then, deal with any residual crud # set a flag so that I can tell easily whether or not IMAGR barfs retval=0 try: print 'trying to IMAGR', sources_use, 'with robust', robust_use print 'using dataset', uvdata imagr = AIPSTask('imagr', version = aips_version_std) imagr.indata = uvdata imagr.outname = out_data.name imagr.outseq = out_data.seq imagr.outdisk = out_data.disk if(type(sources_use) == type('string')): imagr.sources[1] = sources_use else: imagr.sources[1:] = sources_use if(timerange != None): if(timerange[0] == None): imagr.timerang = timerange else: imagr.timerang[1:] = timerange imagr.freqid = 1 if(gainuse >= 1): imagr.docalib = 2 imagr.gainuse = gainuse imagr.flagver=flagver if(bpver > 0): imagr.doband = 3 imagr.bpver = bpver imagr.nchav = uvdata.header['naxis'][2] #imagr.chinc = 1 imagr.stokes = stokes imagr.uvwtfn = 'U' imagr.robust = robust_use imagr.cellsize[1:] = [cellsize, cellsize] imagr.imsize[1:] = [imsize, imsize] imagr.nfield=0 imagr.xtype = 5 imagr.ytype = 5 imagr.niter = niter imagr.gain = 0.03 imagr.flux=flux if(clbox != None): imagr.nboxes = 1 if(clbox[0] != None): imagr.clbox[1] = AIPSList(clbox) else: imagr.clbox[1] = clbox imagr.overlap=1 imagr.minpatch = max(imsize / 3,101) #imagr.imagrprm[8] = -0.1e-4 imagr.dotv = dotv if(imagr.dotv > 0): tv_local.clear() #imagr.inp() retval = imagr() finally: # clean up after ourselves beam = AIPSImage(out_data.name, 'IBM001', out_data.disk, out_data.seq) if beam.exists(): beam.zap() # Did IMAGR barf? If no, then retval will be None if(retval != None): # yes, a barf if(robust_use > -4): # try again, with a smaller robust level robust_use = robust_use - 2 if(robust_use > 4): robust_use = 4 warnings.warn('IMAGR barfed on us, so clean up and try again with ROBUST =%.1f'% robust_use) if(out_data.exists()): out_data.zap() return run_imagr(uvdata, cellsize, imsize, sources_use, out_data, gainuse, niter,timerange,flagver,bpver,stokes,flux,dotv,clbox,robust_use) else: # robust level too low warnings.warn('IMAGR barfed on us, but I cannot make robust smaller than %.1f'% robust_use) # send stuff back to the user if we made an image if out_data.exists(): return out_data return None ################################################################################ # print out a CC file def print_cc(imgdata, middle_name): # print out the CC table outfile = dirname +'/'+ experiment +'.'+ middle_name + '.ccdat' if os.path.isfile(outfile): warnings.warn("File '%s' already exists. Deleting."%outfile) os.remove(outfile) prtcc = AIPSTask('prtcc', version = aips_version_std) prtcc.indata = imgdata prtcc.docrt = -1 prtcc.outprint = outfile prtcc() ################################################################################ # use CALIB to self-cal def run_calib(uvdata, refant, solint, solsub, sources_use='', gainuse=0,timerange=[0,0,0,0],flagver=0,bpver=0, image_data=None, ncomp=None, flux=None, smodel=None, aparm=None, soltype='', solmode='P!A', snver=1, save_data=0, out_data=None): """Runs CALIB to perform Self-Calibration""" if(out_data == None): out_data = AIPSUVData(uvdata.name, 'CALIB', uvdata.disk, 1) if out_data.exists(): out_data.zap() try: print 'trying to CALIB', sources_use, 'in dataset', uvdata calib = AIPSTask('calib', version = aips_version_std) calib.indata = uvdata calib.outdata = out_data if(type(sources_use) == type('string')): calib.calsour[1] = sources_use else: calib.calsour[1:] = sources_use if(timerange != None): if(timerange[0] == None): calib.timerang = timerange else: calib.timerang[1:] = timerange calib.freqid = 1 if(gainuse >= 1): calib.docalib = 2 calib.gainuse = gainuse calib.flagver=flagver if(bpver > 0): calib.doband = 3 calib.bpver = bpver calib.refant = refant calib.solint= solint / 60.0 # convert from s to min calib.solsub = solsub calib.soltype = soltype calib.solmode = solmode if(aparm != None): if(aparm[0] == None): calib.aparm = aparm else: calib.aparm[1:] = aparm if (type(image_data) != type(None)): # We have an image, hopefully with clean components, to use calib.in2data = image_data if (ncomp != None): if(type(ncomp) == type(1)): calib.ncomp[1] = ncomp else: calib.ncomp[1:] = ncomp calib.nmaps = len(ncomp) if(flux != None): calib.flux = flux else: # no image, so smodel needs to be checked if(smodel != None): if((type(smodel) == type(1)) or (type(smodel) == type(1.0))): calib.smodel[1] = smodel else: calib.smodel[1:] = smodel calib.snver = snver #calib.inp() calib() except: # clean up after ourselves if out_data.exists(): out_data.zap() raise if(save_data): if out_data.exists(): return out_data return None ################################################################################ # use FRING to self-cal def run_fring(uvdata, refant, search_list, solint, solsub, sources_use='', gainuse=0, timerange=[0,0,0,0], flagver=0, bpver=0, image_data=None, ncomp=None, flux=None, smodel=None, aparm=None, dparm=None, snver=1): """Runs FRING to perform Self-Calibration on a multi-source dataset""" print 'trying to FRING', sources_use, 'in dataset', uvdata fring = AIPSTask('fring', version = aips_version_std) fring.indata = uvdata if(type(sources_use) == type('string')): fring.calsour[1] = sources_use else: fring.calsour[1:] = sources_use if(timerange != None): if(timerange[0] == None): fring.timerang = timerange else: fring.timerang[1:] = timerange fring.freqid = 1 if(gainuse >= 1): fring.docalib = 2 fring.gainuse = gainuse fring.flagver=flagver if(bpver > 0): fring.doband = 3 fring.bpver = bpver fring.refant = refant if(search_list != None): if(search_list[0] == None): fring.search=search_list else: fring.search[1:]=search_list fring.solint= solint / 60.0 # convert from s to min fring.solsub = solsub if(aparm != None): if(aparm[0] == None): fring.aparm = aparm else: fring.aparm[1:] = aparm if(dparm != None): if(dparm[0] == None): fring.dparm = dparm else: fring.dparm[1:] = dparm if (type(image_data) != type(None)): # We have an image, hopefully with clean components, to use fring.in2data = image_data if (ncomp != None): if(type(ncomp) == type(1)): fring.ncomp[1] = ncomp else: fring.ncomp[1:] = ncomp fring.nmaps = len(ncomp) if(flux != None): fring.flux = flux else: # no image, so smodel needs to be checked if(smodel != None): if((type(smodel) == type(1)) or (type(smodel) == type(1.0))): fring.smodel[1] = smodel else: fring.smodel[1:] = smodel fring.snver = snver #fring.inp() fring() ################################################################################ # run CLCAL in the basic modes def run_clcal(uvdata, gainver_in=0, gainver_out=0, refant=None, sources='', calsources='', timerange=[0,0,0,0], opcode='CALI', interpol='', samptype='', bparm=None, doblank=0, dobtween=0, smotype='', snver_first=0, snver_last=0): """Run CLCAL for basic things to do""" clcal = AIPSTask('clcal', version = aips_version_std) clcal.indata = uvdata if(type(sources) == type('string')): clcal.sources[1] = sources else: clcal.sources[1:] = sources if(type(calsources) == type('string')): clcal.calsour[1] = calsources else: clcal.calsour[1:] = calsources if(timerange != None): if(timerange[0] == None): clcal.timerang = timerange else: clcal.timerang[1:] = timerange clcal.refant=refant clcal.opcode = opcode clcal.interpol = interpol clcal.samptype = samptype if(bparm != None): if(bparm[0] == None): clcal.bparm = bparm else: clcal.bparm[1:] = bparm clcal.doblank = doblank clcal.dobtween = dobtween clcal.snver = snver_first clcal.invers= snver_last clcal.gainver = gainver_in clcal.gainuse = gainver_out clcal() ################################################################################ # use split to slip out the data def split_data(multi_uvdata, thisclass, thissequence, sources_use='', gainuse = 0, bpver=0, flagver=0): """Split out sources (or all sources) from a multisource file multi_uvdata I the multi-source UV data file thisclass I name of the class to write out split files to sources_use I a string, or a list of strings of the source names to split gainuse I CL table to use (0 for highest) bpver I bandpass correction? flagver I which flag table to use returns a list of all of the split out files """ # if we are doing all sources, then check that if(sources_use == ''): sources_use = multi_uvdata.sources if(type(sources_use) == type('string')): sources_use = [sources_use] # declare a list to hold all of the split files split_file_list = [] # Now, check is we need to zap any of these files for source in sources_use: #print "looking for file '" +source.upper() +"'"+thisclass+"'", multi_uvdata.disk, thissequence s_data = AIPSUVData(source.upper(), thisclass, multi_uvdata.disk, thissequence) if s_data.exists(): warnings.warn('Hey, I found an old split file \"%s\"'% s_data.__str__()) warnings.warn( 'Zapping ....') s_data.zap() else: #print 'No such file found' pass split_file_list.append(s_data) # now do the splitting split = AIPSTask('split', version = aips_version_std) split.indata = multi_uvdata split.sources[1:] = sources_use split.docalib = 2 split.gainuse = gainuse split.flagver=flagver if(bpver > 0): split.doband = 3 split.bpver = bpver split.outclass = thisclass split.outseq = thissequence split.outdisk = multi_uvdata.disk split.aparm[1] = 2 split.aparm[4] = 1 split() # check that all of the files were split properly for s_data in split_file_list: if not s_data.exists(): msg = "Error: cannot find uvdata for split file \"%s\"" \ % (s_data) raise RuntimeError, msg return split_file_list ################################################################################ # bounds check def bounds_check(data_array, blc, trc): # This function assumes that blc is actually the bottom left corner bad_message = 'fully out of bounds' if(blc[0] < 0): if(trc[0] < 0): raise Exception, bad_message elif(trc[0] >= data_array.size(0)): trc[0] = data_array.size(0)-1 blc[0] = 0 elif(blc[0] >= data_array.size(0)): raise Exception, bad_message if(blc[1] < 0): if(trc[1] < 0): raise Exception, bad_message elif(trc[1] >= data_array.size(1)): trc[1] = data_array.size(1)-1 blc[1] = 0 elif(blc[1] >= data_array.size(1)): raise Exception, bad_message ################################################################################ # Some stuff for masking areas def mask_image_rectangle(data_array, blc, trc, inside=1, mask=1,\ C_indices=1): """This function works with a 2D numarray.ma data_array to moodify the mask. It works on a rectangular area from blc=[x,y] to trc=[x,y]. By default, it uses C style (0-based) indices. If you want 1-based indices, set C_indices to 0. Note that this also affects the ordering of the indices. AIPS normally wants you to specify blc as blc=[x,y], but ParselTongue and other Python-based astronomical software has adopted C-based array notation, which would more naturally be blc=[y,x] for the way the data is stored. By default, this function will generate the mask to operate on the inside (edge inclusive) of the rectangle described by blc and trc. If you want to mask the outside (edge EXclusive), then set inside=0. By default, this function will set the mask to be masked whereever the mask already existed, plus where the new rectangle mask indicates. If, on the other hand, you wish to UNMASK regions which may have been previously masked, then set mask=0. For instance, suppose that you have data_array which has size (512,512) and has an existing mask which is masking off the entire array. You want to unmask the central region so that you can work with the pixels at the center. Then call mask_image_rectangle(data_array, [256-10,256-10], [256+10,256+10], inside=1,\ mask=0) """ blc = copy.copy(blc) trc = copy.copy(trc) # If FORTRAN indices, then fix if not C_indices: blc = [blc[1]-1,blc[0]-1] trc = [trc[1]-1,trc[0]-1] # If TRC is actually BLC, swap if(blc[0] > trc[0]): blc[0], trc[0] = trc[0], blc[0] if(blc[1] > trc[1]): blc[1], trc[1] = trc[1], blc[1] # bounds checking try: bounds_check(data_array, blc, trc) except: # if there was an exception raised, then the box is totally # off the image, and there is nothing to do return data_array # now make the mask m = numarray.ma.make_mask_none(data_array.shape) for i in range(blc[0],trc[0]+1): jarray = numarray.arrayrange(0,data_array.size(1)) jarray = numarray.where(numarray.logical_and((jarray >= blc[1]),\ (jarray <= trc[1])), 1, 0) m[i,:] = jarray if (not inside)^(not mask): m = numarray.logical_not(m) if not mask: m = m * numarray.ma.getmaskarray(data_array) else: m = numarray.ma.mask_or(m, numarray.ma.getmask(data_array)) return numarray.ma.array(data_array, mask=m, fill_value=data_array.fill_value()) ################################################################################ def mask_image_circle(data_array, radius, center, inside=1, mask=1,\ C_indices=1): """This function works with a 2D numarray.ma data_array to moodify the mask. It works on a circular area centered at center with radius radius. By default, it uses C style (0-based) indices. If you want 1-based indices, set C_indices to 0. Note that this also affects the ordering of the indices. AIPS normally wants you to specify corner as center=[x,y], but ParselTongue and other Python-based astronomical software has adopted C-based array notation, which would more naturally be center=[y,x] for the way the data is stored. By default, this function will generate the mask to operate on the inside (edge inclusive) of the circle described by radius and center. If you want to mask the outside (edge EXclusive), then set inside=0. By default, this function will set the mask to be masked whereever the mask already existed, plus where the new rectangle mask indicates. If, on the other hand, you wish to UNMASK regions which may have been previously masked, then set mask=0. For instance, suppose that you have data_array which has size (512,512) and has an existing mask which is masking off the entire array. You want to unmask the central region so that you can work with the pixels at the center. Then call mask_image_circle(data_array, 10, [256,256], inside=1,\ mask=0) """ center = copy.copy(center) # If FORTRAN indices, then fix if not C_indices: center = [center[1]-1,center[0]-1] # now make the mask m = numarray.zeros(data_array.shape, type=numarray.Float64) x = numarray.arrayrange(data_array.size(1), type=numarray.Float64) - center[1] y = numarray.arrayrange(data_array.size(0), type=numarray.Float64) - center[0] x = x*x y = y*y radius2 = radius*radius for i in xrange(0,data_array.size(1)): m[:,i] += y for i in xrange(0,data_array.size(0)): m[i,:] += x if (not inside)^(not mask): m = numarray.where(m <= radius2, 0, 1) else: m = numarray.where(m <= radius2, 1, 0) if not mask: m = m * numarray.ma.getmaskarray(data_array) else: m = numarray.ma.mask_or(m, numarray.ma.getmask(data_array)) return numarray.ma.array(data_array, mask=m, fill_value=data_array.fill_value()) ################################################################################ def masked_get_min_max(data_array): """This function takes in a masked array. It then searches for the minimum and maximum pixel values, and records their positions. When finished, it returns the min and max, and their positions. """ data = numarray.ma.filled(data_array, 0) mask = (numarray.ma.getmask(data_array) == 0) (min, max, min_pos, max_pos) = numarray.nd_image.extrema(data,mask) min_pos = numarray.array(min_pos) max_pos = numarray.array(max_pos) return (min, max, min_pos, max_pos) ################################################################################ def get_small_source_stats(data_array, source_radius, noise_inner_radius, \ noise_outer_radius, edge_width=10, warning_level=1): """Get statistis from an image for a \"small\" source. This function will inspect an image to locate a single, strong, mostly circular-ish source. The source will be found from the strongest pixel in the image, excluding edge pixels. This function will the conpute the integrated source brightness within source_radius pixels of the max pixel, and the RMS noise level for pixels at least noise_radius pixels away from max, again, excluding edge pixels. data_array I the incoming image. Note that data_array is expected to be a MaskedArray from NumArray, and may already have some pixels masked. source_radius I the radius from the center of the source to sum pixels noise_inner_radius I inner radius of an annulus to calculate the noise level noise_outer_radius I outer radius of an annulus to calculate the noise level edge_width I width in pixels to avoid at the edge of the image warning_level I how much should the user be warned. 0 None 1 Some It returns information as (max, max_pos, sum, rms), where max is the maximum pixel value max_pos is the position list (array) [y,x] of the maximum pixel position sum is the integrated brightness of the source (may be affected by previous masking of pixels) rms is the RMS value off-source for pixels in radius range noise_inner_radius to noise_outer_radius """ # First, make a version with the edge pixels masked off blc=[edge_width-1,edge_width-1] trc=[data_array.size(0)-edge_width,data_array.size(1)-edge_width] base_array = mask_image_rectangle(data_array, blc, trc, inside=0, mask=1) # Now find the min, max stuff (min, max, min_pos, max_pos) = masked_get_min_max(base_array) # Check for how far this is from center if(warning_level > 0): dimensions = numarray.array(base_array.shape, type=numarray.Float64) blc = dimensions * 0.25 trc = dimensions * 0.75 if(numarray.any(max_pos < blc) or numarray.any(max_pos > trc)): warnings.warn("'%s': Max position out of center of image"\ %get_small_source_stats.func_name) # Now flag everything except the source region source_array = mask_image_circle(base_array, source_radius, max_pos,\ inside=0, mask=1) # Get the integrated value. integ_value = numarray.ma.sum(numarray.ma.ravel(source_array)) del source_array # Now flag the source and get the image rms rms_array = mask_image_circle(base_array, noise_inner_radius, max_pos,\ inside=1, mask=1) rms_array = mask_image_circle(rms_array, noise_outer_radius, max_pos,\ inside=0, mask=1) rms = numarray.nd_image.standard_deviation(numarray.ma.filled(rms_array, 0), (numarray.ma.getmask(rms_array) == 0)) # return (max, max_pos, integ_value, rms) ################################################################################ def get_small_source_stats_main(aips_image, \ source_radius, noise_inner_radius, \ noise_outer_radius, edge_width=10, \ print_level=1): """User version for statistis from an image for a \"small\" source. This function takes an AIPSImage, and returns physical (Jy/beam and Jy) values. Note that the returned position is a C position!!! This function will inspect an image to locate a single, strong, mostly circular-ish source. The source will be found from the strongest pixel in the image, excluding edge pixels. This function will the conpute the integrated source brightness within source_radius pixels of the max pixel, and the RMS noise level for pixels at least noise_radius pixels away from max, again, excluding edge pixels. aips_image I the incoming image as an AIPSImage source_radius I the radius from the center of the source to sum pixels noise_inner_radius I inner radius of an annulus to calculate the noise level noise_outer_radius I outer radius of an annulus to calculate the noise level edge_width I width in pixels to avoid at the edge of the image print_level I how much should the user be informed. 0 None 1 Some It returns information as (max, max_pos, sum, rms), where max is the maximum pixel value, in specific intensity (Jy/beam) max_pos is the position list (array) [y,x] of the maximum pixel position sum is the integrated brightness of the source (may be affected by previous masking of pixels). Note that this function returns the sum in flux density (Jy) units, rms is the RMS value off-source for pixels in radius range noise_inner_radius to noise_outer_radius (Jy/beam) """ # Get a Wizardry image to work on the pixels w = Wizardry.AIPSData.AIPSImage(aips_image.name, aips_image.klass, \ aips_image.disk, aips_image.seq) # make a masked array x = numarray.ma.array(w.pixels.copy(), \ mask=numarray.ma.make_mask_none((aips_image.header.naxis[1],aips_image.header.naxis[0]))) (image_max, image_max_pos, image_sum, image_rms) = \ get_small_source_stats(x, source_radius, noise_inner_radius, noise_outer_radius, print_level) # Now correct the sum for the beam size factor = aips_image.header.bmaj * aips_image.header.bmin / \ math.fabs(aips_image.header.cdelt[1] * aips_image.header.cdelt[0]) image_sum = image_sum / factor if(print_level > 0): print 'Got information from image of:' print 'Max:', image_max print 'Pos:', image_max_pos print 'Sum:', image_sum print 'RMS:', image_rms return (image_max, image_max_pos, image_sum, image_rms) ################################################################################ def get_jmfit_statistics(aips_image, max_position, middle_name ): """Run JMFIT and get the statistics in a dictionary aips_image, # an AIPSImage object max_position, # C notation, as an array (numarray) middle_name # middle part of a file(name) to have JMFIT write to """ outfile = dirname + '/' + experiment + '.' + middle_name + '.jmfit' # if the file exists, delete it if os.path.isfile(outfile): os.remove(outfile) # form a FORTRAN position for AIPS image_max_pos_FORTRAN = max_position[-1:None:-1] + 1 # do the actual JMFIT call jmfit_1(aips_image, outfile, center=image_max_pos_FORTRAN) # JMFIT finds lots of stuff, but doesn't return it. It is only available # in the output file. So get it as a dictionary jm_dict = {} try: # open the file and read everything in fp = open(outfile, "r") jm_data = fp.read() fp.close() # Ok, where is the start of the main results section? main_index = jm_data.index("Solution from JMFIT") # Ok, get the peak stuff this_index = jm_data.index("Peak intensity", main_index) words = jm_data[this_index + len('Peak intensity ='):].split(None,6) jm_dict['Peak_Int'] = float(words[0]) jm_dict['Peak_Int_Err'] = float(words[2]) # Ok, get the Integrated stuff this_index = jm_data.index("Integral intensity", main_index) words = jm_data[this_index + len('Integral intensity='):].split(None,5) jm_dict['Int_Int'] = float(words[0]) jm_dict['Int_Int_Err'] = float(words[2]) # Ok, get the X position this_index = jm_data.index("X-position", main_index) words = jm_data[this_index:].split(None,5) jm_dict['X_pos'] = float(words[2]) jm_dict['X_pos_Err'] = float(words[4]) # Ok, get the Y position this_index = jm_data.index("Y-position", main_index) words = jm_data[this_index:].split(None,5) jm_dict['Y_pos'] = float(words[2]) jm_dict['Y_pos_Err'] = float(words[4]) # Ok, get the RA position this_index = jm_data.index("RA", main_index) words = jm_data[this_index:].split(None,6) jm_dict['RA'] = float(words[3]) + (float(words[2]) \ + (float(words[1])) * 60.0)*60.0 jm_dict['RA_Err'] = float(words[5]) jm_dict['RA_Str'] = words[1] + ':' + words[2] + ':' + words[3] # Ok, get the Dec position. Be carefule of DECember this_index = jm_data.index(" DEC", main_index) words = jm_data[this_index:].split(None,6) jm_dict['DEC'] = float(words[3]) + (float(words[2]) \ + (math.fabs(float(words[1]))) * 60.0)*60.0 if(words[1][0] == '-'): jm_dict['DEC'] = -jm_dict['DEC'] jm_dict['DEC_Err'] = float(words[5]) jm_dict['DEC_Str'] = words[1] + ':' + words[2] + ':' + words[3] # Ok, get the deconvolved Major Axis main_index = jm_data.index("Deconvolution of component in asec",main_index) this_index = jm_data.index("Major ax", main_index) words = jm_data[this_index:].split(None,5) jm_dict['Comp_Maj_Ax'] = [float(words[2]), float(words[3]), float(words[4])] # Ok, get the deconvolved Minor Axis this_index = jm_data.index("Minor ax", main_index) words = jm_data[this_index:].split(None,5) jm_dict['Comp_Min_Ax'] = [float(words[2]), float(words[3]), float(words[4])] # Ok, get the deconvolved Position Angle this_index = jm_data.index("Pos ang", main_index) words = jm_data[this_index:].split(None,5) jm_dict['Comp_Pos_Ang'] = [float(words[2]), float(words[3]), float(words[4])] except: msg = "Error while trying to read JMFIT results" raise return jm_dict ################################################################################ # Get the data from the archive, unless it is already here def getoutarchive(): inurl = 'http://archive.jive.nl/exp/' for idi in range(1,nfits+1): file = experiment.lower() + '_1'\ +'_1.IDI'+str(idi) outfile = dirname + '/' + file command = 'wget -nv -O '+outfile+' '\ +'--http-user='+experiment+' '\ +'--http-passwd='+passwd+' '\ +inurl+experiment.upper()+'_'+date+'/fits/'+file #print command runwget = 0 if not os.path.exists(outfile): os.system(command) runwget = 1 elif os.path.getsize(outfile) < 1: os.system(command) runwget = 1 else: print outfile, 'exists!' if (runwget): output = 'Fetched data from the Archive using wget:\n' output += '>' + command else: output = 'No need to get data from archive this time\n' print output for exten in ['uvflg', 'chflag', 'antab']: file = experiment.lower() + '.' + exten outfile = dirname + '/' + file command = 'wget -nv -O '+outfile+' '\ +'--http-user='+experiment+' '\ +'--http-passwd='+passwd+' '\ +inurl+experiment.upper()+'_'+date+'/pipe/'+file print command runwget = 0 if not os.path.exists(outfile): os.system(command) runwget = 1 elif os.path.getsize(outfile) < 1: os.system(command) runwget = 1 else: print outfile, 'exists!' if (runwget): output = 'Fetched data from the Archive using wget:\n' output += '>' + command else: output = 'No need to get data from archive this time\n' print output return ################################################################################ def check_data_file(extension): file = dirname + '/' + experiment + '_' + '.' + extension if not os.path.isfile(file): file = dirname + '/' + experiment + '.' + extension assert(os.path.isfile(file)) return file ################################################################################ def calibration_step_1(): getoutarchive() ################################################################################ # Now load the data into AIPS def calibration_step_2(): print "Creating a new AIPS data set as", uvdata if uvdata.exists(): warnings.warn('zapping old uvdata data') uvdata.zap() fitld = AIPSTask('fitld', version='31DEC04') fitld.infile = dirname + '/' + fits_file fitld.outdata = uvdata fitld.ncount = nfits fitld.douvcomp = 0 fitld.doconcat = 1 fitld.clint = cl_table_interval / 60.0 # convert from s to min fitld.wtthresh = 0.05 fitld.digicor=-1 fitld.inp() fitld() ################################################################################ # Now sort def calibration_step_3(msortdata): """sort the uvdata""" print 'msortdata.exists =', msortdata.exists() if msortdata.exists(): warnings.warn('zapping old msort data \"%s\"'%msortdata) msortdata.zap() #s = raw_input('return') msort = AIPSTask('msort', version=aips_version_try) msort.indata = uvdata msort.outdata = msortdata msort() print "Deleting the original dataset", uvdata uvdata.zap() ################################################################################ # get some useful information out of the UVData def calibration_step_3_after(): """find information about the observations, and stuff that into globals""" global refant, antenna_dict, antenna_array,\ calibrator_instrument_search_num, all_sources, fring_calibrators # get the reference antenna number refant = get_antenna_number(refantname) antenna_dict = get_antenna_dictionary() antenna_array = [None] antenna_array[1:] = uvdata.antennas calibrator_instrument_search_num = range(0,len(calibrator_instrument_search)) for antenna in range(0,len(calibrator_instrument_search)): calibrator_instrument_search_num[antenna] = antenna_dict[calibrator_instrument_search[antenna]] # print uvdata.antennas # print antenna_dict # print antenna_array all_sources = uvdata.sources assert(len(all_sources) > 0) print 'Found sources ', all_sources if(fring_calibrators == None): fring_calibrators = all_sources assert(len(fring_calibrators) > 0) for calibrator in fring_calibrators: assert(calibrator in all_sources) stokes = uvdata.stokes assert(len(stokes) > 0) print 'This dataset has Stokes', stokes polarizations = uvdata.polarizations assert(len(polarizations) > 0) print 'This dataset has polarizations', polarizations ################################################################################ # A priori data flagging and indexing def calibration_step_4(): # Delete any old flag tables. uvdata.zap_table('AIPS FG', -1) uvflg = AIPSTask('uvflg', version = aips_version_std) uvflg.msgkill = -6 uvflg.indata = uvdata uvflg.flagver = 1 if uvflg_file: uvflg.infile = uvflg_file uvflg() if chflg_file: uvflg.infile = chflg_file uvflg() ################################################################################ # index the data def calibration_step_5(): indxr = AIPSTask('indxr', version=aips_version_try) indxr.indata = uvdata indxr.cparm[2] = 22.0 indxr.cparm[3] = cl_table_interval / 60.0 # convert from s to min indxr() ################################################################################ # run listr def calibration_step_6(): listr = AIPSTask('listr', version = aips_version_std) listr.indata = uvdata listr.optype = 'SCAN' outfile = dirname + '/' + experiment + '.scan' if(os.path.isfile(outfile)): warnings.warn("File '%s' already exists. Deleting."%outfile) os.remove(outfile) listr.outprint = outfile listr.docrt = 0 listr() ################################################################################ # run dtsum def calibration_step_7(): dtsum = AIPSTask('dtsum', version = aips_version_std) dtsum.indata = uvdata dtsum.aparm[1] = 1.0 outfile = dirname + '/' + experiment + '.dtsum' if(os.path.isfile(outfile)): warnings.warn("File '%s' already exists. Deleting."%outfile) os.remove(outfile) dtsum.outprint = outfile dtsum.docrt = 0 dtsum() if(pause_after_plots): print "Have a look at the possm plots" possm_1(all_sources, 'HALF', 1, -1) possm_2(all_sources, 'HALF', 1, -1, refant) ################################################################################ # now begin the initial amplitude calibration def calibration_step_8(): uvdata.zap_table('AIPS SN', -1) while uvdata.table_highver('AIPS CL') > 1: uvdata.zap_table('AIPS CL', 0) while uvdata.table_highver('AIPS TY') > 1: uvdata.zap_table('AIPS TY', 0) while uvdata.table_highver('AIPS GC') > 1: uvdata.zap_table('AIPS GC', 0) if antab_file: antab = AIPSTask('antab', version=aips_version_try) antab.infile = antab_file antab.indata = uvdata antab.tyver = 1 antab.gcver = 1 antab.offset = 0.3 antab() elif not uvdata.table_highver('AIPS TY') > 0: raise RuntimeError, "AIPS TY table missing" elif not uvdata.table_highver('AIPS GC') > 0: raise RuntimeError, "AIPS GC table missing" apcal = AIPSTask('apcal', version=aips_version_try) apcal.indata = uvdata apcal.freqid = 1 apcal.tyver = 1 apcal.gcver = 1 apcal.snver = 1 apcal.prtlev = 1 apcal() clcal = AIPSTask('clcal', version = aips_version_std) clcal.indata = uvdata clcal.opcode = 'CALI' # interpol should be 'SELF' when the density of calibration # points is high. Unfortunately, this is a fast phase referencing # observation clcal.interpol = 'BOX' clcal.doblank = 1 clcal.samptype = 'BOX' clcal.bparm[1] = 1E-5 clcal.refant = refant clcal.snver = 1 clcal.gainver = Gainuse_Global clcal.gainuse = Gainuse_Global + 1 # first do CM, because CM has no real calibration clcal.interpol = '2PT' clcal.antennas[1] = antenna_dict['CM'] clcal() # Now do everything else clcal.interpol = 'SELF' clcal.antennas[1] = -antenna_dict['CM'] clcal() print 'You should check the results here' snplot('TY', 1, 'TSYS', opcode='ALSI', snuvdata=uvdata) snplot('CL', Gainuse_Global+1, 'AMP', opcode='ALSI', snuvdata=uvdata) ################################################################################ # now do the paralactic angle change def calibration_step_9(): while uvdata.table_highver('AIPS CL') > Gainuse_Global: uvdata.zap_table('AIPS CL', 0) assert(uvdata.table_highver('AIPS CL') == Gainuse_Global) clcor = AIPSTask('clcor', version = aips_version_std) clcor.indata = uvdata clcor.gainver = Gainuse_Global clcor.gainuse = Gainuse_Global + 1 clcor.opcode = 'PANG' clcor.clcorprm[1] = 1 clcor() ################################################################################ # now do some manual flagging def calibration_step_10(): print 'Now do some manual flagging with EDITR' fg_highver = uvdata.table_highver('AIPS FG') if(fg_highver != 0): # copy FG highver to FG highver+1 tacop = AIPSTask('tacop', version=aips_version_std) tacop.indata = uvdata tacop.outdata = uvdata tacop.inext = 'FG' tacop.inver = fg_highver tacop.outver = fg_highver+1 tacop.ncount = 1 tacop() fg_highver = fg_highver+1 for source in all_sources: # print out the telescopes #for antenna in uvdata.antennas: # print antenna_dict[antenna], antenna print 'Working on source', source # create and editr object editr1(sources_use=source,gainuse=Gainuse_Global,flagver=fg_highver, uvdata_use=uvdata) # Now write out the flag table to keep it safe outfile = dirname + '/' + 'flag%d.dat'%fg_highver if(os.path.isfile(outfile)): warnings.warn("File '%s' already exists. Deleting."%outfile) os.remove(outfile) tbout = AIPSTask('tbout') tbout.indata = uvdata tbout.inext = 'FG' tbout.inver = fg_highver tbout.outfile = outfile tbout.docrt = 256 tbout() # check that the user is happy if(pause_after_plot): pause_prog("At the end of EDITR, are you happy?") # now run possm again for source in all_sources: # print out the telescopes #for antenna in uvdata.antennas: # print antenna_dict[antenna], antenna print 'Working on source', source # create and editr object possm_2(source, 'LL', 3, -1, refant,fg_highver) # uvdata.zap_table('AIPS FG', 3) # tacop = AIPSTask('tacop') # tacop.indata = uvdata # tacop.outdata = uvdata # tacop.inext = 'FG' # tacop.inver = 1 # tacop.outver = 3 # tacop.ncount = 1 # tacop() # del tacop # for source in all_sources: # # print out the telescopes # #for antenna in uvdata.antennas: # # print antenna_dict[antenna], antenna # print 'Working on source', source # # create and editr object # editr1(sources_use=source,gainuse=3,flagver=6,uvdata_use=uvdata) # tbout = AIPSTask('tbout') # tbout.indata = uvdata # tbout.inext = 'FG' # tbout.inver = 6 # tbout.outfile = dirname + '/' + 'flag.dat' # tbout.docrt = 256 # tbout() # del tbout ################################################################################ # load a table in def load_ascii_table(aips_data, infile): if not(os.path.isfile(infile)): # No table to load!!! warnings.warn("I can't load file '%s', no such file."%infile) raise RuntimeError, "ASCII table file missing" tbin = AIPSTask('tbin') tbin.outdata = aips_data tbin.infile = infile tbin() del tbin # # need to remove the bad ifs from CM data # uvflg = AIPSTask('uvflg') # uvflg.indata = uvdata # uvflg.antenns[1] = antenna_dict['CM'] # uvflg.flagver = 6 # uvflg.opcode = 'FLAG' # uvflg.reason = 'bad response' # uvflg.dohist = 1 # uvflg.bif = 2 # uvflg.eif = 2 # uvflg() # uvflg.bif = 4 # uvflg.eif = 4 # uvflg() #final_flagver = 7 # for source in all_sources: # # print out the telescopes # #for antenna in uvdata.antennas: # # print antenna_dict[antenna], antenna # print 'Working on source', source # # create and editr object # possm_2(source, 'LL', 3, -1, refant,final_flagver) ################################################################################ # FRING on the bright sources def calibration_step_11(): print """Fringe fit the data. First make sure that old SN and CL tables are not lying around""" uvdata.zap_table('AIPS SN', -1) while uvdata.table_highver('AIPS CL') > Gainuse_Global: uvdata.zap_table('AIPS CL', 0) assert(uvdata.table_highver('AIPS CL') == Gainuse_Global) fring = AIPSTask('fring', version = aips_version_std) fring.indata = uvdata fring.outdata = uvdata if(type(calibrator_instrument) == type('string')): fring.calsour[1] = calibrator_instrument else: fring.calsour[1:] = calibrator_instrument fring.freqid = 1 fring.docalib = 2 fring.gainuse = Gainuse_Global fring.flagver = final_flagver fring.timerang[1:] = calibrator_instrument_timer fring.subarray = 1 fring.refant = calibrator_instrument_search_num[0] if(len(calibrator_instrument_search_num) > 1): fring.search[1:] = calibrator_instrument_search_num[1:] fring.solint = calibrator_instrument_solint / 60.0 # convert from s to min fring.aparm[1] = 4 fring.aparm[6] = 1 fring.aparm[7] = 11 fring.aparm[9] = 1 fring.dparm[1] = 2 fring.dparm[4] = t_int # set dparm[8] to 1 to zero rates fring.dparm[8] = 0 fring.snver = 1 fring() clcal = AIPSTask('clcal', version = aips_version_std) clcal.indata = uvdata if(type(calibrator_instrument) == type('string')): clcal.calsour[1] = calibrator_instrument else: clcal.calsour[1:] = calibrator_instrument clcal.subarray = 1 clcal.opcode = 'CALI' clcal.interpol = 'AMBG' clcal.refant = calibrator_instrument_search_num[0] clcal.snver = 1 clcal.gainver = Gainuse_Global clcal.gainuse = Gainuse_Global + 1 clcal() snplot('SN',1, 'PHAS', opcode='ALSI',timerange = calibrator_instrument_timer, snuvdata=uvdata, nplots=len(uvdata.antennas)) snplot('SN',1, 'DELA', opcode='ALSI',timerange = calibrator_instrument_timer, snuvdata=uvdata, nplots=len(uvdata.antennas)) snplot('CL',Gainuse_Global+1, 'PHAS', opcode='ALSI', timerange = calibrator_instrument_timer, snuvdata=uvdata, nplots=len(uvdata.antennas)) snplot('CL',Gainuse_Global+1, 'PHAS', opcode='ALSI', snuvdata=uvdata, nplots=len(uvdata.antennas)) if(pause_after_plots): possm_2(calibrator_instrument, 'LL', Gainuse_Global, -1,refant,final_flagver) for source in all_sources: # print out the telescopes #for antenna in uvdata.antennas: # print antenna_dict[antenna], antenna print 'Working on source', source # create and editr object possm_2(source, 'LL', Gainuse_Global+1, -1, refant,final_flagver) ################################################################################ # bandpass calibration def calibration_step_12(): print 'Bandpass calibration' uvdata.zap_table('AIPS SN', -1) uvdata.zap_table('AIPS BP', -1) while uvdata.table_highver('AIPS CL') > Gainuse_Global: uvdata.zap_table('AIPS CL', 0) assert(uvdata.table_highver('AIPS CL') == Gainuse_Global) bpass = AIPSTask('bpass', version = aips_version_std) bpass.indata = uvdata if(type(calibrator_instrument) == type('string')): bpass.calsour[1] = calibrator_instrument else: bpass.calsour[1:] = calibrator_instrument bpass.freqid = 1 bpass.docalib = 2 bpass.gainuse = Gainuse_Global bpass.flagver = final_flagver bpass.timerang[1:] = calibrator_instrument_timer bpass.refant = calibrator_instrument_search_num[0] bpass.solint = 0 #bpass.soltype = 'L1' bpass.outver = Bandpass_Global + 1 #bpass.smooth[1:] = [1,0] #bpass.bpassprm[2] = 13 bpass.bpassprm[5] = 1 bpass.bpassprm[9] = 1 bpass.bpassprm[10] = 3 #bpass.ichansel[1] = [None, 3,14,1,1] #bpass.ichansel[2] = [None, 3,14,1,2] #bpass.ichansel[3] = [None, 3,14,1,3] #bpass.ichansel[4] = [None, 3,14,1,4] #bpass.ichansel[5] = [None, 3,-1,1,5] bpass() bplot = AIPSTask('bplot', version = aips_version_std) bplot.indata = uvdata bplot.bpver = Bandpass_Global + 1 bplot.codetype = 'AMP' bplot.do3col = 1 bplot.dotv = 1 if(bplot.dotv > 0): tv_local.clear() bplot() if(pause_after_plots): pause_prog('Finished BPLOT AMP. Continue?') bplot.codetype = 'PHAS' bplot() if(pause_after_plots): pause_prog('Finished BPLOT PHAS. Continue?') possm_2(calibrator_instrument, 'LL', Gainuse_Global, -1, refant, final_flagver,bpver=Bandpass_Global) # tbout = AIPSTask('tbout') # tbout.indata = uvdata # tbout.inext = 'BP' # tbout.inver = 1 # tbout.outfile = dirname + '/' + 'bpass.dat' # tbout.docrt = 256 # tbout() # del tbout ################################################################################ # Ok, now fring the remaining sources # I want to manually make some good images to get # a feel for how the sources look, and to get some nice clean components # for going back to fring if the sources are complicated def calibration_step_13(): print 'Phase Calibrator Fring' uvdata.zap_table('AIPS SN', -1) assert(uvdata.table_highver('AIPS BP') == 1) while uvdata.table_highver('AIPS CL') > Gainuse_Global: uvdata.zap_table('AIPS CL', 0) assert(uvdata.table_highver('AIPS CL') == Gainuse_Global) for source in fring_calibrators: print 'Working on source', source uvdata.zap_table('AIPS SN', -1) fring = AIPSTask('fring', version = aips_version_std) fring.indata = uvdata fring.outdata = uvdata fring.calsour[1] = source fring.freqid = 1 fring.docalib = 2 fring.gainuse = Gainuse_Global fring.doband = 3 fring.bpver = Bandpass_Global fring.flagver = final_flagver fring.subarray = 0 fring.refant = calibrator_instrument_search_num[0] if(len(calibrator_instrument_search_num) > 1): fring.search[1:] = calibrator_instrument_search_num[1:] fring.solint = solint / 60.0 # convert from s to min fring.aparm[1] = 3 fring.aparm[6] = 0 fring.aparm[7] = 10 fring.aparm[9] = 1 fring.dparm[1] = 2 fring.dparm[4] = t_int fring.dparm[8] = 0 fring.snver = 1 fring() clcal = AIPSTask('clcal', version = aips_version_std) clcal.indata = uvdata clcal.calsour[1] = source clcal.sources[1] = source clcal.subarray = 1 clcal.opcode = 'CALI' clcal.interpol = 'SIMP' clcal.refant = calibrator_instrument_search_num[0] clcal.snver = 1 clcal.gainver = Gainuse_Global clcal.gainuse = Gainuse_Global + 1 clcal() print 'Here come the phase results' snplot('SN', 1, 'PHAS', opcode='ALSI', snuvdata=uvdata, nplots=len(uvdata.antennas)) snplot('SN', 1, 'DELA', opcode='ALSI', snuvdata=uvdata, nplots=len(uvdata.antennas)) snplot('SN', 1, 'RATE', opcode='ALSI', snuvdata=uvdata, nplots=len(uvdata.antennas)) snplot('CL', Gainuse_Global+1, 'PHAS', opcode='ALSI',sources_use=source, snuvdata=uvdata,nplots=len(uvdata.antennas)) # for source in all_sources: # # print out the telescopes # #for antenna in uvdata.antennas: # # print antenna_dict[antenna], antenna # print 'Working on source', source # # create and editr object # possm_2(source, 'LL', Gainuse_Global, -1, refant,final_flagver,bpver=Bandpass_Global) ################################################################################ def set_cellsize_info(uv_data): """get the proper cellsize fo the observations somehow""" global cellsize obsfreq = uv_data.header['crval'][2] obslambda = 2.99792458E8/obsfreq baseline = 7000E3 # meters => 7000 km cellsize = obslambda/baseline # radians cellsize = cellsize*206264.806 # arcsec cellsize = cellsize/4. # I like 4 cells per resolution element ################################################################################ def check_table_version_valid(aips_data, table, version): """Check that version is a valid table, as 'AIPS CL', 'AIPS BP', etc.""" # Check for ANY version of the table. 0 implies 0 or more if(version == 0): return 1 # make sure there is at least ONE table #if(aips_data.table_highver(table) == 0): # raise AttributeError, "No table"+table for t in aips_data.tables: if((t[0] == version) and (t[1] == table)): break else: return 0 return 1 ################################################################################ def check_table_version_exists(aips_data, table, version): """Check whether a table version exists, as 'AIPS CL', 'AIPS BP', etc. if version == 0, then check for ANY of that table, """ # If version is zero, then accept any of that table if(version == 0): for t in aips_data.tables: if(t[1] == table): break else: return 0 return 1 # Otherwise, check for the version number too for t in aips_data.tables: if((t[0] == version) and (t[1] == table)): break else: return 0 return 1 ################################################################################ # Find the highest numbered antenna def find_antenna_highnumber(aips_data): """Find the highest numbered antenna in the dataset """ highnum = 0 for an in xrange(1,aips_data.table_highver('AIPS AN')+1): # check that this exists if not(check_table_version_exists(aips_data, 'AIPS AN', an)): continue # Ok, grab the AN table an_table = aips_data.table('AIPS AN', an) for row in an_table: if(highnum < row.nosta): highnum = row.nosta return highnum ################################################################################ # Count visibilities def count_visibilities(aips_data_orig): """Count the number of visibilities per antenna in a dataset aips_data_orig I a non-Wizardry object to use """ aips_data = Wizardry.AIPSData.AIPSUVData(aips_data_orig.name, aips_data_orig.klass, aips_data_orig.disk, aips_data_orig.seq) # What is the highest valued antenna number? high_ant = find_antenna_highnumber(aips_data_orig) print high_ant antennas = numarray.zeros((high_ant+1)) for row in aips_data: baseline = row.baseline a1 = baseline[0] a2 = baseline[1] if(a1 != a2): antennas[a1] = antennas[a1] +1 antennas[a2] = antennas[a2] +1 return antennas ################################################################################ # Add the standard antenna number of visibilities to a dictionary def add_visibility_count_to_dict(d, vis): """ Add the standard antenna number of visibilities to a dictionary d B the dictionary to add stuff to vis I the number of visibilities, as a list or array This function uses antenna_array and antenna_dict to get antenna info """ max = len(vis) if(max > len(antenna_array)): max = len(antenna_array) for i in xrange(max): if(antenna_array[i] == None): continue d['Ant_Vis_Count_'+antenna_array[i]] = vis[i] return ################################################################################ def make_selfcal_calib_data_in(main_data_file, source_name, calib_class, sequence_in): """Create an AIPSData variable with the appropriate name, class, seq""" if((sequence_in == None) or (sequence_in <= 0)): # use the main data file calib_in = main_data_file else: if((source_name == None) or (source_name == '')): source_name = main_data_file.name if((calib_class == None) or (calib_class == '')): calib_class = 'CALIB' calib_in = AIPSUVData(source_name, calib_class, main_data_file.disk, sequence_in) assert(calib_in.exists()) # Clean up. Remove all SN files. if(calib_in.table_highver('AIPS SN') > 0): calib_in.zap_table('AIPS SN', -1) #assert(check_table_version_valid(calib_in, 'AIPS CL', gainuse_start)) #assert(check_table_version_valid(calib_in, 'AIPS BP', bpver)) return calib_in ################################################################################ def make_selfcal_calib_data_in_nocheck(main_data_file, source_name, calib_class, sequence_in): """Create an AIPSData variable with the appropriate name, class, seq""" if((sequence_in == None) or (sequence_in <= 0)): # use the main data file calib_in = main_data_file else: if((source_name == None) or (source_name == '')): source_name = main_data_file.name if((calib_class == None) or (calib_class == '')): calib_class = 'CALIB' calib_in = AIPSUVData(source_name, calib_class, main_data_file.disk, sequence_in) return calib_in ################################################################################ def make_selfcal_calib_data_out(calib_data_in, source_name, calib_class, sequence_out): """Create an AIPSData variable with the appropriate name, class, seq""" if((source_name == None) or (source_name == '')): source_name = calib_data_in.name if((calib_class == None) or (calib_class == '')): calib_class = 'CALIB' if((sequence_out == None) or (sequence_out <= 0)): sequence_out = 1 calib_out = AIPSUVData(source_name, calib_class, calib_data_in.disk, sequence_out) if(calib_out.exists()): if __debug__: #warnings.warn("'%s': CALIB file \"%s\" already exists. Hit Enter to delete."\ # %(make_selfcal_calib_data_out.func_name,calib_out)) #pause_prog('Continue?') pass calib_out.zap() return calib_out ################################################################################ def make_selfcal_imagr_outfile(calib_data_in, source_name, sequence_out): """Create an AIPSImage variable with the appropriate name, class, seq""" if((source_name == None) or (source_name == '')): source_name = calib_data_in.name if((sequence_out == None) or (sequence_out <= 0)): sequence_out = 1 # Make the Beam image. image_class = 'IBM001' image_data = AIPSImage(source_name, image_class, calib_data_in.disk, sequence_out) # If this exists, delete it, as we may be using a different imsize if(image_data.exists()): if __debug__: warnings.warn("'%s': Beam file \"%s\" already exists. Hit Enter to delete."\ %(make_selfcal_imagr_outfile.func_name,image_data)) if(pause_before_overwrite): pause_prog('Continue?') image_data.zap() # Now make the image image image_class = 'ICL001' image_data = AIPSImage(source_name, image_class, calib_data_in.disk, sequence_out) # If this exists, delete it, as we may be using a different imsize if(image_data.exists()): if __debug__: warnings.warn("'%s': Image file \"%s\" already exists. Hit Enter to delete."\ %(make_selfcal_imagr_outfile.func_name,image_data)) if(pause_before_overwrite): pause_prog('Continue?') image_data.zap() return image_data ################################################################################ def check_selfcal_flagver(aips_data, flagver): if(check_table_version_valid(aips_data, 'AIPS FG', flagver)): return flagver else: return 0 ################################################################################ def is_solmode_doing_amplitude_single(solmode): """is solmode doing aplitude calibration for single source file?""" flag = 0 if((solmode == None) or (solmode == '')): pass elif((solmode[0] == 'A') or (solmode[0] == 'G')): flag=1 return flag ################################################################################ def make_default_selfcal_calib_parameters_dict(imsize, solint, C_center, source_size, combine_data, weak_source): """Setup default CALIB and IMAGR parameters for selfcal runs solint I the solution interval, in seconds C_center I the center location for a clean box for IMAGR, as a C array source_size I the radial size of the source, in pixels combine_data I Should the polarizations and IFs be combined in the fit? weak_source I Should the RMS noise cutoff be lowered? """ cal_par = {} cal_par['run_imagr'] = 1 # run IMAGR to make an image to selfcal cal_par['imagr_flux_fraction'] = 0.001# Don't trust things below 0.1% cal_par['imagr_noise_cutoff'] = 1.0 # Don't trust things below 1-sigma cal_par['niter'] = 300 # normally do up to 200 clean iters cal_par['dotv'] = 0 # For now, turn off the tv cal_par['clbox'] = [None,-1,source_size,C_center[1]+1,C_center[0]+1] cal_par['sequence_in'] = 0 # use the original data cal_par['sequence_out'] = 1 # and make sequence 1 cal_par['solint'] = solint # user specified cal_par['solsub'] = 2 # step by half a solint cal_par['ncomp'] = 100 # take first 100 components + or - cal_par['flux'] = 0.0 # take all flux valued components cal_par['soltype'] = '' # use the default least squares cal_par['solmode'] = 'P!A' # phase only self-cal cal_par['aparm'] = [None, 3,0,0,0,0,2,4,0,1] # default aparm # This sets the min number antennas 3 # some print info level 2 # minimum RMS 4 # keep data from failed solutions if(combine_data): cal_par['aparm'][3] = 1 cal_par['aparm'][5] = 1 if(weak_source): cal_par['aparm'][7] = 3 cal_par['show_SN_table'] = 0 # do not run SNPLT for now return cal_par ################################################################################ def make_default_selfcal_calib_parameters_list(num_selfcals, imsize, solint, solint_amp, C_center, source_size, combine_data, weak_source, do_ampl_cal): """Make an array of dictionaries to control CALIB and IMAGR for selfcal num_selfcals I how many self-cals do you want? solint I the solution interval, in seconds for phase-only calibration solint_amp I the solution intercal in seconds for amplitude selfcal C_center I the center location for a clean box for IMAGR, as a C array source_size I the radial size of the source, in pixels combine_data I Should the polarizations and IFs be combined in the fit? weak_source I Should the RMS noise cutoff be lowered? do_ampl_cal I should an amplitude calibration be done? """ cal_par = range(0,num_selfcals) # make a list # Now, for each element in the list, assign the default dict to it for i in xrange(0,num_selfcals): cal_par[i] = \ make_default_selfcal_calib_parameters_dict(imsize,solint,C_center, source_size, combine_data, weak_source) # Now go back and modify niter, ncomp, etc. for the first few iterations if(num_selfcals > 2): cal_par[0]['niter'] = 50 cal_par[0]['ncomp'] = 5 cal_par[0]['clbox'] = [None,-1,6,C_center[1]+1,C_center[0]+1] cal_par[0]['solint'] = cal_par[0]['solint']*2 cal_par[0]['aparm'][7] = 3 if(num_selfcals > 3): cal_par[1]['niter'] = 75 cal_par[1]['ncomp'] = 10 cal_par[1]['clbox'] = [None,-1,10,C_center[1]+1,C_center[0]+1] cal_par[1]['solint'] = cal_par[0]['solint']*1.5 if(num_selfcals > 4): cal_par[2]['niter'] = 100 cal_par[2]['ncomp'] = 15 cal_par[2]['clbox'] = [None,-1,15,C_center[1]+1,C_center[0]+1] if(num_selfcals > 5): cal_par[3]['niter'] = 150 cal_par[3]['ncomp'] = 20 # show the first and last SN tables cal_par[0]['show_SN_table'] = 1 cal_par[num_selfcals-1]['show_SN_table'] = 1 # check amplitude calibration stuff if(do_ampl_cal): cal_par[num_selfcals-1]['solint'] = solint_amp cal_par[num_selfcals-1]['solmode'] = 'A&P' if(num_selfcals > 1): cal_par[num_selfcals-1]['sequence_in'] = 1 # the phase calibrated set cal_par[num_selfcals-1]['sequence_out'] = 2# and make sequence 2 cal_par[num_selfcals-1]['run_imagr'] = 0 cal_par[num_selfcals-2]['show_SN_table'] = 1 if((do_ampl_cal) and (num_selfcals > 8)): cal_par[num_selfcals-3]['solint'] = solint_amp cal_par[num_selfcals-3]['solmode'] = 'A&P' cal_par[num_selfcals-3]['sequence_in'] = 1 # use the phase calibrated set cal_par[num_selfcals-3]['sequence_out'] = 2# and make sequence 2 cal_par[num_selfcals-3]['run_imagr'] = 0 # that's all from here return cal_par ################################################################################ # Run a series of IMAGR clean and CALIB self-cal's on a Single Source Dataset def self_cal_source(aips_data, source_name='', flagver = 0, bpver = -1, refant = 1, timerange=[0,0,0,0], imsize = 512, cellsize=1, robust=8, previous_image_max = 0.0, previous_image_rms = 0.0, calib_parameters = None, calib_class_name = 'CALIB', use_initial_point_source_calibration=1, image_sequence_number = 1, save_final_image = 0, save_final_calib_output=0, processing_name = '' ): """Runs a series of IMAGR/CALIB steps to self-cal a source source_name # I Used mostly for possible generation of file names # may leave blank calib_parameters # I A list of dictionaries giving the CALIB parameters. There is no gainuse variable, as this is assumed to be a single source file which has been split off and the calibration (CL) table applied. Similarly, this routine does not pay attention to a possible BP table for bandpass corrections, which it assumes have also been applied. """ assert(type(calib_parameters) == type([0])), "Need a list here" # If there are existing SN tables, warn the user that they are about # to be deleted if((__debug__) and (aips_data.table_highver('AIPS SN') > 0)): warnings.warn("'%s': SN tables exist. Hit Enter to delete."\ %self_cal_source.func_name) if(pause_before_overwrite): pause_prog('Continue?') # AIPS stupidly forces SNVER to be 0 for single-source files. So # I have to check which version is to be used myself snver_use = 0 # check that the specified flag table is present assert(check_table_version_valid(aips_data, 'AIPS FG', flagver)) # Create the IMAGR file variable imagr_image = make_selfcal_imagr_outfile(aips_data, source_name, image_sequence_number) # IMAGR will need to run on a dataset calib_out = aips_data # check for valid image max and RMS if(previous_image_max <= 0.0): previous_image_max = 1 if(previous_image_rms <= 0.0): previous_image_rms = 0 # Run through the list of SELFCAL steps for step in xrange(-1,len(calib_parameters)): if __debug__: print 'Entering step', step, 'of the SELFCAL process' # Get the right list of calibration parameters. If this is step -1 # and we are not doing a point source CALIB run, then skip if(step < 0): step = 0 if not(use_initial_point_source_calibration): continue cal_par = calib_parameters[step] if __debug__: print cal_par # Now, check whether or not we are making a new image if( not(use_initial_point_source_calibration) and (cal_par['run_imagr']) ): flagver_use = check_selfcal_flagver(calib_out, flagver) flux_level = max(math.fabs(previous_image_max*cal_par['imagr_flux_fraction']), previous_image_rms*cal_par['imagr_noise_cutoff']) run_imagr(calib_out, cellsize, imsize, source_name, imagr_image, calib_out.table_highver('AIPS SN'), cal_par['niter'], timerange, flagver_use, bpver, 'I', # Stokes flux_level, cal_par['dotv'], cal_par['clbox'], robust ) # Get the image statistics for next time (previous_image_max, junk1, junk2, previous_image_rms) \ = get_small_source_stats_main(imagr_image, imsize/4, imsize/4, imsize/2) if __debug__: # print out the CC table print_cc(imagr_image, processing_name) # End of running imagr part # set up the data files calib_in = make_selfcal_calib_data_in(aips_data, source_name, calib_class_name, cal_par['sequence_in']) if( not(calib_out == calib_in) and not(calib_out == aips_data) ): if(calib_out.exists()): calib_out.zap() calib_out = make_selfcal_calib_data_out(calib_in, source_name, calib_class_name, cal_par['sequence_out']) #print aips_data #print calib_in #print calib_out flagver_use = check_selfcal_flagver(calib_in, flagver) if not(use_initial_point_source_calibration): run_calib(calib_in, refant, cal_par['solint'], cal_par['solsub'], source_name, 0, # gainuse timerange, flagver_use, bpver, imagr_image, cal_par['ncomp'], cal_par['flux'], None, # smodel cal_par['aparm'], cal_par['soltype'], cal_par['solmode'], snver_use, 1, # save_data calib_out ) else: run_calib(calib_in, refant, cal_par['solint'], cal_par['solsub'], source_name, 0, # gainuse timerange, flagver_use, bpver, None, # image_data None, # ncomp None, # flux previous_image_max, # smodel cal_par['aparm'], cal_par['soltype'], cal_par['solmode'], snver_use, 1, # save_data calib_out ) use_initial_point_source_calibration = 0 calib_out.zap() calib_out = calib_in if(cal_par['show_SN_table']): snplot('SN', calib_in.table_highver('AIPS SN'), 'PHAS',opcode='ALSI', snuvdata=calib_in, nplots=len(calib_in.antennas)) #pause_prog('Check SN table') if(is_solmode_doing_amplitude_single(cal_par['solmode'])): snplot('SN',calib_in.table_highver('AIPS SN'),'AMP',opcode='ALSI', snuvdata=calib_in, nplots=len(calib_in.antennas)) #pause_prog('Check SN table') # end of for step in xrange over calib steps # clean up all of the CALIB files for step in xrange(0,len(calib_parameters)): calib_in =\ make_selfcal_calib_data_in_nocheck(aips_data, source_name, calib_class_name, calib_parameters[step]['sequence_in']) if( not(calib_out == calib_in) and not(calib_in == aips_data) ): if(calib_in.exists()): calib_in.zap() # Clean up after ourselves. If we are not keeping the last CALIB output # file, then delete it, unless it is the original data file for some reason. if not(save_final_calib_output): if(calib_out == aips_data): pass else: if(calib_out.exists()): calib_out.zap() if not(save_final_image): if(imagr_image.exists()): imagr_image.zap() # return some stuff to the caller, in case they want them return calib_out,imagr_image,previous_image_max,previous_image_rms ################################################################################ def process_single_source_file_to_images(split_uvdata, source_name, processing_name, robust, timerange=[0,0,0,0], refant=refant, solint_phase = 60.0, # seconds solint_amp = 2.0*60.0*60.0, # seconds imsize=None, # pixels cellsize=None, # arcsec source_size=None, # arcsec starting_sequence=1, Num_SelfCal_Iterations=1, use_initial_point_source_calibration=1, save_image_plots=0, save_output_data=0): """ make a non-self-caled image and a slef-caled image, and get statistics split_uvdata B The UV dataset for a single source file, will have SN tables when done source_name I the name of the source, which may be different from the split_uvdata name processing_name I a string to identify the information on output robust I the robust value for IMAGR timerange I the timerange to use for this imaging and SelfCal session solint_phase I time in seconds for phase-only self-cal solint_amp I the time in seconds for amplitude and phase self-cal imsize I size in pixels of the images to make (scalar) cellsize I size in arcseconds of the pixels (scalar) source_size I size of the source, as a radius, in arcseconds. This is used to calculate the mask sizes to determine the rough source sum and RMS levels starting_sequence I the sequence number of the non-selfcaled image to make. The self-caled image will be starting_sequence+1 Num_SelfCal_Iterations I How many iterations? use_initial_point_source_calibration I Do a point source selfcal first? save_image_plots I save the kntr runs as plots? save_output_data I save the map and calib files? 0 save nothing 1 save the self-calibrated image file 2 save both the uncalibrated and self-calibrated files 3 save both map files and the calibrated data file This function returns (statistics_noself, statistics, image_noself, imgdata, cal_data) as two dictionaries with the statistics and source flux density information and the AIPSImage objects for the images. Warning! the images may have been zapped if you asked for that! """ # a check for special needs if(robust > 2): if((source_name == 'J1218-0119') or (source_name == 'J1159-2228')): robust = 2 # less natural weighting to get this source to work # convert the source_size in " to pixels source_size = source_size / cellsize # Make an image with no Self-Cal image_noself = AIPSImage(source_name, 'ICL001', split_uvdata.disk, starting_sequence) image_noself = run_imagr(split_uvdata, cellsize, imsize, source_name, image_noself, 0, # gainver 350, # niter timerange, 0, # flagver, use highest -1, # bpver 'I', # Stokes 0.0, # flux level 0, # dotv [None, 10,10,imsize-10,imsize-10], # clbox robust ) (image_max, image_max_pos, image_sum, image_rms) = \ get_small_source_stats_main(image_noself, source_size, source_size + 20, source_size + 180) # Get the JMFIT statistics imaging_name = processing_name + '.' + image_noself.seq.__str__() statistics_noself = get_jmfit_statistics(image_noself, image_max_pos, imaging_name) # Add on the plain image statistics from above statistics_noself['Pixel_Max'] = image_max statistics_noself['Pixel_Sum'] = image_sum statistics_noself['Pixel_RMS'] = image_rms statistics_noself['Pixel_Pos'] = image_max_pos # Add on visibility counting info add_visibility_count_to_dict(statistics_noself, count_visibilities(split_uvdata)) print_dict_to_file(statistics_noself, imaging_name, 1, 1) # Run KNTR run_kntr2(image_noself, image_rms, image_max_pos, kntr_imsize) if(save_image_plots): run_kntr2(image_noself, image_rms, image_max_pos, kntr_imsize,dotv=-1) lwpla_last_plot(image_noself, imaging_name) #pause_prog('Done with Object. Continue?') # get the cal list cal_par = \ make_default_selfcal_calib_parameters_list(Num_SelfCal_Iterations, imsize,solint_phase, solint_amp, [imsize/2,imsize/2], source_size, 0,0,1) # do the self-cal cal_data,imgdata,image_max,image_rms,visibility_count = \ self_cal_source(split_uvdata, source_name=source_name, flagver = 0, bpver = -1, refant = refant, timerange=timerange, imsize = imsize, cellsize=cellsize, robust=robust, previous_image_max = image_max, previous_image_rms = image_rms, calib_parameters = cal_par, calib_class_name = 'CALIB', use_initial_point_source_calibration=1, image_sequence_number = starting_sequence+1, save_final_image = 1, save_final_calib_output=1, processing_name=processing_name ) # Make a nice image from this imgdata = run_imagr(cal_data, cellsize, imsize, source_name, imgdata, cal_data.table_highver('AIPS SN'), 350, timerange, 0, # flagver, use highest -1, # bpver 'I', # Stokes image_rms*0.75, # flux level 0, # dotv [None, 10,10,imsize-10,imsize-10], robust ) # Get information about the peak position and the image RMS (image_max, image_max_pos, image_sum, image_rms) = \ get_small_source_stats_main(imgdata, source_size, source_size + 20, source_size + 180) # Get the JMFIT statistics imaging_name = processing_name + '.' + imgdata.seq.__str__() statistics = get_jmfit_statistics(imgdata, image_max_pos, imaging_name) # Add on the plain image statistics from above statistics['Pixel_Max'] = image_max statistics['Pixel_Sum'] = image_sum statistics['Pixel_RMS'] = image_rms statistics['Pixel_Pos'] = image_max_pos # Add on visibility counting info add_visibility_count_to_dict(statistics_noself, count_visibilities(cal_data)) print_dict_to_file(statistics, imaging_name, 1, 1) # Run KNTR run_kntr2(imgdata, image_rms, image_max_pos, kntr_imsize) if(save_image_plots): run_kntr2(imgdata, image_rms, image_max_pos, kntr_imsize,dotv=-1) lwpla_last_plot(imgdata, imaging_name) # clean up and remove everything if(save_output_data < 3): cal_data.zap() if(save_output_data < 2): image_noself.zap() if(save_output_data < 1): imgdata.zap() # return the statistics return (statistics_noself, statistics, image_noself, imgdata, cal_data) ################################################################################ # Ok, now start imaging things def calibration_step_14(default_robust): for test_source in fring_calibrators: #for test_source in ['J1218-0119']: #for test_source in ['3C273B']: # some things to keep in mind splitclass = 'SPLIT' sequence = 1 robust = copy.copy(default_robust) solint_amp = 2.0*60.0*60.0 # the Calibrator 3C273B is whooping bright (~40 Jy), which causes # havoc with the receivers, so allow the amplitude solint to change # fast. if(test_source == '3C273B'): solint_amp = 120. # split this sucker out split_uvdata = split_data(uvdata, splitclass, sequence, test_source, Gainuse_Global, Bandpass_Global, final_flagver) # Now process the data to make an image and self-cal and make an image # I need a string to uniquely describe the process result_name = test_source + '.None.' + processing_option.__str__() \ + ".%.0f"%robust (statistics_noself, statistics, image_noself, imgdata, cal_data) = \ process_single_source_file_to_images( split_uvdata[0], test_source, result_name, robust, [0,0,0,0], # timerange refant, solint, # seconds solint_amp, # seconds imsize, # pixels cellsize, # arcsec 0.035, # source size arcsec 1, # starting sequence number Num_SelfCal_Iterations, use_initial_point_source_calibration=1, save_image_plots=1, save_output_data=0 ) # Now zap the split data split_uvdata[0].zap() if(pause_after_plots): pause_prog('Done with Object. Continue?') ################################################################################ # Get the day of year from the Year, month, day def get_day_of_year(year, month, day): """Get the day-of-year integer from the year/month/day year I YYYY month I MM starting from January == 1 day I DD starting from 1 == 1 """ day_of_year_list = [0,31,59,90,120,151,181,212,243,273,304,334] doy = day_of_year_list[month-1] + day if(month>2): if((year&0x3)==0): if((year % 100 != 0) or (year % 400 == 0)): doy = doy+1 return doy ################################################################################ # Get the day of year from the Year, month, day for the start of observations def get_observation_year_month_day(aips_data): date_string = aips_data.header.date_obs date_list = date_string.split('-') year = int(date_list[0]) month = int(date_list[1]) day = int(date_list[2]) return (year, month, day) ################################################################################ # Get the day of year from the Year, month, day for the start of observations def get_observation_day_of_year(aips_data): (year, month, day) = get_observation_year_month_day(aips_data) return get_day_of_year(year, month, day) ################################################################################ # Get the number of days of observations def get_number_days_observations(aips_data): """Find out how long the observations are. Assumes the data has been sorted """ nx_table = aips_data.table('AIPS NX', 0) num_rows = len(nx_table) # The NX table is FORTRANish # Check the first observation first_row = nx_table[0] start_time = first_row['TIME'][0] - 0.5 * first_row['TIME INTERVAL'][0] if(start_time < 0.0): raise RuntimeError, "Starting observation before starting date: Run FXTIM!" # Now check the last scan last_row = nx_table[num_rows-1] end_time = last_row['TIME'][0] + 0.5 * last_row['TIME INTERVAL'][0] return int(math.ceil(end_time)) ################################################################################ # Run a simple TECOF def ionosphere_simple_iononex(aips_data, processing_option, aparm=[None,1,0]): group_name = '' if(processing_option == 'IO_COD'): group_name='cod' elif(processing_option == 'IO_COR'): group_name='cor' elif(processing_option == 'IO_ESA'): group_name='esa' elif(processing_option == 'IO_ESR'): group_name='esr' elif(processing_option == 'IO_IGR'): group_name='igr' elif(processing_option == 'IO_IGS'): group_name='igs' elif(processing_option == 'IO_JPL'): group_name='jpl' elif(processing_option == 'IO_JPR'): group_name='jpr' elif(processing_option == 'IO_UPC'): group_name='upc' elif(processing_option == 'IO_UPR'): group_name='upr' else: raise RuntimeError, "unknown ionosphere processing option " + processing_option._str__() year = aips_data.header.date_obs[2:4] doy = get_observation_day_of_year(aips_data) filename = dirname + '/ionex/' + group_name + 'g%3.3d0.%si'%(doy,year) # Now run TECOR tecor = AIPSTask('tecor', version = aips_version_std) tecor.indata = aips_data tecor.infile = filename tecor.nfiles = get_number_days_observations(aips_data) tecor.gainver = Gainuse_Global tecor.gainuse = Gainuse_Global + 1 if(aparm != None): if(aparm[0] == None): tecor.aparm = aparm else: tecor.aparm[1:] = aparm tecor() ################################################################################ # here is the main subroutine to farm out the ionosphere options def process_ionosphere(): """find the ionosphere processiing option, and run the ionosphere stuff.""" # Get rid of previous calibration tables while uvdata.table_highver('AIPS CL') > Gainuse_Global: uvdata.zap_table('AIPS CL', 0) assert(uvdata.table_highver('AIPS CL') == Gainuse_Global) # if(processing_option[:3] == 'IO_'): ionosphere_simple_iononex(uvdata, processing_option) else: # don't know what to do raise RuntimeError, "unknown ionosphere processing option " + processing_option._str__() snplot('CL',Gainuse_Global+1,'DDLY',opcode='ALSI',snuvdata=uvdata, nplots=len(uvdata.antennas)) ################################################################################ # Here is the main processing think function def calibrate_data(runlevel): """Perform the basic calibration steps on the data Most of this function's information comes from the global variables assigned by the user at the topof this file. Other global variables are adjusted in this function as necessary based on the starting processing runlevel. """ global uvflg_file, chflg_file, antab_file, uvdata, \ Gainuse_Global, Gainuse_Before_Fring, Bandpass_Global, final_flagver # reset the calibration values # what value of the CL gain table should be globally used Gainuse_Global = 1 # what value of the bandpass table should be globally used Bandpass_Global = 0 # Pull data from the archive if(runlevel <= 1): calibration_step_1() # Set up the ancillary files uvflg_file = check_data_file('uvflg') chflg_file = check_data_file('chflag') antab_file = check_data_file('antab') # Get the initial data area uvdata = AIPSUVData(experiment.upper(), 'UVDATA', 1, 1) # Load the data into AIPS if(runlevel <= 2): calibration_step_2() # Now sort the data msortdata = AIPSUVData(uvdata.name, 'MSORT', uvdata.disk, uvdata.seq) if(runlevel <= 3): calibration_step_3(msortdata) uvdata = msortdata # Fill in some antenna information calibration_step_3_after() # A priori flagging if(runlevel <= 4): calibration_step_4() # index if(runlevel <= 5): calibration_step_5() # listr if(runlevel <= 6): calibration_step_6() # dtsum and first possm plots if(runlevel <= 7): calibration_step_7() # initial amplitude calibration if(runlevel <= 8): calibration_step_8() Gainuse_Global = Gainuse_Global + 1 # If we are applying an ionospheric correction, do it here if(processing_option != None): if(runlevel <= 9): process_ionosphere() Gainuse_Global = Gainuse_Global + 1 # paralactic angle correction if(runlevel <= 9): calibration_step_9() Gainuse_Global = Gainuse_Global + 1 # EDITR flagging # Now check that we have a valid flag table version if(final_flagver == 0): if(runlevel <= 10): calibration_step_10() final_flagver = uvdata.table_highver('AIPS FG') if not(check_table_version_valid(uvdata, 'AIPS FG', final_flagver)): while not(check_table_version_valid(uvdata, 'AIPS FG', final_flagver)): print 'Could not find FG version', final_flagver print 'Highest seen is', uvdata.table_highver('AIPS FG') print 'Trying to load new table ....' load_ascii_table(uvdata, dirname + '/' + 'flag%d.dat'%final_flagver) # Now FRING all of the FRING/AMPLITUDE calibrators if(runlevel <= 11): calibration_step_11() Gainuse_Global = Gainuse_Global + 1 # Now for the Bandpass Calibration if(runlevel <= 12): calibration_step_12() Bandpass_Global = Bandpass_Global + 1 # Now FRING all of the calibrators Gainuse_Before_Fring = Gainuse_Global if(runlevel <= 13): calibration_step_13() Gainuse_Global = Gainuse_Global + 1 # set up the cellsize information set_cellsize_info(uvdata) # Now image and self-cal all of the calibrators #pause_prog('stop before doing all of them') if(runlevel <= 14): calibration_step_14(8) if(runlevel <= 14): calibration_step_14(0) ################################################################################ # Process Phase Referencing Data def process_phase_reference_data(aips_uvdata, PhaseReferencingInfoDict): """Runs through the list of phase referencing observations to make images """ # make sure the argument is a list if(type(PhaseReferencingInfoDict) != type([0])): PhaseReferencingInfoDict = [PhaseReferencingInfoDict] # declare some local variables in case they need to be cleaned up in a # try block junk = AIPSImage('JUNKJUNKJUNKJUNK', 'JUNKCL', 1,1) split_uvdata = [junk] image_noself = junk imgdata = junk cal_data = junk del junk # Ok, what are the numbers of the gain tables which are important CL_base = Gainuse_Before_Fring CL_first_fring = aips_uvdata.table_highver('AIPS CL') +1 CL_second_fring = CL_first_fring +1 # remove any old CL tables if there was a problem if(check_table_version_valid(aips_uvdata, 'AIPS CL', CL_first_fring)): aips_uvdata.zap_table('AIPS CL', CL_first_fring) if(check_table_version_valid(aips_uvdata, 'AIPS CL', CL_second_fring)): aips_uvdata.zap_table('AIPS CL', CL_second_fring) # Clear any old SN tables if(aips_uvdata.table_highver('AIPS SN') > 0): aips_uvdata.zap_table('AIPS SN', -1) # Set up some FRING information fring_snver = 1 fring_aparm = [None, 0,0,0,0,0,0,0,0,0] fring_dparm = [None, 0,0,0,0,0,0,0,0,0] fring_aparm[1] = 3 # need at least 3 antennas for a solution fring_aparm[6] = 0 # turn most printing off fring_aparm[7] = 0 # default SNR level fring_aparm[9] = 1 # use search_list fring_dparm[1] = 2 # use small number as we only have a point-source model fring_dparm[4] = t_int fring_dparm[8] = 0 # do not zero rates or delays # define the split class splitclass = 'SPLIT' # Now loop over all pairings. for prid in PhaseReferencingInfoDict: # now, for each calibrator-target pair, run through and make # the sequence of images and statistical data try: #################################################################### # Step 1: Run Fring on the calibrator in point-source mode search_list = range(0,len(prid['Calibrator_Search_List'])) for antenna in xrange(len(search_list)): search_list[antenna] = \ antenna_dict[prid['Calibrator_Search_List'][antenna]] this_refant = search_list[0] search_list[0] = None run_fring(aips_uvdata, this_refant, search_list, prid['Calibrator_Solint_Fring'], 0, # solsub prid['Calibrator_Name'], CL_base, prid['Calibrator_Timerange'], final_flagver, Bandpass_Global, None, # image None, # ncomp None, # flux None, # smodel fring_aparm, fring_dparm, fring_snver) print 'fring done' run_clcal(aips_uvdata, CL_base, CL_first_fring, this_refant, [prid['Calibrator_Name'],prid['Target_Name']], prid['Calibrator_Name'], None, # do all timerange for these objects 'CALI', 'AMBG', snver_first=fring_snver, snver_last = fring_snver) # show the phases if(pause_after_plots): snplot('CL',CL_first_fring,'PHAS',opcode='ALSI', snuvdata=aips_uvdata,sources_use=[prid['Calibrator_Name'], prid['Target_Name']], nplots=len(aips_uvdata.antennas)) possm_2([prid['Calibrator_Name'],prid['Target_Name']], 'LL', CL_first_fring, -1, this_refant, final_flagver, bpver=Bandpass_Global) #################################################################### # Step 2: process target source to make images # split this sucker out split_uvdata = \ split_data(aips_uvdata, splitclass, 1, prid['Target_Name'], CL_first_fring, Bandpass_Global, final_flagver) # Run with robust = 8 for natural weighting robust = 8 # Now process the data to make an image and self-cal and make an # image. I need a string to uniquely describe the process result_name = prid['Target_Name'] + '.' + prid['Calibrator_Name'] + \ '.' + processing_option.__str__() + ".%.0f"%robust (statistics_noself, statistics, image_noself, imgdata, cal_data) = \ process_single_source_file_to_images( split_uvdata[0], prid['Target_Name'], result_name, robust, prid['Target_Timerange'], this_refant, prid['Target_Solint'], # seconds prid['Target_Solint_Amp'], # seconds imsize, # pixels cellsize, # arcsec prid['Target_Source_Size'], 1, # starting sequence number prid['Target_SelfCal_Iterations'], use_initial_point_source_calibration=1, save_image_plots=1, save_output_data=0 ) prid['Target_Statistics_PNSCR8'] = statistics_noself prid['Target_Statistics_PWSCR8'] = statistics # Run with robust = 0 for smaller beam size robust = 0 # Now process the data to make an image and self-cal and make an # image. I need a string to uniquely describe the process result_name = prid['Target_Name'] + '.' + prid['Calibrator_Name'] + \ '.' + processing_option.__str__() + ".%.0f"%robust (statistics_noself, statistics, image_noself, imgdata, cal_data) = \ process_single_source_file_to_images( split_uvdata[0], prid['Target_Name'], result_name, robust, prid['Target_Timerange'], this_refant, prid['Target_Solint'], # seconds prid['Target_Solint_Amp'], # seconds imsize, # pixels cellsize, # arcsec prid['Target_Source_Size'], 3, # starting sequence number prid['Target_SelfCal_Iterations'], use_initial_point_source_calibration=1, save_image_plots=1, save_output_data=0 ) prid['Target_Statistics_PNSCR0'] = statistics_noself prid['Target_Statistics_PWSCR0'] = statistics # Now zap the split data split_uvdata[0].zap() #################################################################### # Step 3: split out the calibrator to make an image to FRING with. # Then reFRING split_uvdata = \ split_data(aips_uvdata, splitclass, 1, prid['Calibrator_Name'], CL_first_fring, Bandpass_Global, final_flagver) # Run with robust = 8 for natural weighting robust = 8 # Now process the data to make an image and self-cal and make an # image. I need a string to uniquely describe the process result_name = 'tmp' (statistics_noself, statistics, image_noself, imgdata, cal_data) = \ process_single_source_file_to_images( split_uvdata[0], prid['Calibrator_Name'], result_name, robust, prid['Calibrator_Timerange'], this_refant, prid['Calibrator_Calib_Solint'], # seconds prid['Calibrator_Calib_Solint_Amp'], # seconds imsize, # pixels cellsize, # arcsec prid['Calibrator_Source_Size'], 1, # starting sequence number prid['Calibrator_SelfCal_Iterations'], use_initial_point_source_calibration=1, save_image_plots=0, save_output_data=1 # save just the self-caled image ) # zap the split data before I forget split_uvdata[0].zap() # Now fring using the new image aips_uvdata.zap_table('AIPS SN', -1) run_fring(aips_uvdata, this_refant, search_list, prid['Calibrator_Solint_Fring'], 0, # solsub prid['Calibrator_Name'], CL_base, prid['Calibrator_Timerange'], final_flagver, Bandpass_Global, imgdata, # image prid['Calibrator_Ncomp'], prid['Calibrator_Flux_Fraction']*statistics['Pixel_Max'], None, # smodel fring_aparm, fring_dparm, fring_snver) run_clcal(aips_uvdata, CL_base, CL_second_fring, this_refant, [prid['Calibrator_Name'],prid['Target_Name']], prid['Calibrator_Name'], None, # do all timerange for these objects 'CALI', 'AMBG', snver_first=fring_snver, snver_last = fring_snver) #################################################################### # Step 4: process target source to make images # split this sucker out split_uvdata = \ split_data(aips_uvdata, splitclass, 1, prid['Target_Name'], CL_second_fring, Bandpass_Global, final_flagver) # Run with robust = 8 for natural weighting robust = 8 # Now process the data to make an image and self-cal and make an # image. I need a string to uniquely describe the process result_name = prid['Target_Name'] + '.' + prid['Calibrator_Name'] + \ '.' + processing_option.__str__() + ".%.0f"%robust (statistics_noself, statistics, image_noself, imgdata, cal_data) = \ process_single_source_file_to_images( split_uvdata[0], prid['Target_Name'], result_name, robust, prid['Target_Timerange'], this_refant, prid['Target_Solint'], # seconds prid['Target_Solint_Amp'], # seconds imsize, # pixels cellsize, # arcsec prid['Target_Source_Size'], 5, # starting sequence number prid['Target_SelfCal_Iterations'], use_initial_point_source_calibration=1, save_image_plots=1, save_output_data=0 ) prid['Target_Statistics_INSCR8'] = statistics_noself prid['Target_Statistics_IWSCR8'] = statistics # Run with robust = 0 for smaller beam size robust = 0 # Now process the data to make an image and self-cal and make an # image. I need a string to uniquely describe the process result_name = prid['Target_Name'] + '.' + prid['Calibrator_Name'] + \ '.' + processing_option.__str__() + ".%.0f"%robust (statistics_noself, statistics, image_noself, imgdata, cal_data) = \ process_single_source_file_to_images( split_uvdata[0], prid['Target_Name'], result_name, robust, prid['Target_Timerange'], this_refant, prid['Target_Solint'], # seconds prid['Target_Solint_Amp'], # seconds imsize, # pixels cellsize, # arcsec prid['Target_Source_Size'], 7, # starting sequence number prid['Target_SelfCal_Iterations'], use_initial_point_source_calibration=1, save_image_plots=1, save_output_data=0 ) prid['Target_Statistics_INSCR0'] = statistics_noself prid['Target_Statistics_IWSCR0'] = statistics # Now zap the split data split_uvdata[0].zap() # Remove the new CL tables aips_uvdata.zap_table('AIPS CL', CL_first_fring) aips_uvdata.zap_table('AIPS CL', CL_second_fring) finally: if(split_uvdata[0].exists()): split_uvdata[0].zap() if(image_noself.exists()): image_noself.zap() if(imgdata.exists()): imgdata.zap() if(cal_data.exists()): cal_data.zap() ################################################################################ ################################################################################ ################################################################################ ################ MAIN CODE ##################################################### ################################################################################ ################################################################################ # This calibrates the data and sets up global variables saying what # calibration lies where calibrate_data(restart_runlevel) pause_prog('stop here!') # clean up CL tables while(uvdata.table_highver('AIPS CL') > Gainuse_Global): uvdata.zap_table('AIPS CL', uvdata.table_highver('AIPS CL')) # load in the phase referencing information import n05l2pr process_phase_reference_data(uvdata, n05l2pr.PhaseReferencingInfoDictList)