| 1 | import os, shutil
|
|---|
| 2 | import numpy
|
|---|
| 3 | import numpy.fft as FFT
|
|---|
| 4 | import math
|
|---|
| 5 |
|
|---|
| 6 | from asap.scantable import scantable
|
|---|
| 7 | from asap.parameters import rcParams
|
|---|
| 8 | from asap.logging import asaplog, asaplog_post_dec
|
|---|
| 9 | from asap.selector import selector
|
|---|
| 10 | from asap.asapgrid import asapgrid2
|
|---|
| 11 | from asap._asap import SBSeparator
|
|---|
| 12 |
|
|---|
| 13 | class sbseparator:
|
|---|
| 14 | """
|
|---|
| 15 | The sbseparator class is defined to separate SIGNAL and IMAGE
|
|---|
| 16 | sideband spectra observed by frequency-switching technique.
|
|---|
| 17 | It also helps supressing emmission of IMAGE sideband.
|
|---|
| 18 | *** WARNING *** THIS MODULE IS EXPERIMENTAL
|
|---|
| 19 | Known issues:
|
|---|
| 20 | - Frequencies of IMAGE sideband cannot be reconstructed from
|
|---|
| 21 | information in scantable in sideband sparation. Frequency
|
|---|
| 22 | setting of SIGNAL sideband is stored in output table for now.
|
|---|
| 23 | - Flag information (stored in FLAGTRA) is ignored.
|
|---|
| 24 |
|
|---|
| 25 | Example:
|
|---|
| 26 | # Create sideband separator instance whith 3 input data
|
|---|
| 27 | sbsep = sbseparator(['test1.asap', 'test2.asap', 'test3.asap'])
|
|---|
| 28 | # Set reference IFNO and tolerance to select data
|
|---|
| 29 | sbsep.set_frequency(5, 30, frame='TOPO')
|
|---|
| 30 | # Set direction tolerance to select data in unit of radian
|
|---|
| 31 | sbsep.set_dirtol(1.e-5)
|
|---|
| 32 | # Set rejection limit of solution
|
|---|
| 33 | sbsep.set_limit(0.2)
|
|---|
| 34 | # Solve image sideband as well
|
|---|
| 35 | sbsep.set_both(True)
|
|---|
| 36 | # Invoke sideband separation
|
|---|
| 37 | sbsep.separate('testout.asap', overwrite = True)
|
|---|
| 38 | """
|
|---|
| 39 | def __init__(self, infiles):
|
|---|
| 40 | self._separator = SBSeparator(infiles)
|
|---|
| 41 |
|
|---|
| 42 | ################## temp functions #################
|
|---|
| 43 | # def rfft(self, invec):
|
|---|
| 44 | # print "Python method cppRfft"
|
|---|
| 45 | # self._separator._cpprfft(invec)
|
|---|
| 46 | # return
|
|---|
| 47 | ################## temp functions #################
|
|---|
| 48 |
|
|---|
| 49 | def set_frequency(self, baseif, freqtol, frame=""):
|
|---|
| 50 | """
|
|---|
| 51 | Set IFNO and frequency tolerance to select data to process.
|
|---|
| 52 |
|
|---|
| 53 | Parameters:
|
|---|
| 54 | - reference IFNO to process in the first table in the list
|
|---|
| 55 | - frequency tolerance from reference IF to select data (string)
|
|---|
| 56 | frame : frequency frame to select IF
|
|---|
| 57 | """
|
|---|
| 58 | if type(freqtol) in (float, int):
|
|---|
| 59 | freqtol = str(freqtol)
|
|---|
| 60 | elif isinstance(freqtol, dict):
|
|---|
| 61 | try:
|
|---|
| 62 | freqtol = str(freqtol['value']) + freqtol['unit']
|
|---|
| 63 | except:
|
|---|
| 64 | raise ValueError, "Invalid frequency tolerance."
|
|---|
| 65 | self._separator.set_freq(baseif, freqtol, frame)
|
|---|
| 66 |
|
|---|
| 67 |
|
|---|
| 68 | def set_dirtol(self, dirtol=["2arcsec", "2arcsec"]):
|
|---|
| 69 | """
|
|---|
| 70 | Set tolerance of direction to select data
|
|---|
| 71 | """
|
|---|
| 72 | if isinstance(dirtol, str):
|
|---|
| 73 | dirtol = [dirtol]
|
|---|
| 74 |
|
|---|
| 75 | self._separator.set_dirtol(dirtol)
|
|---|
| 76 |
|
|---|
| 77 |
|
|---|
| 78 | def set_shift(self, imageshift=[]):
|
|---|
| 79 | """
|
|---|
| 80 | Set shift mode and channel shift of image band.
|
|---|
| 81 |
|
|---|
| 82 | imageshift : a list of number of channels shifted in image
|
|---|
| 83 | side band of each scantable.
|
|---|
| 84 | If the shifts are not set, they are assumed to be
|
|---|
| 85 | equal to those of signal side band, but in opposite
|
|---|
| 86 | direction as usual by LO1 offsetting of DSB receivers.
|
|---|
| 87 | """
|
|---|
| 88 | if not imageshift:
|
|---|
| 89 | imageshift = []
|
|---|
| 90 | self._separator.set_shift(imageshift)
|
|---|
| 91 |
|
|---|
| 92 |
|
|---|
| 93 | @asaplog_post_dec
|
|---|
| 94 | def set_both(self, flag=False):
|
|---|
| 95 | """
|
|---|
| 96 | Resolve both image and signal sideband when True.
|
|---|
| 97 | """
|
|---|
| 98 | self._separator.solve_both(flag)
|
|---|
| 99 | if flag:
|
|---|
| 100 | asaplog.push("Both signal and image sidebands will be solved and stored in separate tables.")
|
|---|
| 101 | else:
|
|---|
| 102 | asaplog.push("Only signal sideband will be solved and stored in an table.")
|
|---|
| 103 |
|
|---|
| 104 | @asaplog_post_dec
|
|---|
| 105 | def set_limit(self, threshold=0.2):
|
|---|
| 106 | """
|
|---|
| 107 | Set rejection limit of solution.
|
|---|
| 108 | """
|
|---|
| 109 | self._separator.set_limit(threshold)
|
|---|
| 110 |
|
|---|
| 111 |
|
|---|
| 112 | @asaplog_post_dec
|
|---|
| 113 | def set_solve_other(self, flag=False):
|
|---|
| 114 | """
|
|---|
| 115 | Calculate spectra by subtracting the solution of the other sideband
|
|---|
| 116 | when True.
|
|---|
| 117 | """
|
|---|
| 118 | self._separator.subtract_other(flag)
|
|---|
| 119 | if flag:
|
|---|
| 120 | asaplog.push("Expert mode: solution are obtained by subtraction of the other sideband.")
|
|---|
| 121 |
|
|---|
| 122 |
|
|---|
| 123 | def set_lo1(self,lo1, frame="TOPO", reftime=-1, refdir=""):
|
|---|
| 124 | """
|
|---|
| 125 | Set LO1 frequency to calculate frequency of image sideband.
|
|---|
| 126 |
|
|---|
| 127 | lo1 : LO1 frequency
|
|---|
| 128 | frame : the frequency frame of LO1
|
|---|
| 129 | reftime : the reference time used in frequency frame conversion.
|
|---|
| 130 | refdir : the reference direction used in frequency frame conversion.
|
|---|
| 131 | """
|
|---|
| 132 | self._separator.set_lo1(lo1val, frame, reftime, refdir)
|
|---|
| 133 |
|
|---|
| 134 |
|
|---|
| 135 | def set_lo1root(self, name):
|
|---|
| 136 | """
|
|---|
| 137 | Set MS name which stores LO1 frequency of signal side band.
|
|---|
| 138 | It is used to calculate frequency of image sideband.
|
|---|
| 139 |
|
|---|
| 140 | name : MS name which contains 'ASDM_SPECTRALWINDOW' and
|
|---|
| 141 | 'ASDM_RECEIVER' tables.
|
|---|
| 142 | """
|
|---|
| 143 | self._separator.set_lo1root(name)
|
|---|
| 144 |
|
|---|
| 145 | @asaplog_post_dec
|
|---|
| 146 | def separate(self, outname="", overwrite=False):
|
|---|
| 147 | """
|
|---|
| 148 | Invoke sideband separation.
|
|---|
| 149 |
|
|---|
| 150 | outname : a name of output scantable
|
|---|
| 151 | overwrite : overwrite existing table
|
|---|
| 152 | """
|
|---|
| 153 | out_default = "sbseparated.asap"
|
|---|
| 154 | if len(outname) == 0:
|
|---|
| 155 | outname = out_default
|
|---|
| 156 | asaplog.post()
|
|---|
| 157 | asaplog.push("The output file name is not specified.")
|
|---|
| 158 | asaplog.push("Using default name '%s'" % outname)
|
|---|
| 159 | asaplog.post("WARN")
|
|---|
| 160 |
|
|---|
| 161 | if os.path.exists(outname):
|
|---|
| 162 | if overwrite:
|
|---|
| 163 | asaplog.push("removing the old file '%s'" % outname)
|
|---|
| 164 | shutil.rmtree(outname)
|
|---|
| 165 | else:
|
|---|
| 166 | asaplog.post()
|
|---|
| 167 | asaplog.push("Output file '%s' already exists." % outname)
|
|---|
| 168 | asaplog.post("ERROR")
|
|---|
| 169 | return False
|
|---|
| 170 |
|
|---|
| 171 | self._separator.separate(outname)
|
|---|
| 172 |
|
|---|
| 173 | ######################################################################
|
|---|
| 174 | # @asaplog_post_dec
|
|---|
| 175 | # def separate0(self, outname="", overwrite=False):
|
|---|
| 176 | # """
|
|---|
| 177 | # Invoke sideband separation.
|
|---|
| 178 |
|
|---|
| 179 | # outname : a name of output scantable
|
|---|
| 180 | # overwrite : overwrite existing table
|
|---|
| 181 | # """
|
|---|
| 182 | # # List up valid scantables and IFNOs to convolve.
|
|---|
| 183 | # #self._separator.separate()
|
|---|
| 184 | # self._setup_shift()
|
|---|
| 185 | # #self._preprocess_tables()
|
|---|
| 186 |
|
|---|
| 187 | # ### TEMPORAL ###
|
|---|
| 188 | # self._separator._get_asistb_from_scantb(self.tables[0])
|
|---|
| 189 | # ################
|
|---|
| 190 |
|
|---|
| 191 | # nshift = len(self.tables)
|
|---|
| 192 | # signaltab = self._grid_outtable(self.tables[0])
|
|---|
| 193 | # if self.getboth:
|
|---|
| 194 | # imagetab = signaltab.copy()
|
|---|
| 195 |
|
|---|
| 196 | # rejrow = []
|
|---|
| 197 | # for irow in xrange(signaltab.nrow()):
|
|---|
| 198 | # currpol = signaltab.getpol(irow)
|
|---|
| 199 | # currbeam = signaltab.getbeam(irow)
|
|---|
| 200 | # currdir = signaltab.get_directionval(irow)
|
|---|
| 201 | # spec_array, tabidx = self._get_specarray(polid=currpol,\
|
|---|
| 202 | # beamid=currbeam,\
|
|---|
| 203 | # dir=currdir)
|
|---|
| 204 | # #if not spec_array:
|
|---|
| 205 | # if len(tabidx) == 0:
|
|---|
| 206 | # asaplog.post()
|
|---|
| 207 | # asaplog.push("skipping row %d" % irow)
|
|---|
| 208 | # rejrow.append(irow)
|
|---|
| 209 | # continue
|
|---|
| 210 | # signal = self._solve_signal(spec_array, tabidx)
|
|---|
| 211 | # signaltab.set_spectrum(signal, irow)
|
|---|
| 212 |
|
|---|
| 213 | # # Solve image side side band
|
|---|
| 214 | # if self.getboth:
|
|---|
| 215 | # image = self._solve_image(spec_array, tabidx)
|
|---|
| 216 | # imagetab.set_spectrum(image, irow)
|
|---|
| 217 |
|
|---|
| 218 | # # TODO: Need to remove rejrow form scantables here
|
|---|
| 219 | # signaltab.flag_row(rejrow)
|
|---|
| 220 | # if self.getboth:
|
|---|
| 221 | # imagetab.flag_row(rejrow)
|
|---|
| 222 |
|
|---|
| 223 | # if outname == "":
|
|---|
| 224 | # outname = "sbsepareted.asap"
|
|---|
| 225 | # signalname = outname + ".signalband"
|
|---|
| 226 | # if os.path.exists(signalname):
|
|---|
| 227 | # if not overwrite:
|
|---|
| 228 | # raise Exception, "Output file '%s' exists." % signalname
|
|---|
| 229 | # else:
|
|---|
| 230 | # shutil.rmtree(signalname)
|
|---|
| 231 | # signaltab.save(signalname)
|
|---|
| 232 |
|
|---|
| 233 | # if self.getboth:
|
|---|
| 234 | # # Warnings
|
|---|
| 235 | # asaplog.post()
|
|---|
| 236 | # asaplog.push("Saving IMAGE sideband.")
|
|---|
| 237 | # #asaplog.push("Note, frequency information of IMAGE sideband cannot be properly filled so far. (future development)")
|
|---|
| 238 | # #asaplog.push("Storing frequency setting of SIGNAL sideband in FREQUENCIES table for now.")
|
|---|
| 239 | # #asaplog.post("WARN")
|
|---|
| 240 |
|
|---|
| 241 | # imagename = outname + ".imageband"
|
|---|
| 242 | # if os.path.exists(imagename):
|
|---|
| 243 | # if not overwrite:
|
|---|
| 244 | # raise Exception, "Output file '%s' exists." % imagename
|
|---|
| 245 | # else:
|
|---|
| 246 | # shutil.rmtree(imagename)
|
|---|
| 247 | # # Update frequency information
|
|---|
| 248 | # self._separator.set_imgtable(imagetab)
|
|---|
| 249 | # self._separator.solve_imgfreq()
|
|---|
| 250 | # imagetab.save(imagename)
|
|---|
| 251 |
|
|---|
| 252 |
|
|---|
| 253 | # def _solve_signal(self, data, tabidx=None):
|
|---|
| 254 | # if not tabidx:
|
|---|
| 255 | # tabidx = range(len(data))
|
|---|
| 256 |
|
|---|
| 257 | # tempshift = []
|
|---|
| 258 | # dshift = []
|
|---|
| 259 | # if self.solveother:
|
|---|
| 260 | # for idx in tabidx:
|
|---|
| 261 | # tempshift.append(-self.imageShift[idx])
|
|---|
| 262 | # dshift.append(self.signalShift[idx] - self.imageShift[idx])
|
|---|
| 263 | # else:
|
|---|
| 264 | # for idx in tabidx:
|
|---|
| 265 | # tempshift.append(-self.signalShift[idx])
|
|---|
| 266 | # dshift.append(self.imageShift[idx] - self.signalShift[idx])
|
|---|
| 267 |
|
|---|
| 268 | # shiftdata = numpy.zeros(data.shape, numpy.float)
|
|---|
| 269 | # for i in range(len(data)):
|
|---|
| 270 | # shiftdata[i] = self._shiftSpectrum(data[i], tempshift[i])
|
|---|
| 271 | # ifftdata = self._Deconvolution(shiftdata, dshift, self.rejlimit)
|
|---|
| 272 | # result_image = self._combineResult(ifftdata)
|
|---|
| 273 | # if not self.solveother:
|
|---|
| 274 | # return result_image
|
|---|
| 275 | # result_signal = self._subtractOtherSide(shiftdata, dshift, result_image)
|
|---|
| 276 | # return result_signal
|
|---|
| 277 |
|
|---|
| 278 |
|
|---|
| 279 | # def _solve_image(self, data, tabidx=None):
|
|---|
| 280 | # if not tabidx:
|
|---|
| 281 | # tabidx = range(len(data))
|
|---|
| 282 |
|
|---|
| 283 | # tempshift = []
|
|---|
| 284 | # dshift = []
|
|---|
| 285 | # if self.solveother:
|
|---|
| 286 | # for idx in tabidx:
|
|---|
| 287 | # tempshift.append(-self.signalShift[idx])
|
|---|
| 288 | # dshift.append(self.imageShift[idx] - self.signalShift[idx])
|
|---|
| 289 | # else:
|
|---|
| 290 | # for idx in tabidx:
|
|---|
| 291 | # tempshift.append(-self.imageShift[idx])
|
|---|
| 292 | # dshift.append(self.signalShift[idx] - self.imageShift[idx])
|
|---|
| 293 |
|
|---|
| 294 | # shiftdata = numpy.zeros(data.shape, numpy.float)
|
|---|
| 295 | # for i in range(len(data)):
|
|---|
| 296 | # shiftdata[i] = self._shiftSpectrum(data[i], tempshift[i])
|
|---|
| 297 | # ifftdata = self._Deconvolution(shiftdata, dshift, self.rejlimit)
|
|---|
| 298 | # result_image = self._combineResult(ifftdata)
|
|---|
| 299 | # if not self.solveother:
|
|---|
| 300 | # return result_image
|
|---|
| 301 | # result_signal = self._subtractOtherSide(shiftdata, dshift, result_image)
|
|---|
| 302 | # return result_signal
|
|---|
| 303 |
|
|---|
| 304 | # @asaplog_post_dec
|
|---|
| 305 | # def _grid_outtable(self, table):
|
|---|
| 306 | # # Generate gridded table for output (Just to get rows)
|
|---|
| 307 | # gridder = asapgrid2(table)
|
|---|
| 308 | # gridder.setIF(self.baseif)
|
|---|
| 309 |
|
|---|
| 310 | # cellx = str(self.dirtol[0])+"rad"
|
|---|
| 311 | # celly = str(self.dirtol[1])+"rad"
|
|---|
| 312 | # dirarr = numpy.array(table.get_directionval()).transpose()
|
|---|
| 313 | # mapx = dirarr[0].max() - dirarr[0].min()
|
|---|
| 314 | # mapy = dirarr[1].max() - dirarr[1].min()
|
|---|
| 315 | # centy = 0.5 * (dirarr[1].max() + dirarr[1].min())
|
|---|
| 316 | # nx = max(1, numpy.ceil(mapx*numpy.cos(centy)/self.dirtol[0]))
|
|---|
| 317 | # ny = max(1, numpy.ceil(mapy/self.dirtol[0]))
|
|---|
| 318 |
|
|---|
| 319 | # asaplog.push("Regrid output scantable with cell = [%s, %s]" % \
|
|---|
| 320 | # (cellx, celly))
|
|---|
| 321 | # gridder.defineImage(nx=nx, ny=ny, cellx=cellx, celly=celly)
|
|---|
| 322 | # gridder.setFunc(func='box', convsupport=1)
|
|---|
| 323 | # gridder.setWeight(weightType='uniform')
|
|---|
| 324 | # gridder.grid()
|
|---|
| 325 | # return gridder.getResult()
|
|---|
| 326 |
|
|---|
| 327 | # @asaplog_post_dec
|
|---|
| 328 | # def _get_specarray(self, polid=None, beamid=None, dir=None):
|
|---|
| 329 | # ntable = len(self.tables)
|
|---|
| 330 | # spec_array = numpy.zeros((ntable, self.nchan), numpy.float)
|
|---|
| 331 | # nspec = 0
|
|---|
| 332 | # asaplog.push("Start data selection by POL=%d, BEAM=%d, direction=[%f, %f]" % (polid, beamid, dir[0], dir[1]))
|
|---|
| 333 | # tabidx = []
|
|---|
| 334 | # for itab in range(ntable):
|
|---|
| 335 | # tab = self.tables[itab]
|
|---|
| 336 | # # Select rows by POLNO and BEAMNO
|
|---|
| 337 | # try:
|
|---|
| 338 | # tab.set_selection(pols=[polid], beams=[beamid])
|
|---|
| 339 | # if tab.nrow() > 0: tabidx.append(itab)
|
|---|
| 340 | # except: # no selection
|
|---|
| 341 | # asaplog.post()
|
|---|
| 342 | # asaplog.push("table %d - No spectrum ....skipping the table" % (itab))
|
|---|
| 343 | # asaplog.post("WARN")
|
|---|
| 344 | # continue
|
|---|
| 345 |
|
|---|
| 346 | # # Select rows by direction
|
|---|
| 347 | # spec = numpy.zeros(self.nchan, numpy.float)
|
|---|
| 348 | # selrow = []
|
|---|
| 349 | # for irow in range(tab.nrow()):
|
|---|
| 350 | # currdir = tab.get_directionval(irow)
|
|---|
| 351 | # if (abs(currdir[0] - dir[0]) > self.dirtol[0]) or \
|
|---|
| 352 | # (abs(currdir[1] - dir[1]) > self.dirtol[1]):
|
|---|
| 353 | # continue
|
|---|
| 354 | # selrow.append(irow)
|
|---|
| 355 | # if len(selrow) == 0:
|
|---|
| 356 | # asaplog.post()
|
|---|
| 357 | # asaplog.push("table %d - No spectrum ....skipping the table" % (itab))
|
|---|
| 358 | # asaplog.post("WARN")
|
|---|
| 359 | # continue
|
|---|
| 360 | # else:
|
|---|
| 361 | # seltab = tab.copy()
|
|---|
| 362 | # seltab.set_selection(selector(rows=selrow))
|
|---|
| 363 |
|
|---|
| 364 | # if seltab.nrow() > 1:
|
|---|
| 365 | # asaplog.push("table %d - More than a spectrum selected. averaging rows..." % (itab))
|
|---|
| 366 | # tab = seltab.average_time(scanav=False, weight="tintsys")
|
|---|
| 367 | # else:
|
|---|
| 368 | # tab = seltab
|
|---|
| 369 |
|
|---|
| 370 | # spec_array[nspec] = tab._getspectrum()
|
|---|
| 371 | # nspec += 1
|
|---|
| 372 |
|
|---|
| 373 | # if nspec != ntable:
|
|---|
| 374 | # asaplog.post()
|
|---|
| 375 | # #asaplog.push("Some tables has no spectrum with POL=%d BEAM=%d. averaging rows..." % (polid, beamid))
|
|---|
| 376 | # asaplog.push("Could not find corresponding rows in some tables.")
|
|---|
| 377 | # asaplog.push("Number of spectra selected = %d (table: %d)" % (nspec, ntable))
|
|---|
| 378 | # if nspec < 2:
|
|---|
| 379 | # asaplog.push("At least 2 spectra are necessary for convolution")
|
|---|
| 380 | # asaplog.post("ERROR")
|
|---|
| 381 | # return False, tabidx
|
|---|
| 382 |
|
|---|
| 383 | # return spec_array[0:nspec], tabidx
|
|---|
| 384 |
|
|---|
| 385 |
|
|---|
| 386 | # @asaplog_post_dec
|
|---|
| 387 | # def _setup_shift(self):
|
|---|
| 388 | # ### define self.tables, self.signalShift, and self.imageShift
|
|---|
| 389 | # if not self.intables:
|
|---|
| 390 | # asaplog.post()
|
|---|
| 391 | # raise RuntimeError, "Input data is not defined."
|
|---|
| 392 | # #if self.baseif < 0:
|
|---|
| 393 | # # asaplog.post()
|
|---|
| 394 | # # raise RuntimeError, "Reference IFNO is not defined."
|
|---|
| 395 |
|
|---|
| 396 | # byname = False
|
|---|
| 397 | # #if not self.intables:
|
|---|
| 398 | # if isinstance(self.intables[0], str):
|
|---|
| 399 | # # A list of file name is given
|
|---|
| 400 | # if not os.path.exists(self.intables[0]):
|
|---|
| 401 | # asaplog.post()
|
|---|
| 402 | # raise RuntimeError, "Could not find '%s'" % self.intables[0]
|
|---|
| 403 |
|
|---|
| 404 | # stab = scantable(self.intables[0],average=False)
|
|---|
| 405 | # ntab = len(self.intables)
|
|---|
| 406 | # byname = True
|
|---|
| 407 | # else:
|
|---|
| 408 | # stab = self.intables[0]
|
|---|
| 409 | # ntab = len(self.intables)
|
|---|
| 410 |
|
|---|
| 411 | # if len(stab.getbeamnos()) > 1:
|
|---|
| 412 | # asaplog.post()
|
|---|
| 413 | # asaplog.push("Mult-beam data is not supported by this module.")
|
|---|
| 414 | # asaplog.post("ERROR")
|
|---|
| 415 | # return
|
|---|
| 416 |
|
|---|
| 417 | # valid_ifs = stab.getifnos()
|
|---|
| 418 | # if self.baseif < 0:
|
|---|
| 419 | # self.baseif = valid_ifs[0]
|
|---|
| 420 | # asaplog.post()
|
|---|
| 421 | # asaplog.push("IFNO is not selected. Using the first IF in the first scantable. Reference IFNO = %d" % (self.baseif))
|
|---|
| 422 |
|
|---|
| 423 | # if not (self.baseif in valid_ifs):
|
|---|
| 424 | # asaplog.post()
|
|---|
| 425 | # errmsg = "IF%d does not exist in the first scantable" % \
|
|---|
| 426 | # self.baseif
|
|---|
| 427 | # raise RuntimeError, errmsg
|
|---|
| 428 |
|
|---|
| 429 | # asaplog.push("Start selecting tables and IFNOs to solve.")
|
|---|
| 430 | # asaplog.push("Checking frequency of the reference IF")
|
|---|
| 431 | # unit_org = stab.get_unit()
|
|---|
| 432 | # coord = stab._getcoordinfo()
|
|---|
| 433 | # frame_org = coord[1]
|
|---|
| 434 | # stab.set_unit("Hz")
|
|---|
| 435 | # if len(self.freqframe) > 0:
|
|---|
| 436 | # stab.set_freqframe(self.freqframe)
|
|---|
| 437 | # stab.set_selection(ifs=[self.baseif])
|
|---|
| 438 | # spx = stab._getabcissa()
|
|---|
| 439 | # stab.set_selection()
|
|---|
| 440 | # basech0 = spx[0]
|
|---|
| 441 | # baseinc = spx[1]-spx[0]
|
|---|
| 442 | # self.nchan = len(spx)
|
|---|
| 443 | # # frequency tolerance
|
|---|
| 444 | # if isinstance(self.freqtol, dict) and self.freqtol['unit'] == "Hz":
|
|---|
| 445 | # vftol = abs(self.freqtol['value'])
|
|---|
| 446 | # else:
|
|---|
| 447 | # vftol = abs(baseinc * float(self.freqtol))
|
|---|
| 448 | # self.freqtol = dict(value=vftol, unit="Hz")
|
|---|
| 449 | # # tolerance of frequency increment
|
|---|
| 450 | # inctol = abs(baseinc/float(self.nchan))
|
|---|
| 451 | # asaplog.push("Reference frequency setup (Table = 0, IFNO = %d): nchan = %d, chan0 = %f Hz, incr = %f Hz" % (self.baseif, self.nchan, basech0, baseinc))
|
|---|
| 452 | # asaplog.push("Allowed frequency tolerance = %f Hz ( %f channels)" % (vftol, vftol/baseinc))
|
|---|
| 453 | # poltype0 = stab.poltype()
|
|---|
| 454 |
|
|---|
| 455 | # self.tables = []
|
|---|
| 456 | # self.signalShift = []
|
|---|
| 457 | # if self.dsbmode:
|
|---|
| 458 | # self.imageShift = []
|
|---|
| 459 |
|
|---|
| 460 | # for itab in range(ntab):
|
|---|
| 461 | # asaplog.push("Table %d:" % itab)
|
|---|
| 462 | # tab_selected = False
|
|---|
| 463 | # if itab > 0:
|
|---|
| 464 | # if byname:
|
|---|
| 465 | # stab = scantable(self.intables[itab],average=False)
|
|---|
| 466 | # else:
|
|---|
| 467 | # stab = self.intables[itab]
|
|---|
| 468 | # unit_org = stab.get_unit()
|
|---|
| 469 | # coord = stab._getcoordinfo()
|
|---|
| 470 | # frame_org = coord[1]
|
|---|
| 471 | # stab.set_unit("Hz")
|
|---|
| 472 | # if len(self.freqframe) > 0:
|
|---|
| 473 | # stab.set_freqframe(self.freqframe)
|
|---|
| 474 |
|
|---|
| 475 | # # Check POLTYPE should be equal to the first table.
|
|---|
| 476 | # if stab.poltype() != poltype0:
|
|---|
| 477 | # asaplog.post()
|
|---|
| 478 | # raise Exception, "POLTYPE should be equal to the first table."
|
|---|
| 479 | # # Multiple beam data may not handled properly
|
|---|
| 480 | # if len(stab.getbeamnos()) > 1:
|
|---|
| 481 | # asaplog.post()
|
|---|
| 482 | # asaplog.push("table contains multiple beams. It may not be handled properly.")
|
|---|
| 483 | # asaplog.push("WARN")
|
|---|
| 484 |
|
|---|
| 485 | # for ifno in stab.getifnos():
|
|---|
| 486 | # stab.set_selection(ifs=[ifno])
|
|---|
| 487 | # spx = stab._getabcissa()
|
|---|
| 488 | # if (len(spx) != self.nchan) or \
|
|---|
| 489 | # (abs(spx[0]-basech0) > vftol) or \
|
|---|
| 490 | # (abs(spx[1]-spx[0]-baseinc) > inctol):
|
|---|
| 491 | # continue
|
|---|
| 492 | # tab_selected = True
|
|---|
| 493 | # seltab = stab.copy()
|
|---|
| 494 | # seltab.set_unit(unit_org)
|
|---|
| 495 | # seltab.set_freqframe(frame_org)
|
|---|
| 496 | # self.tables.append(seltab)
|
|---|
| 497 | # self.signalShift.append((spx[0]-basech0)/baseinc)
|
|---|
| 498 | # if self.dsbmode:
|
|---|
| 499 | # self.imageShift.append(-self.signalShift[-1])
|
|---|
| 500 | # asaplog.push("- IF%d selected: sideband shift = %16.12e channels" % (ifno, self.signalShift[-1]))
|
|---|
| 501 | # stab.set_selection()
|
|---|
| 502 | # stab.set_unit(unit_org)
|
|---|
| 503 | # stab.set_freqframe(frame_org)
|
|---|
| 504 | # if not tab_selected:
|
|---|
| 505 | # asaplog.post()
|
|---|
| 506 | # asaplog.push("No data selected in Table %d" % itab)
|
|---|
| 507 | # asaplog.post("WARN")
|
|---|
| 508 |
|
|---|
| 509 | # asaplog.push("Total number of IFs selected = %d" % len(self.tables))
|
|---|
| 510 | # if len(self.tables) < 2:
|
|---|
| 511 | # asaplog.post()
|
|---|
| 512 | # raise RuntimeError, "At least 2 IFs are necessary for convolution!"
|
|---|
| 513 |
|
|---|
| 514 | # if not self.dsbmode and len(self.imageShift) != len(self.signalShift):
|
|---|
| 515 | # asaplog.post()
|
|---|
| 516 | # errmsg = "User defined channel shift of image sideband has %d elements, while selected IFNOs are %d" % (len(self.imageShift), len(self.signalShift))
|
|---|
| 517 | # errmsg += "\nThe frequency tolerance (freqtol) you set may be too small."
|
|---|
| 518 | # raise RuntimeError, errmsg
|
|---|
| 519 |
|
|---|
| 520 | # self.signalShift = numpy.array(self.signalShift)
|
|---|
| 521 | # self.imageShift = numpy.array(self.imageShift)
|
|---|
| 522 | # self.nshift = len(self.tables)
|
|---|
| 523 |
|
|---|
| 524 | # @asaplog_post_dec
|
|---|
| 525 | # def _preprocess_tables(self):
|
|---|
| 526 | # ### temporary method to preprocess data
|
|---|
| 527 | # ### Do time averaging for now.
|
|---|
| 528 | # for itab in range(len(self.tables)):
|
|---|
| 529 | # self.tables[itab] = self.tables[itab].average_time(scanav=False, weight="tintsys")
|
|---|
| 530 |
|
|---|
| 531 |
|
|---|
| 532 | # ########################################################################
|
|---|
| 533 | # def _Deconvolution(self, data_array, shift, threshold=0.00000001):
|
|---|
| 534 | # FObs = []
|
|---|
| 535 | # Reject = 0
|
|---|
| 536 | # nshift, nchan = data_array.shape
|
|---|
| 537 | # nspec = nshift*(nshift-1)/2
|
|---|
| 538 | # ifftObs = numpy.zeros((nspec, nchan), numpy.float)
|
|---|
| 539 | # for i in range(nshift):
|
|---|
| 540 | # F = FFT.fft(data_array[i])
|
|---|
| 541 | # FObs.append(F)
|
|---|
| 542 | # z = 0
|
|---|
| 543 | # for i in range(nshift):
|
|---|
| 544 | # for j in range(i+1, nshift):
|
|---|
| 545 | # Fobs = (FObs[i]+FObs[j])/2.0
|
|---|
| 546 | # dX = (shift[j]-shift[i])*2.0*math.pi/float(self.nchan)
|
|---|
| 547 | # #print 'dX,i,j=',dX,i,j
|
|---|
| 548 | # for k in range(1,self.nchan):
|
|---|
| 549 | # if math.fabs(math.sin(dX*k)) > threshold:
|
|---|
| 550 | # Fobs[k] += ((FObs[i][k]-FObs[j][k])/2.0/(1.0-math.cos(dX*k))*math.sin(dX*k))*1.0j
|
|---|
| 551 | # else: Reject += 1
|
|---|
| 552 | # ifftObs[z] = FFT.ifft(Fobs)
|
|---|
| 553 | # z += 1
|
|---|
| 554 | # print 'Threshold=%s Reject=%d' % (threshold, Reject)
|
|---|
| 555 | # return ifftObs
|
|---|
| 556 |
|
|---|
| 557 | # def _combineResult(self, ifftObs):
|
|---|
| 558 | # nspec = len(ifftObs)
|
|---|
| 559 | # sum = ifftObs[0]
|
|---|
| 560 | # for i in range(1,nspec):
|
|---|
| 561 | # sum += ifftObs[i]
|
|---|
| 562 | # return(sum/float(nspec))
|
|---|
| 563 |
|
|---|
| 564 | # def _subtractOtherSide(self, data_array, shift, Data):
|
|---|
| 565 | # sum = numpy.zeros(len(Data), numpy.float)
|
|---|
| 566 | # numSP = len(data_array)
|
|---|
| 567 | # for i in range(numSP):
|
|---|
| 568 | # SPsub = data_array[i] - Data
|
|---|
| 569 | # sum += self._shiftSpectrum(SPsub, -shift[i])
|
|---|
| 570 | # return(sum/float(numSP))
|
|---|
| 571 |
|
|---|
| 572 | # def _shiftSpectrum(self, data, Shift):
|
|---|
| 573 | # Out = numpy.zeros(self.nchan, numpy.float)
|
|---|
| 574 | # w2 = Shift % 1
|
|---|
| 575 | # w1 = 1.0 - w2
|
|---|
| 576 | # for i in range(self.nchan):
|
|---|
| 577 | # c1 = int((Shift + i) % self.nchan)
|
|---|
| 578 | # c2 = (c1 + 1) % self.nchan
|
|---|
| 579 | # Out[c1] += data[i] * w1
|
|---|
| 580 | # Out[c2] += data[i] * w2
|
|---|
| 581 | # return Out.copy()
|
|---|
