Ignore:
Timestamp:
07/29/10 19:13:46 (14 years ago)
Author:
Kana Sugimoto
Message:

New Development: Yes

JIRA Issue: No (test merging alma branch)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s):

Description:


Location:
branches/mergetest
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • branches/mergetest

  • branches/mergetest/python/asapmath.py

    r1591 r1779  
    11from asap.scantable import scantable
    22from asap import rcParams
    3 from asap import print_log_dec
     3from asap import print_log, print_log_dec
    44from asap import selector
    5 
    6 @print_log_dec
     5from asap import asaplog
     6from asap import asaplotgui
     7
     8#@print_log_dec
    79def average_time(*args, **kwargs):
    810    """
    … …  
    4648    if kwargs.has_key('align'):
    4749        align = kwargs.get('align')
     50    compel = False
     51    if kwargs.has_key('compel'):
     52        compel = kwargs.get('compel')
    4853    varlist = vars()
    4954    if isinstance(args[0],list):
    … …  
    6469            msg = "Please give a list of scantables"
    6570            if rcParams['verbose']:
    66                 print msg
     71                #print msg
     72                asaplog.push(msg)
     73                print_log('ERROR')
    6774                return
    6875            else:
    … …  
    8794        del merged
    8895    else:
    89         s = scantable(stm._average(alignedlst, mask, weight.upper(), scanav))
     96        #s = scantable(stm._average(alignedlst, mask, weight.upper(), scanav))
     97        s = scantable(stm._new_average(alignedlst, compel, mask, weight.upper(), scanav))
    9098    s._add_history("average_time",varlist)
     99    print_log()
    91100    return s
    92101
    … …  
    111120    s = scantable(stm._quotient(source, reference, preserve))
    112121    s._add_history("quotient",varlist)
     122    print_log()
    113123    return s
    114124
    115 @print_log_dec
     125#@print_log_dec
    116126def dototalpower(calon, caloff, tcalval=0.0):
    117127    """
    … …  
    129139    s = scantable(stm._dototalpower(calon, caloff, tcalval))
    130140    s._add_history("dototalpower",varlist)
     141    print_log()
    131142    return s
    132143
    133 @print_log_dec
     144#@print_log_dec
    134145def dosigref(sig, ref, smooth, tsysval=0.0, tauval=0.0):
    135146    """
    … …  
    149160    s = scantable(stm._dosigref(sig, ref, smooth, tsysval, tauval))
    150161    s._add_history("dosigref",varlist)
     162    print_log()
    151163    return s
    152164
    153 @print_log_dec
    154 def calps(scantab, scannos, smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0):
     165#@print_log_dec
     166def calps(scantab, scannos, smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0, verify=False):
    155167    """
    156168    Calibrate GBT position switched data
    … …  
    176188    varlist = vars()
    177189    # check for the appropriate data
    178     s = scantab.get_scan('*_ps*')
    179     if s is None:
     190##    s = scantab.get_scan('*_ps*')
     191##     if s is None:
     192##         msg = "The input data appear to contain no position-switch mode data."
     193##         if rcParams['verbose']:
     194##             #print msg
     195##             asaplog.push(msg)
     196##             print_log('ERROR')
     197##             return
     198##         else:
     199##             raise TypeError(msg)
     200    s = scantab.copy()
     201    from asap._asap import srctype
     202    sel = selector()
     203    sel.set_types( srctype.pson )
     204    try:
     205        scantab.set_selection( sel )
     206    except Exception, e:
    180207        msg = "The input data appear to contain no position-switch mode data."
    181208        if rcParams['verbose']:
    182             print msg
     209            #print msg
     210            asaplog.push(msg)
     211            print_log('ERROR')
    183212            return
    184213        else:
    185214            raise TypeError(msg)
     215    s.set_selection()
     216    sel.reset()
    186217    ssub = s.get_scan(scannos)
    187218    if ssub is None:
    188219        msg = "No data was found with given scan numbers!"
    189220        if rcParams['verbose']:
    190             print msg
     221            #print msg
     222            asaplog.push(msg)
     223            print_log('ERROR')
    191224            return
    192225        else:
    193226            raise TypeError(msg)
    194     ssubon = ssub.get_scan('*calon')
    195     ssuboff = ssub.get_scan('*[^calon]')
     227    #ssubon = ssub.get_scan('*calon')
     228    #ssuboff = ssub.get_scan('*[^calon]')
     229    sel.set_types( [srctype.poncal,srctype.poffcal] )
     230    ssub.set_selection( sel )
     231    ssubon = ssub.copy()
     232    ssub.set_selection()
     233    sel.reset()
     234    sel.set_types( [srctype.pson,srctype.psoff] )
     235    ssub.set_selection( sel )
     236    ssuboff = ssub.copy()
     237    ssub.set_selection()
     238    sel.reset()
    196239    if ssubon.nrow() != ssuboff.nrow():
    197240        msg = "mismatch in numbers of CAL on/off scans. Cannot calibrate. Check the scan numbers."
    198241        if rcParams['verbose']:
    199             print msg
     242            #print msg
     243            asaplog.push(msg)
     244            print_log('ERROR')
    200245            return
    201246        else:
    202247            raise TypeError(msg)
    203248    cals = dototalpower(ssubon, ssuboff, tcalval)
    204     sig = cals.get_scan('*ps')
    205     ref = cals.get_scan('*psr')
     249    #sig = cals.get_scan('*ps')
     250    #ref = cals.get_scan('*psr')
     251    sel.set_types( srctype.pson )
     252    cals.set_selection( sel )
     253    sig = cals.copy()
     254    cals.set_selection()
     255    sel.reset()
     256    sel.set_types( srctype.psoff )
     257    cals.set_selection( sel )
     258    ref = cals.copy()
     259    cals.set_selection()
     260    sel.reset()
    206261    if sig.nscan() != ref.nscan():
    207262        msg = "mismatch in numbers of on/off scans. Cannot calibrate. Check the scan numbers."
    208263        if rcParams['verbose']:
    209             print msg
     264            #print msg
     265            asaplog.push(msg)
     266            print_log('ERROR')
    210267            return
    211268        else:
    … …  
    217274            msg = "Need to supply a valid tau to use the supplied Tsys"
    218275            if rcParams['verbose']:
    219                 print msg
     276                #print msg
     277                asaplog.push(msg)
     278                print_log('ERROR')
    220279                return
    221280            else:
    … …  
    236295    #ress = dosigref(sig, ref, smooth, tsysval)
    237296    ress = dosigref(sig, ref, smooth, tsysval, tauval)
     297    ###
     298    if verify:
     299        # get data
     300        import numpy
     301        precal={}
     302        postcal=[]
     303        keys=['ps','ps_calon','psr','psr_calon']
     304        types=[srctype.pson,srctype.poncal,srctype.psoff,srctype.poffcal]
     305        ifnos=list(ssub.getifnos())
     306        polnos=list(ssub.getpolnos())
     307        sel=selector()
     308        for i in range(2):
     309            #ss=ssuboff.get_scan('*'+keys[2*i])
     310            ll=[]
     311            for j in range(len(ifnos)):
     312                for k in range(len(polnos)):
     313                    sel.set_ifs(ifnos[j])
     314                    sel.set_polarizations(polnos[k])
     315                    sel.set_types(types[2*i])
     316                    try:
     317                        #ss.set_selection(sel)
     318                        ssuboff.set_selection(sel)
     319                    except:
     320                        continue
     321                    #ll.append(numpy.array(ss._getspectrum(0)))
     322                    ll.append(numpy.array(ssuboff._getspectrum(0)))
     323                    sel.reset()
     324                    ssuboff.set_selection()
     325            precal[keys[2*i]]=ll
     326            #del ss
     327            #ss=ssubon.get_scan('*'+keys[2*i+1])
     328            ll=[]
     329            for j in range(len(ifnos)):
     330                for k in range(len(polnos)):
     331                    sel.set_ifs(ifnos[j])
     332                    sel.set_polarizations(polnos[k])
     333                    sel.set_types(types[2*i+1])
     334                    try:
     335                        #ss.set_selection(sel)
     336                        ssubon.set_selection(sel)
     337                    except:
     338                        continue
     339                    #ll.append(numpy.array(ss._getspectrum(0)))
     340                    ll.append(numpy.array(ssubon._getspectrum(0)))
     341                    sel.reset()
     342                    ssubon.set_selection()
     343            precal[keys[2*i+1]]=ll
     344            #del ss
     345        for j in range(len(ifnos)):
     346            for k in range(len(polnos)):
     347                sel.set_ifs(ifnos[j])
     348                sel.set_polarizations(polnos[k])
     349                try:
     350                    ress.set_selection(sel)
     351                except:
     352                    continue
     353                postcal.append(numpy.array(ress._getspectrum(0)))
     354                sel.reset()
     355                ress.set_selection()
     356        del sel
     357        # plot
     358        print_log()
     359        asaplog.push('Plot only first spectrum for each [if,pol] pairs to verify calibration.')
     360        print_log('WARN')
     361        p=asaplotgui.asaplotgui()
     362        #nr=min(6,len(ifnos)*len(polnos))
     363        nr=len(ifnos)*len(polnos)
     364        titles=[]
     365        btics=[]
     366        if nr<4:
     367            p.set_panels(rows=nr,cols=2,nplots=2*nr,ganged=False)
     368            for i in range(2*nr):
     369                b=False
     370                if i >= 2*nr-2:
     371                    b=True
     372                btics.append(b)
     373        elif nr==4:
     374            p.set_panels(rows=2,cols=4,nplots=8,ganged=False)
     375            for i in range(2*nr):
     376                b=False
     377                if i >= 2*nr-4:
     378                    b=True
     379                btics.append(b)
     380        elif nr<7:
     381            p.set_panels(rows=3,cols=4,nplots=2*nr,ganged=False)
     382            for i in range(2*nr):
     383                if i >= 2*nr-4:
     384                    b=True
     385                btics.append(b)
     386        else:
     387            print_log()
     388            asaplog.push('Only first 6 [if,pol] pairs are plotted.')
     389            print_log('WARN')
     390            nr=6
     391            for i in range(2*nr):
     392                b=False
     393                if i >= 2*nr-4:
     394                    b=True
     395                btics.append(b)
     396            p.set_panels(rows=3,cols=4,nplots=2*nr,ganged=False)
     397        for i in range(nr):
     398            p.subplot(2*i)
     399            p.color=0
     400            title='raw data IF%s POL%s' % (ifnos[int(i/len(polnos))],polnos[i%len(polnos)])
     401            titles.append(title)
     402            #p.set_axes('title',title,fontsize=40)
     403            ymin=1.0e100
     404            ymax=-1.0e100
     405            nchan=s.nchan()
     406            edge=int(nchan*0.01)
     407            for j in range(4):
     408                spmin=min(precal[keys[j]][i][edge:nchan-edge])
     409                spmax=max(precal[keys[j]][i][edge:nchan-edge])
     410                ymin=min(ymin,spmin)
     411                ymax=max(ymax,spmax)
     412            for j in range(4):
     413                if i==0:
     414                    p.set_line(label=keys[j])
     415                else:
     416                    p.legend()
     417                p.plot(precal[keys[j]][i])
     418            p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax))
     419            if not btics[2*i]:
     420                p.axes.set_xticks([])
     421            p.subplot(2*i+1)
     422            p.color=0
     423            title='cal data IF%s POL%s' % (ifnos[int(i/len(polnos))],polnos[i%len(polnos)])
     424            titles.append(title)
     425            #p.set_axes('title',title)
     426            p.legend()
     427            ymin=postcal[i][edge:nchan-edge].min()
     428            ymax=postcal[i][edge:nchan-edge].max()
     429            p.plot(postcal[i])
     430            p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax))
     431            if not btics[2*i+1]:
     432                p.axes.set_xticks([])
     433        for i in range(2*nr):
     434            p.subplot(i)
     435            p.set_axes('title',titles[i],fontsize='medium')
     436        x=raw_input('Accept calibration ([y]/n): ' )
     437        if x.upper() == 'N':
     438            p.unmap()
     439            del p
     440            return scabtab
     441        p.unmap()
     442        del p
     443    ###
    238444    ress._add_history("calps", varlist)
     445    print_log()
    239446    return ress
    240447
    241 @print_log_dec
    242 def calnod(scantab, scannos=[], smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0):
     448#@print_log_dec
     449def calnod(scantab, scannos=[], smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0, verify=False):
    243450    """
    244451    Do full (but a pair of scans at time) processing of GBT Nod data
    … …  
    255462    varlist = vars()
    256463    from asap._asap import stmath
     464    from asap._asap import srctype
    257465    stm = stmath()
    258466    stm._setinsitu(False)
    259467
    260468    # check for the appropriate data
    261     s = scantab.get_scan('*_nod*')
    262     if s is None:
     469##     s = scantab.get_scan('*_nod*')
     470##     if s is None:
     471##         msg = "The input data appear to contain no Nod observing mode data."
     472##         if rcParams['verbose']:
     473##             #print msg
     474##             asaplog.push(msg)
     475##             print_log('ERROR')
     476##             return
     477##         else:
     478##             raise TypeError(msg)
     479    s = scantab.copy()
     480    sel = selector()
     481    sel.set_types( srctype.nod )
     482    try:
     483        s.set_selection( sel )
     484    except Exception, e:
    263485        msg = "The input data appear to contain no Nod observing mode data."
    264486        if rcParams['verbose']:
    265             print msg
     487            #print msg
     488            asaplog.push(msg)
     489            print_log('ERROR')
    266490            return
    267491        else:
    268492            raise TypeError(msg)
     493    sel.reset()
     494    del sel
     495    del s
    269496
    270497    # need check correspondance of each beam with sig-ref ...
    … …  
    300527            msg = "Need to supply a valid tau to use the supplied Tsys"
    301528            if rcParams['verbose']:
    302                 print msg
     529                #print msg
     530                asaplog.push(msg)
     531                print_log('ERROR')
    303532                return
    304533            else:
    … …  
    307536            scantab.recalc_azel()
    308537    resspec = scantable(stm._donod(scantab, pairScans, smooth, tsysval,tauval,tcalval))
     538    ###
     539    if verify:
     540        # get data
     541        import numpy
     542        precal={}
     543        postcal=[]
     544        keys=['','_calon']
     545        types=[srctype.nod,srctype.nodcal]
     546        ifnos=list(scantab.getifnos())
     547        polnos=list(scantab.getpolnos())
     548        sel=selector()
     549        ss = scantab.copy()
     550        for i in range(2):
     551            #ss=scantab.get_scan('*'+keys[i])
     552            ll=[]
     553            ll2=[]
     554            for j in range(len(ifnos)):
     555                for k in range(len(polnos)):
     556                    sel.set_ifs(ifnos[j])
     557                    sel.set_polarizations(polnos[k])
     558                    sel.set_scans(pairScans[0])
     559                    sel.set_types(types[i])
     560                    try:
     561                        ss.set_selection(sel)
     562                    except:
     563                        continue
     564                    ll.append(numpy.array(ss._getspectrum(0)))
     565                    sel.reset()
     566                    ss.set_selection()
     567                    sel.set_ifs(ifnos[j])
     568                    sel.set_polarizations(polnos[k])
     569                    sel.set_scans(pairScans[1])
     570                    sel.set_types(types[i])
     571                    try:
     572                        ss.set_selection(sel)
     573                    except:
     574                        ll.pop()
     575                        continue
     576                    ll2.append(numpy.array(ss._getspectrum(0)))
     577                    sel.reset()
     578                    ss.set_selection()
     579            key='%s%s' %(pairScans[0],keys[i])
     580            precal[key]=ll
     581            key='%s%s' %(pairScans[1],keys[i])
     582            precal[key]=ll2
     583            #del ss
     584        keys=precal.keys()
     585        for j in range(len(ifnos)):
     586            for k in range(len(polnos)):
     587                sel.set_ifs(ifnos[j])
     588                sel.set_polarizations(polnos[k])
     589                sel.set_scans(pairScans[0])
     590                try:
     591                    resspec.set_selection(sel)
     592                except:
     593                    continue
     594                postcal.append(numpy.array(resspec._getspectrum(0)))
     595                sel.reset()
     596                resspec.set_selection()
     597        del sel
     598        # plot
     599        print_log()
     600        asaplog.push('Plot only first spectrum for each [if,pol] pairs to verify calibration.')
     601        print_log('WARN')
     602        p=asaplotgui.asaplotgui()
     603        #nr=min(6,len(ifnos)*len(polnos))
     604        nr=len(ifnos)*len(polnos)
     605        titles=[]
     606        btics=[]
     607        if nr<4:
     608            p.set_panels(rows=nr,cols=2,nplots=2*nr,ganged=False)
     609            for i in range(2*nr):
     610                b=False
     611                if i >= 2*nr-2:
     612                    b=True
     613                btics.append(b)
     614        elif nr==4:
     615            p.set_panels(rows=2,cols=4,nplots=8,ganged=False)
     616            for i in range(2*nr):
     617                b=False
     618                if i >= 2*nr-4:
     619                    b=True
     620                btics.append(b)
     621        elif nr<7:
     622            p.set_panels(rows=3,cols=4,nplots=2*nr,ganged=False)
     623            for i in range(2*nr):
     624                if i >= 2*nr-4:
     625                    b=True
     626                btics.append(b)
     627        else:
     628            print_log()
     629            asaplog.push('Only first 6 [if,pol] pairs are plotted.')
     630            print_log('WARN')
     631            nr=6
     632            for i in range(2*nr):
     633                b=False
     634                if i >= 2*nr-4:
     635                    b=True
     636                btics.append(b)
     637            p.set_panels(rows=3,cols=4,nplots=2*nr,ganged=False)
     638        for i in range(nr):
     639            p.subplot(2*i)
     640            p.color=0
     641            title='raw data IF%s POL%s' % (ifnos[int(i/len(polnos))],polnos[i%len(polnos)])
     642            titles.append(title)
     643            #p.set_axes('title',title,fontsize=40)
     644            ymin=1.0e100
     645            ymax=-1.0e100
     646            nchan=scantab.nchan()
     647            edge=int(nchan*0.01)
     648            for j in range(4):
     649                spmin=min(precal[keys[j]][i][edge:nchan-edge])
     650                spmax=max(precal[keys[j]][i][edge:nchan-edge])
     651                ymin=min(ymin,spmin)
     652                ymax=max(ymax,spmax)
     653            for j in range(4):
     654                if i==0:
     655                    p.set_line(label=keys[j])
     656                else:
     657                    p.legend()
     658                p.plot(precal[keys[j]][i])
     659            p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax))
     660            if not btics[2*i]:
     661                p.axes.set_xticks([])
     662            p.subplot(2*i+1)
     663            p.color=0
     664            title='cal data IF%s POL%s' % (ifnos[int(i/len(polnos))],polnos[i%len(polnos)])
     665            titles.append(title)
     666            #p.set_axes('title',title)
     667            p.legend()
     668            ymin=postcal[i][edge:nchan-edge].min()
     669            ymax=postcal[i][edge:nchan-edge].max()
     670            p.plot(postcal[i])
     671            p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax))
     672            if not btics[2*i+1]:
     673                p.axes.set_xticks([])
     674        for i in range(2*nr):
     675            p.subplot(i)
     676            p.set_axes('title',titles[i],fontsize='medium')
     677        x=raw_input('Accept calibration ([y]/n): ' )
     678        if x.upper() == 'N':
     679            p.unmap()
     680            del p
     681            return scabtab
     682        p.unmap()
     683        del p
     684    ###
    309685    resspec._add_history("calnod",varlist)
     686    print_log()
    310687    return resspec
    311688
    312 @print_log_dec
    313 def calfs(scantab, scannos=[], smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0):
     689#@print_log_dec
     690def calfs(scantab, scannos=[], smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0, verify=False):
    314691    """
    315692    Calibrate GBT frequency switched data.
    … …  
    333710    varlist = vars()
    334711    from asap._asap import stmath
     712    from asap._asap import srctype
    335713    stm = stmath()
    336714    stm._setinsitu(False)
    … …  
    348726
    349727    resspec = scantable(stm._dofs(s, scannos, smooth, tsysval,tauval,tcalval))
     728    ###
     729    if verify:
     730        # get data
     731        ssub = s.get_scan(scannos)
     732        #ssubon = ssub.get_scan('*calon')
     733        #ssuboff = ssub.get_scan('*[^calon]')
     734        sel = selector()
     735        sel.set_types( [srctype.foncal,srctype.foffcal] )
     736        ssub.set_selection( sel )
     737        ssubon = ssub.copy()
     738        ssub.set_selection()
     739        sel.reset()
     740        sel.set_types( [srctype.fson,srctype.fsoff] )
     741        ssub.set_selection( sel )
     742        ssuboff = ssub.copy()
     743        ssub.set_selection()
     744        sel.reset()
     745        import numpy
     746        precal={}
     747        postcal=[]
     748        keys=['fs','fs_calon','fsr','fsr_calon']
     749        types=[srctype.fson,srctype.foncal,srctype.fsoff,srctype.foffcal]
     750        ifnos=list(ssub.getifnos())
     751        polnos=list(ssub.getpolnos())
     752        for i in range(2):
     753            #ss=ssuboff.get_scan('*'+keys[2*i])
     754            ll=[]
     755            for j in range(len(ifnos)):
     756                for k in range(len(polnos)):
     757                    sel.set_ifs(ifnos[j])
     758                    sel.set_polarizations(polnos[k])
     759                    sel.set_types(types[2*i])
     760                    try:
     761                        #ss.set_selection(sel)
     762                        ssuboff.set_selection(sel)
     763                    except:
     764                        continue
     765                    ll.append(numpy.array(ss._getspectrum(0)))
     766                    sel.reset()
     767                    #ss.set_selection()
     768                    ssuboff.set_selection()
     769            precal[keys[2*i]]=ll
     770            #del ss
     771            #ss=ssubon.get_scan('*'+keys[2*i+1])
     772            ll=[]
     773            for j in range(len(ifnos)):
     774                for k in range(len(polnos)):
     775                    sel.set_ifs(ifnos[j])
     776                    sel.set_polarizations(polnos[k])
     777                    sel.set_types(types[2*i+1])
     778                    try:
     779                        #ss.set_selection(sel)
     780                        ssubon.set_selection(sel)
     781                    except:
     782                        continue
     783                    ll.append(numpy.array(ss._getspectrum(0)))
     784                    sel.reset()
     785                    #ss.set_selection()
     786                    ssubon.set_selection()
     787            precal[keys[2*i+1]]=ll
     788            #del ss
     789        #sig=resspec.get_scan('*_fs')
     790        #ref=resspec.get_scan('*_fsr')
     791        sel.set_types( srctype.fson )
     792        resspec.set_selection( sel )
     793        sig=resspec.copy()
     794        resspec.set_selection()
     795        sel.reset()
     796        sel.set_type( srctype.fsoff )
     797        resspec.set_selection( sel )
     798        ref=resspec.copy()
     799        resspec.set_selection()
     800        sel.reset()
     801        for k in range(len(polnos)):
     802            for j in range(len(ifnos)):
     803                sel.set_ifs(ifnos[j])
     804                sel.set_polarizations(polnos[k])
     805                try:
     806                    sig.set_selection(sel)
     807                    postcal.append(numpy.array(sig._getspectrum(0)))
     808                except:
     809                    ref.set_selection(sel)
     810                    postcal.append(numpy.array(ref._getspectrum(0)))
     811                sel.reset()
     812                resspec.set_selection()
     813        del sel
     814        # plot
     815        print_log()
     816        asaplog.push('Plot only first spectrum for each [if,pol] pairs to verify calibration.')
     817        print_log('WARN')
     818        p=asaplotgui.asaplotgui()
     819        #nr=min(6,len(ifnos)*len(polnos))
     820        nr=len(ifnos)/2*len(polnos)
     821        titles=[]
     822        btics=[]
     823        if nr>3:
     824            print_log()
     825            asaplog.push('Only first 3 [if,pol] pairs are plotted.')
     826            print_log('WARN')
     827            nr=3
     828        p.set_panels(rows=nr,cols=3,nplots=3*nr,ganged=False)
     829        for i in range(3*nr):
     830            b=False
     831            if i >= 3*nr-3:
     832                b=True
     833            btics.append(b)
     834        for i in range(nr):
     835            p.subplot(3*i)
     836            p.color=0
     837            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)])
     838            titles.append(title)
     839            #p.set_axes('title',title,fontsize=40)
     840            ymin=1.0e100
     841            ymax=-1.0e100
     842            nchan=s.nchan()
     843            edge=int(nchan*0.01)
     844            for j in range(4):
     845                spmin=min(precal[keys[j]][i][edge:nchan-edge])
     846                spmax=max(precal[keys[j]][i][edge:nchan-edge])
     847                ymin=min(ymin,spmin)
     848                ymax=max(ymax,spmax)
     849            for j in range(4):
     850                if i==0:
     851                    p.set_line(label=keys[j])
     852                else:
     853                    p.legend()
     854                p.plot(precal[keys[j]][i])
     855            p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax))
     856            if not btics[3*i]:
     857                p.axes.set_xticks([])
     858            p.subplot(3*i+1)
     859            p.color=0
     860            title='sig data IF%s POL%s' % (ifnos[2*int(i/len(polnos))],polnos[i%len(polnos)])
     861            titles.append(title)
     862            #p.set_axes('title',title)
     863            p.legend()
     864            ymin=postcal[2*i][edge:nchan-edge].min()
     865            ymax=postcal[2*i][edge:nchan-edge].max()
     866            p.plot(postcal[2*i])
     867            p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax))
     868            if not btics[3*i+1]:
     869                p.axes.set_xticks([])
     870            p.subplot(3*i+2)
     871            p.color=0
     872            title='ref data IF%s POL%s' % (ifnos[2*int(i/len(polnos))+1],polnos[i%len(polnos)])
     873            titles.append(title)
     874            #p.set_axes('title',title)
     875            p.legend()
     876            ymin=postcal[2*i+1][edge:nchan-edge].min()
     877            ymax=postcal[2*i+1][edge:nchan-edge].max()
     878            p.plot(postcal[2*i+1])
     879            p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax))
     880            if not btics[3*i+2]:
     881                p.axes.set_xticks([])
     882        for i in range(3*nr):
     883            p.subplot(i)
     884            p.set_axes('title',titles[i],fontsize='medium')
     885        x=raw_input('Accept calibration ([y]/n): ' )
     886        if x.upper() == 'N':
     887            p.unmap()
     888            del p
     889            return scabtab
     890        p.unmap()
     891        del p
     892    ###
    350893    resspec._add_history("calfs",varlist)
     894    print_log()
    351895    return resspec
    352896
    353 
    354 @print_log_dec
     897#@print_log_dec
    355898def merge(*args):
    356899    """
    … …  
    380923            msg = "Please give a list of scantables"
    381924            if rcParams['verbose']:
    382                 print msg
     925                #print msg
     926                asaplog.push(msg)
     927                print_log('ERROR')
    383928                return
    384929            else:
    … …  
    386931    s = scantable(stm._merge(lst))
    387932    s._add_history("merge", varlist)
     933    print_log()
    388934    return s
     935
     936def calibrate( scantab, scannos=[], calmode='none', verify=None ):
     937    """
     938    Calibrate data.
     939   
     940    Parameters:
     941        scantab:       scantable
     942        scannos:       list of scan number
     943        calmode:       calibration mode
     944        verify:        verify calibration     
     945    """
     946    antname = scantab.get_antennaname()
     947    if ( calmode == 'nod' ):
     948        asaplog.push( 'Calibrating nod data.' )
     949        print_log()
     950        scal = calnod( scantab, scannos=scannos, verify=verify )
     951    elif ( calmode == 'quotient' ):
     952        asaplog.push( 'Calibrating using quotient.' )
     953        print_log()
     954        scal = scantab.auto_quotient( verify=verify )
     955    elif ( calmode == 'ps' ):
     956        asaplog.push( 'Calibrating %s position-switched data.' % antname )
     957        print_log()
     958        if ( antname.find( 'APEX' ) != -1 ):
     959            scal = apexcal( scantab, scannos, calmode, verify )
     960        elif ( antname.find( 'ALMA' ) != -1 or antname.find( 'OSF' ) != -1 ):
     961            scal = almacal( scantab, scannos, calmode, verify )
     962        else:
     963            scal = calps( scantab, scannos=scannos, verify=verify )
     964    elif ( calmode == 'fs' or calmode == 'fsotf' ):
     965        asaplog.push( 'Calibrating %s frequency-switched data.' % antname )
     966        print_log()
     967        if ( antname.find( 'APEX' ) != -1 ):
     968            scal = apexcal( scantab, scannos, calmode, verify )
     969        elif ( antname.find( 'ALMA' ) != -1 or antname.find( 'OSF' ) != -1 ):
     970            scal = almacal( scantab, scannos, calmode, verify )
     971        else:
     972            scal = calfs( scantab, scannos=scannos, verify=verify )
     973    elif ( calmode == 'otf' ):
     974        asaplog.push( 'Calibrating %s On-The-Fly data.' % antname )
     975        print_log()
     976        scal = almacal( scantab, scannos, calmode, verify )
     977    else:
     978        asaplog.push( 'No calibration.' )
     979        scal = scantab.copy()
     980
     981    return scal
     982
     983def apexcal( scantab, scannos=[], calmode='none', verify=False ):
     984    """
     985    Calibrate APEX data
     986
     987    Parameters:
     988        scantab:       scantable
     989        scannos:       list of scan number
     990        calmode:       calibration mode
     991
     992        verify:        verify calibration     
     993    """
     994    from asap._asap import stmath
     995    stm = stmath()
     996    antname = scantab.get_antennaname()
     997    ssub = scantab.get_scan( scannos )
     998    scal = scantable( stm.cwcal( ssub, calmode, antname ) )
     999    return scal
     1000
     1001def almacal( scantab, scannos=[], calmode='none', verify=False ):
     1002    """
     1003    Calibrate ALMA data
     1004
     1005    Parameters:
     1006        scantab:       scantable
     1007        scannos:       list of scan number
     1008        calmode:       calibration mode
     1009
     1010        verify:        verify calibration     
     1011    """
     1012    from asap._asap import stmath
     1013    stm = stmath()
     1014    ssub = scantab.get_scan( scannos )
     1015    scal = scantable( stm.almacal( ssub, calmode ) )
     1016    return scal
     1017
     1018def splitant(filename, outprefix='',overwrite=False):
     1019    """
     1020    Split Measurement set by antenna name, save data as a scantables,
     1021    and return a list of filename.
     1022    Notice this method can only be available from CASA.
     1023    Prameter
     1024       filename:    the name of Measurement set to be read.
     1025       outprefix:   the prefix of output scantable name.
     1026                    the names of output scantable will be
     1027                    outprefix.antenna1, outprefix.antenna2, ....
     1028                    If not specified, outprefix = filename is assumed.
     1029       overwrite    If the file should be overwritten if it exists.
     1030                    The default False is to return with warning
     1031                    without writing the output. USE WITH CARE.
     1032                   
     1033    """
     1034    # Import the table toolkit from CASA
     1035    try:
     1036       import casac
     1037    except ImportError:
     1038       if rcParams['verbose']:
     1039           #print "failed to load casa"
     1040           print_log()
     1041           asaplog.push("failed to load casa")
     1042           print_log('ERROR')
     1043       else: raise
     1044       return False
     1045    try:
     1046       tbtool = casac.homefinder.find_home_by_name('tableHome')
     1047       tb = tbtool.create()
     1048       tb2 = tbtool.create()
     1049    except:
     1050       if rcParams['verbose']:
     1051           #print "failed to load a table tool:\n", e
     1052           print_log()
     1053           asaplog.push("failed to load table tool")
     1054           print_log('ERROR')
     1055       else: raise
     1056       return False
     1057    # Check the input filename
     1058    if isinstance(filename, str):
     1059        import os.path
     1060        filename = os.path.expandvars(filename)
     1061        filename = os.path.expanduser(filename)
     1062        if not os.path.exists(filename):
     1063            s = "File '%s' not found." % (filename)
     1064            if rcParams['verbose']:
     1065                print_log()
     1066                asaplog.push(s)
     1067                print_log('ERROR')
     1068                return
     1069            raise IOError(s)
     1070        # check if input file is MS
     1071        if not os.path.isdir(filename) \
     1072               or not os.path.exists(filename+'/ANTENNA') \
     1073               or not os.path.exists(filename+'/table.f1'):
     1074            s = "File '%s' is not a Measurement set." % (filename)
     1075            if rcParams['verbose']:
     1076                print_log()
     1077                asaplog.push(s)
     1078                print_log('ERROR')
     1079                return
     1080            raise IOError(s)
     1081    else:
     1082        s = "The filename should be string. "
     1083        if rcParams['verbose']:
     1084            print_log()
     1085            asaplog.push(s)
     1086            print_log('ERROR')
     1087            return
     1088        raise TypeError(s)
     1089    # Check out put file name
     1090    outname=''
     1091    if len(outprefix) > 0: prefix=outprefix+'.'
     1092    else:
     1093        prefix=filename.rstrip('/')
     1094    # Now do the actual splitting.
     1095    outfiles=[]
     1096    tb.open(tablename=filename+'/ANTENNA',nomodify=True)
     1097    nant=tb.nrows()
     1098    antnames=tb.getcol('NAME',0,nant,1)
     1099    antpos=tb.getcol('POSITION',0,nant,1).transpose()
     1100    tb.close()
     1101    tb.open(tablename=filename,nomodify=True)
     1102    ant1=tb.getcol('ANTENNA1',0,-1,1)
     1103    tb.close()
     1104    for antid in set(ant1):
     1105        scan=scantable(filename,average=False,getpt=True,antenna=int(antid))
     1106        outname=prefix+antnames[antid]+'.asap'
     1107        scan.save(outname,format='ASAP',overwrite=overwrite)
     1108        del scan
     1109        outfiles.append(outname)
     1110    del tb, tb2
     1111    return outfiles
     1112
     1113def _array2dOp( scan, value, mode="ADD", tsys=False ):
     1114    """
     1115    This function is workaround on the basic operation of scantable
     1116    with 2 dimensional float list.
     1117
     1118    scan:    scantable operand
     1119    value:   float list operand
     1120    mode:    operation mode (ADD, SUB, MUL, DIV)
     1121    tsys:    if True, operate tsys as well
     1122    """
     1123    nrow = scan.nrow()
     1124    s = None
     1125    if len( value ) == 1:
     1126        from asap._asap import stmath
     1127        stm = stmath()
     1128        s = scantable( stm._arrayop( scan.copy(), value[0], mode, tsys ) )
     1129        del stm
     1130    elif len( value ) != nrow:
     1131        asaplog.push( 'len(value) must be 1 or conform to scan.nrow()' )
     1132        print_log( 'ERROR' )
     1133    else:
     1134        from asap._asap import stmath
     1135        stm = stmath()
     1136        # insitu must be True
     1137        stm._setinsitu( True )
     1138        s = scan.copy()
     1139        sel = selector()
     1140        for irow in range( nrow ):
     1141            sel.set_rows( irow )
     1142            s.set_selection( sel )
     1143            if len( value[irow] ) == 1:
     1144                stm._unaryop( s, value[irow][0], mode, tsys )
     1145            else:
     1146                stm._arrayop( s, value[irow], mode, tsys, 'channel' )
     1147            s.set_selection()
     1148            sel.reset()
     1149        del sel
     1150        del stm
     1151    return s
     1152
     1153           
     1154           
Note: See TracChangeset for help on using the changeset viewer.