source: branches/alma/python/asapmath.py@ 1631

Last change on this file since 1631 was 1631, checked in by Takeshi Nakazato, 15 years ago

New Development: No

JIRA Issue: Yes CAS-1424

Ready to Release: Yes

Interface Changes: Yes

What Interface Changed: Added 'verify' parameter for asapmath.calps, asapmath.calnod, asapmath.calfs, and scantable.auto_quotient

Test Programs: do calibration using sdaverage with verify=True

Put in Release Notes: Yes

Module(s): task sdaverage and sdcal

Description: Describe your changes here...

I have added verify parameter to verify calibration result.
If verify=True, methods plot pre- and post-calibration spectra and asks user
to be acceptable calibration or not.
If not, calibration will not be applied.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 30.1 KB
Line 
1from asap.scantable import scantable
2from asap import rcParams
3from asap import print_log
4from asap import selector
5from asap import asaplog
6from asap import asaplotgui
7
8def average_time(*args, **kwargs):
9 """
10 Return the (time) average of a scan or list of scans. [in channels only]
11 The cursor of the output scan is set to 0
12 Parameters:
13 one scan or comma separated scans or a list of scans
14 mask: an optional mask (only used for 'var' and 'tsys' weighting)
15 scanav: True averages each scan separately.
16 False (default) averages all scans together,
17 weight: Weighting scheme.
18 'none' (mean no weight)
19 'var' (1/var(spec) weighted)
20 'tsys' (1/Tsys**2 weighted)
21 'tint' (integration time weighted)
22 'tintsys' (Tint/Tsys**2)
23 'median' ( median averaging)
24 align: align the spectra in velocity before averaging. It takes
25 the time of the first spectrum in the first scantable
26 as reference time.
27 Example:
28 # return a time averaged scan from scana and scanb
29 # without using a mask
30 scanav = average_time(scana,scanb)
31 # or equivalent
32 # scanav = average_time([scana, scanb])
33 # return the (time) averaged scan, i.e. the average of
34 # all correlator cycles
35 scanav = average_time(scan, scanav=True)
36 """
37 scanav = False
38 if kwargs.has_key('scanav'):
39 scanav = kwargs.get('scanav')
40 weight = 'tint'
41 if kwargs.has_key('weight'):
42 weight = kwargs.get('weight')
43 mask = ()
44 if kwargs.has_key('mask'):
45 mask = kwargs.get('mask')
46 align = False
47 if kwargs.has_key('align'):
48 align = kwargs.get('align')
49 compel = False
50 if kwargs.has_key('compel'):
51 compel = kwargs.get('compel')
52 varlist = vars()
53 if isinstance(args[0],list):
54 lst = args[0]
55 elif isinstance(args[0],tuple):
56 lst = list(args[0])
57 else:
58 lst = list(args)
59
60 del varlist["kwargs"]
61 varlist["args"] = "%d scantables" % len(lst)
62 # need special formatting here for history...
63
64 from asap._asap import stmath
65 stm = stmath()
66 for s in lst:
67 if not isinstance(s,scantable):
68 msg = "Please give a list of scantables"
69 if rcParams['verbose']:
70 #print msg
71 asaplog.push(msg)
72 print_log('ERROR')
73 return
74 else:
75 raise TypeError(msg)
76 if scanav: scanav = "SCAN"
77 else: scanav = "NONE"
78 alignedlst = []
79 if align:
80 refepoch = lst[0].get_time(0)
81 for scan in lst:
82 alignedlst.append(scan.freq_align(refepoch,insitu=False))
83 else:
84 alignedlst = lst
85 if weight.upper() == 'MEDIAN':
86 # median doesn't support list of scantables - merge first
87 merged = None
88 if len(alignedlst) > 1:
89 merged = merge(alignedlst)
90 else:
91 merged = alignedlst[0]
92 s = scantable(stm._averagechannel(merged, 'MEDIAN', scanav))
93 del merged
94 else:
95 #s = scantable(stm._average(alignedlst, mask, weight.upper(), scanav))
96 s = scantable(stm._new_average(alignedlst, compel, mask, weight.upper(), scanav))
97 s._add_history("average_time",varlist)
98 print_log()
99 return s
100
101def quotient(source, reference, preserve=True):
102 """
103 Return the quotient of a 'source' (signal) scan and a 'reference' scan.
104 The reference can have just one scan, even if the signal has many. Otherwise
105 they must have the same number of scans.
106 The cursor of the output scan is set to 0
107 Parameters:
108 source: the 'on' scan
109 reference: the 'off' scan
110 preserve: you can preserve (default) the continuum or
111 remove it. The equations used are
112 preserve: Output = Toff * (on/off) - Toff
113 remove: Output = Toff * (on/off) - Ton
114 """
115 varlist = vars()
116 from asap._asap import stmath
117 stm = stmath()
118 stm._setinsitu(False)
119 s = scantable(stm._quotient(source, reference, preserve))
120 s._add_history("quotient",varlist)
121 print_log()
122 return s
123
124def dototalpower(calon, caloff, tcalval=0.0):
125 """
126 Do calibration for CAL on,off signals.
127 Adopted from GBTIDL dototalpower
128 Parameters:
129 calon: the 'cal on' subintegration
130 caloff: the 'cal off' subintegration
131 tcalval: user supplied Tcal value
132 """
133 varlist = vars()
134 from asap._asap import stmath
135 stm = stmath()
136 stm._setinsitu(False)
137 s = scantable(stm._dototalpower(calon, caloff, tcalval))
138 s._add_history("dototalpower",varlist)
139 print_log()
140 return s
141
142def dosigref(sig, ref, smooth, tsysval=0.0, tauval=0.0):
143 """
144 Calculate a quotient (sig-ref/ref * Tsys)
145 Adopted from GBTIDL dosigref
146 Parameters:
147 sig: on source data
148 ref: reference data
149 smooth: width of box car smoothing for reference
150 tsysval: user specified Tsys (scalar only)
151 tauval: user specified Tau (required if tsysval is set)
152 """
153 varlist = vars()
154 from asap._asap import stmath
155 stm = stmath()
156 stm._setinsitu(False)
157 s = scantable(stm._dosigref(sig, ref, smooth, tsysval, tauval))
158 s._add_history("dosigref",varlist)
159 print_log()
160 return s
161
162def calps(scantab, scannos, smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0, verify=False):
163 """
164 Calibrate GBT position switched data
165 Adopted from GBTIDL getps
166 Currently calps identify the scans as position switched data if they
167 contain '_ps' in the source name. The data must contains 'CAL' signal
168 on/off in each integration. To identify 'CAL' on state, the word, 'calon'
169 need to be present in the source name field.
170 (GBT MS data reading process to scantable automatically append these
171 id names to the source names)
172
173 Parameters:
174 scantab: scantable
175 scannos: list of scan numbers
176 smooth: optional box smoothing order for the reference
177 (default is 1 = no smoothing)
178 tsysval: optional user specified Tsys (default is 0.0,
179 use Tsys in the data)
180 tauval: optional user specified Tau
181 tcalval: optional user specified Tcal (default is 0.0,
182 use Tcal value in the data)
183 """
184 varlist = vars()
185 # check for the appropriate data
186 s = scantab.get_scan('*_ps*')
187 if s is None:
188 msg = "The input data appear to contain no position-switch mode data."
189 if rcParams['verbose']:
190 #print msg
191 asaplog.push(msg)
192 print_log('ERROR')
193 return
194 else:
195 raise TypeError(msg)
196 ssub = s.get_scan(scannos)
197 if ssub is None:
198 msg = "No data was found with given scan numbers!"
199 if rcParams['verbose']:
200 #print msg
201 asaplog.push(msg)
202 print_log('ERROR')
203 return
204 else:
205 raise TypeError(msg)
206 ssubon = ssub.get_scan('*calon')
207 ssuboff = ssub.get_scan('*[^calon]')
208 if ssubon.nrow() != ssuboff.nrow():
209 msg = "mismatch in numbers of CAL on/off scans. Cannot calibrate. Check the scan numbers."
210 if rcParams['verbose']:
211 #print msg
212 asaplog.push(msg)
213 print_log('ERROR')
214 return
215 else:
216 raise TypeError(msg)
217 cals = dototalpower(ssubon, ssuboff, tcalval)
218 sig = cals.get_scan('*ps')
219 ref = cals.get_scan('*psr')
220 if sig.nscan() != ref.nscan():
221 msg = "mismatch in numbers of on/off scans. Cannot calibrate. Check the scan numbers."
222 if rcParams['verbose']:
223 #print msg
224 asaplog.push(msg)
225 print_log('ERROR')
226 return
227 else:
228 raise TypeError(msg)
229
230 #for user supplied Tsys
231 if tsysval>0.0:
232 if tauval<=0.0:
233 msg = "Need to supply a valid tau to use the supplied Tsys"
234 if rcParams['verbose']:
235 #print msg
236 asaplog.push(msg)
237 print_log('ERROR')
238 return
239 else:
240 raise TypeError(msg)
241 else:
242 sig.recalc_azel()
243 ref.recalc_azel()
244 #msg = "Use of user specified Tsys is not fully implemented yet."
245 #if rcParams['verbose']:
246 # print msg
247 # return
248 #else:
249 # raise TypeError(msg)
250 # use get_elevation to get elevation and
251 # calculate a scaling factor using the formula
252 # -> tsys use to dosigref
253
254 #ress = dosigref(sig, ref, smooth, tsysval)
255 ress = dosigref(sig, ref, smooth, tsysval, tauval)
256 ###
257 if verify:
258 # get data
259 import numpy
260 precal={}
261 postcal=[]
262 keys=['ps','ps_calon','psr','psr_calon']
263 ifnos=list(ssub.getifnos())
264 polnos=list(ssub.getpolnos())
265 sel=selector()
266 for i in range(2):
267 ss=ssuboff.get_scan('*'+keys[2*i])
268 ll=[]
269 for j in range(len(ifnos)):
270 for k in range(len(polnos)):
271 sel.set_ifs(ifnos[j])
272 sel.set_polarizations(polnos[k])
273 try:
274 ss.set_selection(sel)
275 except:
276 continue
277 ll.append(numpy.array(ss._getspectrum(0)))
278 sel.reset()
279 ss.set_selection()
280 precal[keys[2*i]]=ll
281 del ss
282 ss=ssubon.get_scan('*'+keys[2*i+1])
283 ll=[]
284 for j in range(len(ifnos)):
285 for k in range(len(polnos)):
286 sel.set_ifs(ifnos[j])
287 sel.set_polarizations(polnos[k])
288 try:
289 ss.set_selection(sel)
290 except:
291 continue
292 ll.append(numpy.array(ss._getspectrum(0)))
293 sel.reset()
294 ss.set_selection()
295 precal[keys[2*i+1]]=ll
296 del ss
297 for j in range(len(ifnos)):
298 for k in range(len(polnos)):
299 sel.set_ifs(ifnos[j])
300 sel.set_polarizations(polnos[k])
301 try:
302 ress.set_selection(sel)
303 except:
304 continue
305 postcal.append(numpy.array(ress._getspectrum(0)))
306 sel.reset()
307 ress.set_selection()
308 del sel
309 # plot
310 print_log()
311 asaplog.push('Plot only first spectrum for each [if,pol] pairs to verify calibration.')
312 print_log('WARN')
313 p=asaplotgui.asaplotgui()
314 #nr=min(6,len(ifnos)*len(polnos))
315 nr=len(ifnos)*len(polnos)
316 titles=[]
317 btics=[]
318 if nr<4:
319 p.set_panels(rows=nr,cols=2,nplots=2*nr,ganged=False)
320 for i in range(2*nr):
321 b=False
322 if i >= 2*nr-2:
323 b=True
324 btics.append(b)
325 elif nr==4:
326 p.set_panels(rows=2,cols=4,nplots=8,ganged=False)
327 for i in range(2*nr):
328 b=False
329 if i >= 2*nr-4:
330 b=True
331 btics.append(b)
332 elif nr<7:
333 p.set_panels(rows=3,cols=4,nplots=2*nr,ganged=False)
334 for i in range(2*nr):
335 if i >= 2*nr-4:
336 b=True
337 btics.append(b)
338 else:
339 print_log()
340 asaplog.push('Only first 6 [if,pol] pairs are plotted.')
341 print_log('WARN')
342 nr=6
343 for i in range(2*nr):
344 b=False
345 if i >= 2*nr-4:
346 b=True
347 btics.append(b)
348 p.set_panels(rows=3,cols=4,nplots=2*nr,ganged=False)
349 for i in range(nr):
350 p.subplot(2*i)
351 p.color=0
352 title='raw data IF%s POL%s' % (ifnos[int(i/len(polnos))],polnos[i%len(polnos)])
353 titles.append(title)
354 #p.set_axes('title',title,fontsize=40)
355 ymin=1.0e100
356 ymax=-1.0e100
357 nchan=s.nchan()
358 edge=int(nchan*0.01)
359 for j in range(4):
360 spmin=min(precal[keys[j]][i][edge:nchan-edge])
361 spmax=max(precal[keys[j]][i][edge:nchan-edge])
362 ymin=min(ymin,spmin)
363 ymax=max(ymax,spmax)
364 for j in range(4):
365 if i==0:
366 p.set_line(label=keys[j])
367 else:
368 p.legend()
369 p.plot(precal[keys[j]][i])
370 p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax))
371 if not btics[2*i]:
372 p.axes.set_xticks([])
373 p.subplot(2*i+1)
374 p.color=0
375 title='cal data IF%s POL%s' % (ifnos[int(i/len(polnos))],polnos[i%len(polnos)])
376 titles.append(title)
377 #p.set_axes('title',title)
378 p.legend()
379 ymin=postcal[i][edge:nchan-edge].min()
380 ymax=postcal[i][edge:nchan-edge].max()
381 p.plot(postcal[i])
382 p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax))
383 if not btics[2*i+1]:
384 p.axes.set_xticks([])
385 for i in range(2*nr):
386 p.subplot(i)
387 p.set_axes('title',titles[i],fontsize='medium')
388 x=raw_input('Accept calibration ([y]/n): ' )
389 if x.upper() == 'N':
390 p.unmap()
391 del p
392 return scabtab
393 p.unmap()
394 del p
395 ###
396 ress._add_history("calps", varlist)
397 print_log()
398 return ress
399
400def calnod(scantab, scannos=[], smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0, verify=False):
401 """
402 Do full (but a pair of scans at time) processing of GBT Nod data
403 calibration.
404 Adopted from GBTIDL's getnod
405 Parameters:
406 scantab: scantable
407 scannos: a pair of scan numbers, or the first scan number of the pair
408 smooth: box car smoothing order
409 tsysval: optional user specified Tsys value
410 tauval: optional user specified tau value (not implemented yet)
411 tcalval: optional user specified Tcal value
412 """
413 varlist = vars()
414 from asap._asap import stmath
415 stm = stmath()
416 stm._setinsitu(False)
417
418 # check for the appropriate data
419 s = scantab.get_scan('*_nod*')
420 if s is None:
421 msg = "The input data appear to contain no Nod observing mode data."
422 if rcParams['verbose']:
423 #print msg
424 asaplog.push(msg)
425 print_log('ERROR')
426 return
427 else:
428 raise TypeError(msg)
429
430 # need check correspondance of each beam with sig-ref ...
431 # check for timestamps, scan numbers, subscan id (not available in
432 # ASAP data format...). Assume 1st scan of the pair (beam 0 - sig
433 # and beam 1 - ref...)
434 # First scan number of paired scans or list of pairs of
435 # scan numbers (has to have even number of pairs.)
436
437 #data splitting
438 scan1no = scan2no = 0
439
440 if len(scannos)==1:
441 scan1no = scannos[0]
442 scan2no = scannos[0]+1
443 pairScans = [scan1no, scan2no]
444 else:
445 #if len(scannos)>2:
446 # msg = "calnod can only process a pair of nod scans at time."
447 # if rcParams['verbose']:
448 # print msg
449 # return
450 # else:
451 # raise TypeError(msg)
452 #
453 #if len(scannos)==2:
454 # scan1no = scannos[0]
455 # scan2no = scannos[1]
456 pairScans = list(scannos)
457
458 if tsysval>0.0:
459 if tauval<=0.0:
460 msg = "Need to supply a valid tau to use the supplied Tsys"
461 if rcParams['verbose']:
462 #print msg
463 asaplog.push(msg)
464 print_log('ERROR')
465 return
466 else:
467 raise TypeError(msg)
468 else:
469 scantab.recalc_azel()
470 resspec = scantable(stm._donod(scantab, pairScans, smooth, tsysval,tauval,tcalval))
471 ###
472 if verify:
473 # get data
474 import numpy
475 precal={}
476 postcal=[]
477 keys=['nod','nod_calon']
478 ifnos=list(scantab.getifnos())
479 polnos=list(scantab.getpolnos())
480 sel=selector()
481 for i in range(2):
482 ss=scantab.get_scan('*'+keys[i])
483 ll=[]
484 ll2=[]
485 for j in range(len(ifnos)):
486 for k in range(len(polnos)):
487 sel.set_ifs(ifnos[j])
488 sel.set_polarizations(polnos[k])
489 sel.set_scans(pairScans[0])
490 try:
491 ss.set_selection(sel)
492 except:
493 continue
494 ll.append(numpy.array(ss._getspectrum(0)))
495 sel.reset()
496 ss.set_selection()
497 sel.set_ifs(ifnos[j])
498 sel.set_polarizations(polnos[k])
499 sel.set_scans(pairScans[1])
500 try:
501 ss.set_selection(sel)
502 except:
503 ll.pop()
504 continue
505 ll2.append(numpy.array(ss._getspectrum(0)))
506 sel.reset()
507 ss.set_selection()
508 key='%s%s' %(pairScans[0],keys[i].lstrip('nod'))
509 precal[key]=ll
510 key='%s%s' %(pairScans[1],keys[i].lstrip('nod'))
511 precal[key]=ll2
512 del ss
513 keys=precal.keys()
514 for j in range(len(ifnos)):
515 for k in range(len(polnos)):
516 sel.set_ifs(ifnos[j])
517 sel.set_polarizations(polnos[k])
518 sel.set_scans(pairScans[0])
519 try:
520 resspec.set_selection(sel)
521 except:
522 continue
523 postcal.append(numpy.array(resspec._getspectrum(0)))
524 sel.reset()
525 resspec.set_selection()
526 del sel
527 # plot
528 print_log()
529 asaplog.push('Plot only first spectrum for each [if,pol] pairs to verify calibration.')
530 print_log('WARN')
531 p=asaplotgui.asaplotgui()
532 #nr=min(6,len(ifnos)*len(polnos))
533 nr=len(ifnos)*len(polnos)
534 titles=[]
535 btics=[]
536 if nr<4:
537 p.set_panels(rows=nr,cols=2,nplots=2*nr,ganged=False)
538 for i in range(2*nr):
539 b=False
540 if i >= 2*nr-2:
541 b=True
542 btics.append(b)
543 elif nr==4:
544 p.set_panels(rows=2,cols=4,nplots=8,ganged=False)
545 for i in range(2*nr):
546 b=False
547 if i >= 2*nr-4:
548 b=True
549 btics.append(b)
550 elif nr<7:
551 p.set_panels(rows=3,cols=4,nplots=2*nr,ganged=False)
552 for i in range(2*nr):
553 if i >= 2*nr-4:
554 b=True
555 btics.append(b)
556 else:
557 print_log()
558 asaplog.push('Only first 6 [if,pol] pairs are plotted.')
559 print_log('WARN')
560 nr=6
561 for i in range(2*nr):
562 b=False
563 if i >= 2*nr-4:
564 b=True
565 btics.append(b)
566 p.set_panels(rows=3,cols=4,nplots=2*nr,ganged=False)
567 for i in range(nr):
568 p.subplot(2*i)
569 p.color=0
570 title='raw data IF%s POL%s' % (ifnos[int(i/len(polnos))],polnos[i%len(polnos)])
571 titles.append(title)
572 #p.set_axes('title',title,fontsize=40)
573 ymin=1.0e100
574 ymax=-1.0e100
575 nchan=scantab.nchan()
576 edge=int(nchan*0.01)
577 for j in range(4):
578 spmin=min(precal[keys[j]][i][edge:nchan-edge])
579 spmax=max(precal[keys[j]][i][edge:nchan-edge])
580 ymin=min(ymin,spmin)
581 ymax=max(ymax,spmax)
582 for j in range(4):
583 if i==0:
584 p.set_line(label=keys[j])
585 else:
586 p.legend()
587 p.plot(precal[keys[j]][i])
588 p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax))
589 if not btics[2*i]:
590 p.axes.set_xticks([])
591 p.subplot(2*i+1)
592 p.color=0
593 title='cal data IF%s POL%s' % (ifnos[int(i/len(polnos))],polnos[i%len(polnos)])
594 titles.append(title)
595 #p.set_axes('title',title)
596 p.legend()
597 ymin=postcal[i][edge:nchan-edge].min()
598 ymax=postcal[i][edge:nchan-edge].max()
599 p.plot(postcal[i])
600 p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax))
601 if not btics[2*i+1]:
602 p.axes.set_xticks([])
603 for i in range(2*nr):
604 p.subplot(i)
605 p.set_axes('title',titles[i],fontsize='medium')
606 x=raw_input('Accept calibration ([y]/n): ' )
607 if x.upper() == 'N':
608 p.unmap()
609 del p
610 return scabtab
611 p.unmap()
612 del p
613 ###
614 resspec._add_history("calnod",varlist)
615 print_log()
616 return resspec
617
618def calfs(scantab, scannos=[], smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0, verify=False):
619 """
620 Calibrate GBT frequency switched data.
621 Adopted from GBTIDL getfs.
622 Currently calfs identify the scans as frequency switched data if they
623 contain '_fs' in the source name. The data must contains 'CAL' signal
624 on/off in each integration. To identify 'CAL' on state, the word, 'calon'
625 need to be present in the source name field.
626 (GBT MS data reading via scantable automatically append these
627 id names to the source names)
628
629 Parameters:
630 scantab: scantable
631 scannos: list of scan numbers
632 smooth: optional box smoothing order for the reference
633 (default is 1 = no smoothing)
634 tsysval: optional user specified Tsys (default is 0.0,
635 use Tsys in the data)
636 tauval: optional user specified Tau
637 """
638 varlist = vars()
639 from asap._asap import stmath
640 stm = stmath()
641 stm._setinsitu(False)
642
643# check = scantab.get_scan('*_fs*')
644# if check is None:
645# msg = "The input data appear to contain no Nod observing mode data."
646# if rcParams['verbose']:
647# print msg
648# return
649# else:
650# raise TypeError(msg)
651 s = scantab.get_scan(scannos)
652 del scantab
653
654 resspec = scantable(stm._dofs(s, scannos, smooth, tsysval,tauval,tcalval))
655 ###
656 if verify:
657 # get data
658 ssub = s.get_scan(scannos)
659 ssubon = ssub.get_scan('*calon')
660 ssuboff = ssub.get_scan('*[^calon]')
661 import numpy
662 precal={}
663 postcal=[]
664 keys=['fs','fs_calon','fsr','fsr_calon']
665 ifnos=list(ssub.getifnos())
666 polnos=list(ssub.getpolnos())
667 sel=selector()
668 for i in range(2):
669 ss=ssuboff.get_scan('*'+keys[2*i])
670 ll=[]
671 for j in range(len(ifnos)):
672 for k in range(len(polnos)):
673 sel.set_ifs(ifnos[j])
674 sel.set_polarizations(polnos[k])
675 try:
676 ss.set_selection(sel)
677 except:
678 continue
679 ll.append(numpy.array(ss._getspectrum(0)))
680 sel.reset()
681 ss.set_selection()
682 precal[keys[2*i]]=ll
683 del ss
684 ss=ssubon.get_scan('*'+keys[2*i+1])
685 ll=[]
686 for j in range(len(ifnos)):
687 for k in range(len(polnos)):
688 sel.set_ifs(ifnos[j])
689 sel.set_polarizations(polnos[k])
690 try:
691 ss.set_selection(sel)
692 except:
693 continue
694 ll.append(numpy.array(ss._getspectrum(0)))
695 sel.reset()
696 ss.set_selection()
697 precal[keys[2*i+1]]=ll
698 del ss
699 sig=resspec.get_scan('*_fs')
700 ref=resspec.get_scan('*_fsr')
701 for k in range(len(polnos)):
702 for j in range(len(ifnos)):
703 sel.set_ifs(ifnos[j])
704 sel.set_polarizations(polnos[k])
705 try:
706 sig.set_selection(sel)
707 postcal.append(numpy.array(sig._getspectrum(0)))
708 except:
709 ref.set_selection(sel)
710 postcal.append(numpy.array(ref._getspectrum(0)))
711 sel.reset()
712 resspec.set_selection()
713 del sel
714 # plot
715 print_log()
716 asaplog.push('Plot only first spectrum for each [if,pol] pairs to verify calibration.')
717 print_log('WARN')
718 p=asaplotgui.asaplotgui()
719 #nr=min(6,len(ifnos)*len(polnos))
720 nr=len(ifnos)/2*len(polnos)
721 titles=[]
722 btics=[]
723 if nr>3:
724 print_log()
725 asaplog.push('Only first 3 [if,pol] pairs are plotted.')
726 print_log('WARN')
727 nr=3
728 p.set_panels(rows=nr,cols=3,nplots=3*nr,ganged=False)
729 for i in range(3*nr):
730 b=False
731 if i >= 3*nr-3:
732 b=True
733 btics.append(b)
734 for i in range(nr):
735 p.subplot(3*i)
736 p.color=0
737 title='raw data IF%s,%s POL%s' % (ifnos[2*int(i/len(polnos))],ifnos[2*int(i/len(polnos))+1],polnos[i%len(polnos)])
738 titles.append(title)
739 #p.set_axes('title',title,fontsize=40)
740 ymin=1.0e100
741 ymax=-1.0e100
742 nchan=s.nchan()
743 edge=int(nchan*0.01)
744 for j in range(4):
745 spmin=min(precal[keys[j]][i][edge:nchan-edge])
746 spmax=max(precal[keys[j]][i][edge:nchan-edge])
747 ymin=min(ymin,spmin)
748 ymax=max(ymax,spmax)
749 for j in range(4):
750 if i==0:
751 p.set_line(label=keys[j])
752 else:
753 p.legend()
754 p.plot(precal[keys[j]][i])
755 p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax))
756 if not btics[3*i]:
757 p.axes.set_xticks([])
758 p.subplot(3*i+1)
759 p.color=0
760 title='sig data IF%s POL%s' % (ifnos[2*int(i/len(polnos))],polnos[i%len(polnos)])
761 titles.append(title)
762 #p.set_axes('title',title)
763 p.legend()
764 ymin=postcal[2*i][edge:nchan-edge].min()
765 ymax=postcal[2*i][edge:nchan-edge].max()
766 p.plot(postcal[2*i])
767 p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax))
768 if not btics[3*i+1]:
769 p.axes.set_xticks([])
770 p.subplot(3*i+2)
771 p.color=0
772 title='ref data IF%s POL%s' % (ifnos[2*int(i/len(polnos))+1],polnos[i%len(polnos)])
773 titles.append(title)
774 #p.set_axes('title',title)
775 p.legend()
776 ymin=postcal[2*i+1][edge:nchan-edge].min()
777 ymax=postcal[2*i+1][edge:nchan-edge].max()
778 p.plot(postcal[2*i+1])
779 p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax))
780 if not btics[3*i+2]:
781 p.axes.set_xticks([])
782 for i in range(3*nr):
783 p.subplot(i)
784 p.set_axes('title',titles[i],fontsize='medium')
785 x=raw_input('Accept calibration ([y]/n): ' )
786 if x.upper() == 'N':
787 p.unmap()
788 del p
789 return scabtab
790 p.unmap()
791 del p
792 ###
793 resspec._add_history("calfs",varlist)
794 print_log()
795 return resspec
796
797def simple_math(left, right, op='add', tsys=True):
798 """
799 Apply simple mathematical binary operations to two
800 scan tables, returning the result in a new scan table.
801 The operation is applied to both the correlations and the TSys data
802 The cursor of the output scan is set to 0
803 Parameters:
804 left: the 'left' scan
805 right: the 'right' scan
806 op: the operation: 'add' (default), 'sub', 'mul', 'div'
807 tsys: if True (default) then apply the operation to Tsys
808 as well as the data
809 """
810 #print "simple_math is deprecated use +=/* instead."
811 asaplog.push( "simple_math is deprecated use +=/* instead." )
812 print_log('WARN')
813
814def merge(*args):
815 """
816 Merge a list of scanatables, or comma-sperated scantables into one
817 scnatble.
818 Parameters:
819 A list [scan1, scan2] or scan1, scan2.
820 Example:
821 myscans = [scan1, scan2]
822 allscans = merge(myscans)
823 # or equivalent
824 sameallscans = merge(scan1, scan2)
825 """
826 varlist = vars()
827 if isinstance(args[0],list):
828 lst = tuple(args[0])
829 elif isinstance(args[0],tuple):
830 lst = args[0]
831 else:
832 lst = tuple(args)
833 varlist["args"] = "%d scantables" % len(lst)
834 # need special formatting her for history...
835 from asap._asap import stmath
836 stm = stmath()
837 for s in lst:
838 if not isinstance(s,scantable):
839 msg = "Please give a list of scantables"
840 if rcParams['verbose']:
841 #print msg
842 asaplog.push(msg)
843 print_log('ERROR')
844 return
845 else:
846 raise TypeError(msg)
847 s = scantable(stm._merge(lst))
848 s._add_history("merge", varlist)
849 print_log()
850 return s
851
852## def apexcal( scantab, calmode ):
853## """
854## Calibrate APEX data.
855
856## Parameters:
857## scantab: scantable
858## calmode: calibration mode
859## """
860## from asap._asap import stmath
861## stm = stmath()
862## if ( calmode == 'ps' ):
863## asaplog.push( 'APEX position-switch calibration' )
864## print_log()
865## elif ( calmode == 'fs' ):
866## asaplog.push( 'APEX frequency-switch calibration' )
867## print_log()
868## elif ( calmode == 'wob' ):
869## asaplog.push( 'APEX wobbler-switch calibration' )
870## print_log()
871## elif ( calmode == 'otf' ):
872## asaplog.push( 'APEX On-The-Fly calibration' )
873## print_log()
874## else:
875## asaplog.push( '%s: unknown calibration mode' % calmode )
876## print_log('ERROR')
Note: See TracBrowser for help on using the repository browser.