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.