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

Last change on this file since 1651 was 1651, checked in by Kana Sugimoto, 15 years ago

New Development: No

JIRA Issue: Yes (CAS-1441)

Ready to Release: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs:

Put in Release Notes: No

Module(s): Module Names change impacts.

Description: A bit improvement in a parameter handling.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 35.5 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
852def calibrate( scantab, scannos=[], calmode='none', verify=None ):
853    """
854    Calibrate data.
855   
856    Parameters:
857        scantab:       scantable
858        scannos:       list of scan number
859        calmode:       calibration mode
860        verify:        verify calibration     
861    """
862    antname = scantab.get_antennaname()
863    if ( calmode == 'nod' ):
864        asaplog.push( 'Calibrating nod data.' )
865        print_log()
866        scal = calnod( scantab, scannos=scannos, verify=verify )
867    elif ( calmode == 'quotient' ):
868        asaplog.push( 'Calibrating using quotient.' )
869        print_log()
870        scal = scantab.auto_quotient( verify=verify )
871    elif ( calmode == 'ps' ):
872        asaplog.push( 'Calibrating %s position-switched data.' % antname )
873        print_log()
874        if ( antname.find( 'APEX' ) != -1 or antname.find( 'ALMA' ) != -1 ):
875            scal = almacal( scantab, scannos, calmode, verify )
876        else:
877            scal = calps( scantab, scannos=scannos, verify=verify )
878    elif ( calmode == 'fs' or calmode == 'fsotf' ):
879        asaplog.push( 'Calibrating %s frequency-switched data.' % antname )
880        print_log()
881        if ( antname.find( 'APEX' ) != -1 or antname.find( 'ALMA' ) != -1 ):
882            scal = almacal( scantab, scannos, calmode, verify )
883        else:
884            scal = calfs( scantab, scannos=scannos, verify=verify )
885    elif ( calmode == 'otf' ):
886        asaplog.push( 'Calibrating %s On-The-Fly data.' % antname )
887        print_log()
888        scal = almacal( scantab, scannos, calmode, verify )
889    else:
890        asaplog.push( 'No calibration.' )
891        scal = scantab.copy()
892
893    return scal
894
895def almacal( scantab, scannos=[], calmode='none', verify=False ):
896    """
897    Calibrate APEX and ALMA data
898
899    Parameters:
900        scantab:       scantable
901        scannos:       list of scan number
902        calmode:       calibration mode
903
904        verify:        verify calibration     
905    """
906    from asap._asap import stmath
907    stm = stmath()
908    antname = scantab.get_antennaname()
909    ssub = scantab.get_scan( scannos )
910    scal = scantable( stm.cwcal( ssub, calmode, antname ) )
911    return scal
912
913### Start mod: 2009.10.16 kana ###
914def splitant(filename, outprefix='',overwrite=False):
915    """
916    Split Measurement set by antenna name, save data as a scantables,
917    and return a list of filename.
918    Notice this method can only be available from CASA.
919    Prameter
920       filename:    the name of Measurement set to be read.
921       outprefix:   the prefix of output scantable name.
922                    the names of output scantable will be
923                    outprefix.antenna1, outprefix.antenna2, ....
924                    If not specified, outprefix = filename is assumed.
925       overwrite    If the file should be overwritten if it exists.
926                    The default False is to return with warning
927                    without writing the output. USE WITH CARE.
928                   
929    """
930    # Import the table toolkit from CASA
931    try:
932        import casac
933    except ImportError:
934        if rcParams['verbose']:
935            #print "failed to load casa"
936            print_log()
937            asaplog.push("failed to load casa")
938            print_log('ERROR')
939        else: raise
940        return False
941    try:
942        tbtool = casac.homefinder.find_home_by_name('tableHome')
943        tb = tbtool.create()
944        tb2 = tbtool.create()
945    except:
946        if rcParams['verbose']:
947            #print "failed to load a table tool:\n", e
948            print_log()
949            asaplog.push("failed to load table tool")
950            print_log('ERROR')
951        else: raise
952        return False
953    # Check the input filename
954    if isinstance(filename, str):
955        import os.path
956        filename = os.path.expandvars(filename)
957        filename = os.path.expanduser(filename)
958        if not os.path.exists(filename):
959            s = "File '%s' not found." % (filename)
960            if rcParams['verbose']:
961                print_log()
962                asaplog.push(s)
963                print_log('ERROR')
964                return
965            raise IOError(s)
966        # check if input file is MS
967        if not os.path.isdir(filename) \
968               or not os.path.exists(filename+'/ANTENNA') \
969               or not os.path.exists(filename+'/table.f1'):
970            s = "File '%s' is not a Measurement set." % (filename)
971            if rcParams['verbose']:
972                print_log()
973                asaplog.push(s)
974                print_log('ERROR')
975                return
976            raise IOError(s)
977    else:
978        s = "The filename should be string. "
979        if rcParams['verbose']:
980            print_log()
981            asaplog.push(s)
982            print_log('ERROR')
983            return
984        raise TypeError(s)
985    # Check out put file name
986    outname=''
987    if len(outprefix) > 0: prefix=outprefix+'.'
988    else:
989        prefix=filename
990    # Now do the actual splitting.
991    outfiles=[]
992    tmpms="temp_antsplit.ms"
993    if os.path.exists(tmpms):
994        ans=raw_input('Temporal file '+tmpms+' exists. Delete it and continue? [y/N]: ')
995        if ans.upper() == 'Y':
996            os.system('rm -rf '+tmpms)
997            asaplog.push('The file '+tmpms+' deleted.')
998        else:
999            asaplog.push('Exit without splitting.')
1000            return
1001    tb.open(tablename=filename+'/ANTENNA',nomodify=True)
1002    nant=tb.nrows()
1003    antnames=tb.getcol('NAME',0,nant,1)
1004    antpos=tb.getcol('POSITION',0,nant,1).transpose()
1005    tb.close()
1006    tb.open(tablename=filename,nomodify=True)
1007    ant1=tb.getcol('ANTENNA1',0,-1,1)
1008    for antid in set(ant1):
1009        qstr="ANTENNA1 == "+str(antid)
1010        stab = tb.queryC(qstr)
1011        ctab = stab.copy(tmpms,deep=True)
1012        stab.close()
1013        ctab.close()
1014        scan=scantable(tmpms,False)
1015        outname=prefix+antnames[antid]+'.asap'
1016        scan.save(outname,format='ASAP',overwrite=overwrite)
1017        # Modify scantable header
1018        tb2.open(tablename=outname,nomodify=False)
1019        tb2.putkeyword(keyword='AntennaName',value=antnames[antid])
1020        tb2.putkeyword(keyword='AntennaPosition',value=antpos[antid])
1021        tb2.flush()
1022        tb2.close()
1023        del scan, ctab, stab
1024        outfiles.append(outname)
1025    tb.close()
1026    del tb, tb2
1027    os.system('rm -rf '+tmpms)
1028    return outfiles
1029### End mod ######################
Note: See TracBrowser for help on using the repository browser.