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