source: trunk/python/scantable.py @ 513

Last change on this file since 513 was 513, checked in by mar637, 19 years ago

*Moved most asapmath functions to scantable member functions. Only ones taking more than one scantable as input remain.

  • added asap._is_sequence_or_number to check input args
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 54.9 KB
Line 
1from asap._asap import sdtable
2from asap import rcParams
3from numarray import ones,zeros
4import sys
5
6class scantable(sdtable):
7    """
8        The ASAP container for scans
9    """
10   
11    def __init__(self, filename, average=None, unit=None):
12        """
13        Create a scantable from a saved one or make a reference
14        Parameters:
15            filename:    the name of an asap table on disk
16                         or
17                         the name of a rpfits/sdfits/ms file
18                         (integrations within scans are auto averaged
19                         and the whole file is read)
20                         or
21                         [advanced] a reference to an existing
22                         scantable
23            average:     average all integrations withinb a scan on read.
24                         The default (True) is taken from .asaprc.
25            unit:         brightness unit; must be consistent with K or Jy.
26                         Over-rides the default selected by the reader
27                         (input rpfits/sdfits/ms) or replaces the value
28                         in existing scantables
29        """       
30        if average is None or type(average) is not bool:
31            average = rcParams['scantable.autoaverage']           
32
33        varlist = vars()
34        self._vb = rcParams['verbose']
35        self._p = None
36
37        if isinstance(filename,sdtable):
38            sdtable.__init__(self, filename)           
39            if unit is not None:
40                self.set_fluxunit(unit)                       
41        else:
42            import os.path
43            if not os.path.exists(filename):
44                print "File '%s' not found." % (filename)
45                return
46            filename = os.path.expandvars(filename)
47            if os.path.isdir(filename):
48                # crude check if asap table
49                if os.path.exists(filename+'/table.info'):
50                    sdtable.__init__(self, filename)
51                    if unit is not None:
52                        self.set_fluxunit(unit)                       
53                else:
54                    print "The given file '%s'is not a valid asap table." % (filename)
55                    return
56            else:           
57                from asap._asap import sdreader
58                ifSel = -1
59                beamSel = -1
60                r = sdreader(filename,ifSel,beamSel)
61                print 'Importing data...'
62                r._read([-1])
63                tbl = r._getdata()
64                if unit is not None:
65                    tbl.set_fluxunit(unit)
66                if average:
67                    from asap._asap import average as _av
68                    print 'Auto averaging integrations...'
69                    tbl2 = _av((tbl,),(),True,'none')
70                    sdtable.__init__(self,tbl2)
71                    del tbl2
72                else:
73                    sdtable.__init__(self,tbl)
74                del r,tbl
75                self._add_history("scantable", varlist)
76
77    def save(self, name=None, format=None, stokes=False, overwrite=False):
78        """
79        Store the scantable on disk. This can be an asap (aips++) Table, SDFITS,
80        Image FITS or MS2 format.
81        Parameters:
82            name:        the name of the outputfile. For format="FITS" this
83                         is the directory file name into which all the files
84                         will be written (default is 'asap_FITS'). For format
85                         "ASCII" this is the root file name (data in 'name'.txt
86                         and header in 'name'_header.txt)
87            format:      an optional file format. Default is ASAP.
88                         Allowed are - 'ASAP' (save as ASAP [aips++] Table),
89                                       'SDFITS' (save as SDFITS file)
90                                       'FITS' (saves each row as a FITS Image)
91                                       'ASCII' (saves as ascii text file)
92                                       'MS2' (saves as an aips++
93                                              MeasurementSet V2)
94            stokes:      Convert to Stokes parameters (only available
95                         currently with FITS and ASCII formats.
96                         Default is False.
97            overwrite:   If the file should be overwritten if it exists.
98                         The default False is to return with warning
99                         without writing the output. USE WITH CARE.
100        Example:
101            scan.save('myscan.asap')
102            scan.save('myscan.sdfits','SDFITS')
103        """
104        from os import path
105        if format is None: format = rcParams['scantable.save']
106        suffix = '.'+format.lower()
107        if name is None or name =="":
108            name = 'scantable'+suffix
109            print "No filename given. Using default name %s..." % name
110        name = path.expandvars(name)
111        if path.isfile(name) or path.isdir(name):
112            if not overwrite:
113                print "File %s already exists." % name
114                return
115        format2 = format.upper()
116        if format2 == 'ASAP':
117            self._save(name)
118        else:
119            from asap._asap import sdwriter as _sw
120            w = _sw(format2)
121            w.write(self, name, stokes)
122        return
123
124    def copy(self):
125        """
126        Return a copy of this scantable.
127        Parameters:
128            none
129        Example:
130            copiedscan = scan.copy()
131        """
132        sd = scantable(sdtable._copy(self))
133        return sd
134
135    def get_scan(self, scanid=None):
136        """
137        Return a specific scan (by scanno) or collection of scans (by
138        source name) in a new scantable.
139        Parameters:
140            scanid:    a (list of) scanno or a source name, unix-style
141                       patterns are accepted for source name matching, e.g.
142                       '*_R' gets all 'ref scans
143        Example:
144            # get all scans containing the source '323p459'
145            newscan = scan.get_scan('323p459')
146            # get all 'off' scans
147            refscans = scan.get_scan('*_R')
148            # get a susbset of scans by scanno (as listed in scan.summary())
149            newscan = scan.get_scan([0,2,7,10])
150        """
151        if scanid is None:
152            print "Please specify a scan no or name to retrieve from the scantable"
153        try:
154            if type(scanid) is str:
155                s = sdtable._getsource(self,scanid)
156                return scantable(s)
157            elif type(scanid) is int:
158                s = sdtable._getscan(self,[scanid])
159                return scantable(s)
160            elif type(scanid) is list:
161                s = sdtable._getscan(self,scanid)
162                return scantable(s)
163            else:
164                print "Illegal scanid type, use 'int' or 'list' if ints."
165        except RuntimeError:
166            print "Couldn't find any match."
167
168    def __str__(self):
169        return sdtable._summary(self,True)
170
171    def summary(self,filename=None, verbose=None):
172        """
173        Print a summary of the contents of this scantable.
174        Parameters:
175            filename:    the name of a file to write the putput to
176                         Default - no file output
177            verbose:     print extra info such as the frequency table
178                         The default (False) is taken from .asaprc
179        """
180        info = sdtable._summary(self, verbose)
181        if verbose is None: verbose = rcParams['scantable.verbosesummary']
182        if filename is not None:
183            if filename is "":
184                filename = 'scantable_summary.txt'
185            from os.path import expandvars, isdir
186            filename = expandvars(filename)
187            if not isdir(filename):
188                data = open(filename, 'w')
189                data.write(info)
190                data.close()
191            else:
192                print "Illegal file name '%s'." % (filename)
193        print info
194           
195    def set_cursor(self, thebeam=0,theif=0,thepol=0):
196        """
197        Set the spectrum for individual operations.
198        Parameters:
199            thebeam,theif,thepol:    a number
200        Example:
201            scan.set_cursor(0,0,1)
202            pol1sig = scan.stats(all=False) # returns std dev for beam=0
203                                            # if=0, pol=1
204        """
205        varlist = vars()
206        self.setbeam(thebeam)
207        self.setpol(thepol)
208        self.setif(theif)
209        self._add_history("set_cursor",varlist)
210        return
211
212    def get_cursor(self):
213        """
214        Return/print a the current 'cursor' into the Beam/IF/Pol cube.
215        Parameters:
216            none
217        Returns:
218            a list of values (currentBeam,currentIF,currentPol)
219        Example:
220            none
221        """
222        i = self.getbeam()
223        j = self.getif()
224        k = self.getpol()
225        if self._vb:
226            print "--------------------------------------------------"
227            print " Cursor position"
228            print "--------------------------------------------------"
229            out = 'Beam=%d IF=%d Pol=%d ' % (i,j,k)
230            print out
231        return i,j,k
232
233    def stats(self, stat='stddev', mask=None, allaxes=None):
234        """
235        Determine the specified statistic of the current beam/if/pol
236        Takes a 'mask' as an optional parameter to specify which
237        channels should be excluded.
238        Parameters:
239            stat:    'min', 'max', 'sumsq', 'sum', 'mean'
240                     'var', 'stddev', 'avdev', 'rms', 'median'
241            mask:    an optional mask specifying where the statistic
242                     should be determined.
243            allaxes: if True apply to all spectra. Otherwise
244                     apply only to the selected (beam/pol/if)spectra only.
245                     The default is taken from .asaprc (True if none)
246        Example:
247            scan.set_unit('channel')
248            msk = scan.create_mask([100,200],[500,600])
249            scan.stats(stat='mean', mask=m)
250        """
251        if allaxes is None: allaxes = rcParams['scantable.allaxes']
252        from asap._asap import stats as _stats
253        from numarray import array,zeros,Float
254        if mask == None:
255            mask = ones(self.nchan())
256        axes = ['Beam','IF','Pol','Time']
257
258        beamSel,IFSel,polSel = (self.getbeam(),self.getif(),self.getpol())
259        if allaxes:
260            n = self.nbeam()*self.nif()*self.npol()*self.nrow()
261            shp = [self.nbeam(),self.nif(),self.npol(),self.nrow()]
262            arr = array(zeros(n),shape=shp,type=Float)
263
264            for i in range(self.nbeam()):
265                self.setbeam(i)
266                for j in range(self.nif()):
267                    self.setif(j)
268                    for k in range(self.npol()):
269                        self.setpol(k)
270                        arr[i,j,k,:] = _stats(self,mask,stat,-1)
271            retval = {'axes': axes, 'data': arr, 'cursor':None}
272            tm = [self._gettime(val) for val in range(self.nrow())]
273            if self._vb:
274                self._print_values(retval,stat,tm)
275            self.setbeam(beamSel)
276            self.setif(IFSel)
277            self.setpol(polSel)
278            return retval
279
280        else:
281            statval = _stats(self,mask,stat,-1)
282            out = ''
283            for l in range(self.nrow()):
284                tm = self._gettime(l)
285                out += 'Time[%s]:\n' % (tm)
286                if self.nbeam() > 1: out +=  ' Beam[%d] ' % (beamSel)
287                if self.nif() > 1: out +=  ' IF[%d] ' % (IFSel)
288                if self.npol() > 1: out +=  ' Pol[%d] ' % (polSel)
289                out += '= %3.3f\n' % (statval[l])
290                out +=  "--------------------------------------------------\n"
291
292            if self._vb:
293                print "--------------------------------------------------"
294                print " ",stat
295                print "--------------------------------------------------"
296                print out
297            retval = {'axes': axes, 'data': array(statval), 'cursor':(i,j,k)}
298            return retval
299
300    def stddev(self,mask=None, allaxes=None):
301        """
302        Determine the standard deviation of the current beam/if/pol
303        Takes a 'mask' as an optional parameter to specify which
304        channels should be excluded.
305        Parameters:
306            mask:    an optional mask specifying where the standard
307                     deviation should be determined.
308            allaxes: optional flag to show all or a cursor selected
309                     spectrum of Beam/IF/Pol. Default is all or taken
310                     from .asaprc
311
312        Example:
313            scan.set_unit('channel')
314            msk = scan.create_mask([100,200],[500,600])
315            scan.stddev(mask=m)
316        """
317        if allaxes is None: allaxes = rcParams['scantable.allaxes']
318        return self.stats(stat='stddev',mask=mask, allaxes=allaxes);
319
320    def get_tsys(self, allaxes=None):
321        """
322        Return the System temperatures.
323        Parameters:
324           allaxes:     if True apply to all spectra. Otherwise
325                        apply only to the selected (beam/pol/if)spectra only.
326                        The default is taken from .asaprc (True if none)
327        Returns:
328            a list of Tsys values.
329        """
330        if allaxes is None: allaxes = rcParams['scantable.allaxes']
331        from numarray import array,zeros,Float
332        axes = ['Beam','IF','Pol','Time']
333
334        if allaxes:
335            n = self.nbeam()*self.nif()*self.npol()*self.nrow()
336            shp = [self.nbeam(),self.nif(),self.npol(),self.nrow()]
337            arr = array(zeros(n),shape=shp,type=Float)
338
339            for i in range(self.nbeam()):
340                self.setbeam(i)
341                for j in range(self.nif()):
342                    self.setif(j)
343                    for k in range(self.npol()):
344                        self.setpol(k)
345                        arr[i,j,k,:] = self._gettsys()
346            retval = {'axes': axes, 'data': arr, 'cursor':None}
347            tm = [self._gettime(val) for val in range(self.nrow())]
348            if self._vb:
349                self._print_values(retval,'Tsys',tm)
350            return retval
351
352        else:
353            i,j,k = (self.getbeam(),self.getif(),self.getpol())
354            statval = self._gettsys()
355            out = ''
356            for l in range(self.nrow()):
357                tm = self._gettime(l)
358                out += 'Time[%s]:\n' % (tm)
359                if self.nbeam() > 1: out +=  ' Beam[%d] ' % (i)
360                if self.nif() > 1: out +=  ' IF[%d] ' % (j)
361                if self.npol() > 1: out +=  ' Pol[%d] ' % (k)
362                out += '= %3.3f\n' % (statval[l])
363                out +=  "--------------------------------------------------\n"
364
365            if self._vb:
366                print "--------------------------------------------------"
367                print " TSys"
368                print "--------------------------------------------------"
369                print out
370            retval = {'axes': axes, 'data': array(statval), 'cursor':(i,j,k)}
371            return retval
372       
373    def get_time(self, row=-1):
374        """
375        Get a list of time stamps for the observations.
376        Return a string for each integration in the scantable.
377        Parameters:
378            row:    row no of integration. Default -1 return all rows
379        Example:
380            none
381        """
382        out = []
383        if row == -1:
384            for i in range(self.nrow()):
385                out.append(self._gettime(i))
386            return out
387        else:
388            if row < self.nrow():
389                return self._gettime(row)
390
391    def set_unit(self, unit='channel'):
392        """
393        Set the unit for all following operations on this scantable
394        Parameters:
395            unit:    optional unit, default is 'channel'
396                     one of '*Hz','km/s','channel', ''
397        """
398        varlist = vars()
399        if unit in ['','pixel', 'channel']:
400            unit = ''
401        inf = list(self._getcoordinfo())
402        inf[0] = unit
403        self._setcoordinfo(inf)
404        if self._p: self.plot()
405        self._add_history("set_unit",varlist)
406
407    def set_instrument(self, instr):
408        """
409        Set the instrument for subsequent processing
410        Parameters:
411            instr:    Select from 'ATPKSMB', 'ATPKSHOH', 'ATMOPRA',
412                      'DSS-43' (Tid), 'CEDUNA', and 'HOBART'
413        """
414        self._setInstrument(instr)
415        self._add_history("set_instument",vars())
416
417    def set_doppler(self, doppler='RADIO'):
418        """
419        Set the doppler for all following operations on this scantable.
420        Parameters:
421            doppler:    One of 'RADIO', 'OPTICAL', 'Z', 'BETA', 'GAMMA'
422        """
423        varlist = vars()
424        inf = list(self._getcoordinfo())
425        inf[2] = doppler
426        self._setcoordinfo(inf)
427        if self._p: self.plot()
428        self._add_history("set_doppler",vars())
429 
430    def set_freqframe(self, frame=None):
431        """
432        Set the frame type of the Spectral Axis.
433        Parameters:
434            frame:   an optional frame type, default 'LSRK'.
435        Examples:
436            scan.set_freqframe('BARY')
437        """
438        if frame is None: frame = rcParams['scantable.freqframe']
439        varlist = vars()
440        valid = ['REST','TOPO','LSRD','LSRK','BARY', \
441                   'GEO','GALACTO','LGROUP','CMB']
442        if frame in valid:
443            inf = list(self._getcoordinfo())
444            inf[1] = frame
445            self._setcoordinfo(inf)
446            self._add_history("set_freqframe",varlist)
447        else:
448            print "Please specify a valid freq type. Valid types are:\n",valid
449           
450    def get_unit(self):
451        """
452        Get the default unit set in this scantable
453        Parameters:
454        Returns:
455            A unit string
456        """
457        inf = self._getcoordinfo()
458        unit = inf[0]
459        if unit == '': unit = 'channel'
460        return unit
461
462    def get_abcissa(self, rowno=0):
463        """
464        Get the abcissa in the current coordinate setup for the currently
465        selected Beam/IF/Pol
466        Parameters:
467            rowno:    an optional row number in the scantable. Default is the
468                      first row, i.e. rowno=0
469        Returns:
470            The abcissa values and it's format string (as a dictionary)
471        """
472        abc = self._getabcissa(rowno)
473        lbl = self._getabcissalabel(rowno)       
474        return abc, lbl
475        #return {'abcissa':abc,'label':lbl}
476
477    def create_mask(self, *args, **kwargs):
478        """
479        Compute and return a mask based on [min,max] windows.
480        The specified windows are to be INCLUDED, when the mask is
481        applied.
482        Parameters:
483            [min,max],[min2,max2],...
484                Pairs of start/end points specifying the regions
485                to be masked
486            invert:     optional argument. If specified as True,
487                        return an inverted mask, i.e. the regions
488                        specified are EXCLUDED
489            row:        create the mask using the specified row for
490                        unit conversions, default is row=0
491                        only necessary if frequency varies over rows.
492        Example:
493            scan.set_unit('channel')
494
495            a)
496            msk = scan.set_mask([400,500],[800,900])
497            # masks everything outside 400 and 500
498            # and 800 and 900 in the unit 'channel'
499
500            b)
501            msk = scan.set_mask([400,500],[800,900], invert=True)
502            # masks the regions between 400 and 500
503            # and 800 and 900 in the unit 'channel'
504           
505        """
506        row = 0
507        if kwargs.has_key("row"):
508            row = kwargs.get("row")
509        data = self._getabcissa(row)
510        u = self._getcoordinfo()[0]
511        if self._vb:
512            if u == "": u = "channel"
513            print "The current mask window unit is", u
514        n = self.nchan()
515        msk = zeros(n)
516        for window in args:
517            if (len(window) != 2 or window[0] > window[1] ):
518                print "A window needs to be defined as [min,max]"
519                return
520            for i in range(n):
521                if data[i] >= window[0] and data[i] < window[1]:
522                    msk[i] = 1
523        if kwargs.has_key('invert'):
524            if kwargs.get('invert'):
525                from numarray import logical_not
526                msk = logical_not(msk)           
527        return msk
528   
529    def get_restfreqs(self):
530        """
531        Get the restfrequency(s) stored in this scantable.
532        The return value(s) are always of unit 'Hz'
533        Parameters:
534            none
535        Returns:
536            a list of doubles
537        """
538        return list(self._getrestfreqs())
539
540    def lines(self):
541        """
542        Print the list of known spectral lines
543        """
544        sdtable._lines(self)
545
546    def set_restfreqs(self, freqs=None, unit='Hz', lines=None, source=None,
547                      theif=None):
548        """
549        Select the restfrequency for the specified source and IF OR
550        replace for all IFs.  If the 'freqs' argument holds a scalar,
551        then that rest frequency will be applied to the selected
552        data (and added to the list of available rest frequencies).
553        In this way, you can set a rest frequency for each
554        source and IF combination.   If the 'freqs' argument holds
555        a vector, then it MUST be of length the number of IFs
556        (and the available restfrequencies will be replaced by
557        this vector).  In this case, *all* data ('source' and
558        'theif' are ignored) have the restfrequency set per IF according
559        to the corresponding value you give in the 'freqs' vector. 
560        E.g. 'freqs=[1e9,2e9]'  would mean IF 0 gets restfreq 1e9 and
561        IF 1 gets restfreq 2e9.
562
563        You can also specify the frequencies via known line names
564        in the argument 'lines'.  Use 'freqs' or 'lines'.  'freqs'
565        takes precedence. See the list of known names via function
566        scantable.lines()
567        Parameters:
568            freqs:   list of rest frequencies
569            unit:    unit for rest frequency (default 'Hz')
570            lines:   list of known spectral lines names (alternative to freqs).
571                     See possible list via scantable.lines()
572            source:  Source name (blank means all)
573            theif:   IF (-1 means all)
574        Example:
575            scan.set_restfreqs(freqs=1.4e9, source='NGC253', theif=2)
576            scan.set_restfreqs(freqs=[1.4e9,1.67e9])
577        """
578        varlist = vars()
579        if source is None:
580            source = ""
581        if theif is None:
582            theif = -1
583        t = type(freqs)
584        if t is int or t is float:
585           freqs = [freqs]
586        if freqs is None:
587           freqs = []
588        t = type(lines)
589        if t is str:
590           lines = [lines]
591        if lines is None:
592           lines = []
593        sdtable._setrestfreqs(self, freqs, unit, lines, source, theif)
594        self._add_history("set_restfreqs", varlist)
595       
596
597
598    def flag_spectrum(self, thebeam, theif, thepol):
599        """
600        This flags a selected spectrum in the scan 'for good'.
601        USE WITH CARE - not reversible.
602        Use masks for non-permanent exclusion of channels.
603        Parameters:
604            thebeam,theif,thepol:    all have to be explicitly
605                                     specified
606        Example:
607            scan.flag_spectrum(0,0,1)
608            flags the spectrum for Beam=0, IF=0, Pol=1
609        """
610        if (thebeam < self.nbeam() and
611            theif < self.nif() and
612            thepol < self.npol()):
613            sdtable.setbeam(self, thebeam)
614            sdtable.setif(self, theif)
615            sdtable.setpol(self, thepol)
616            sdtable._flag(self)
617            self._add_history("flag_spectrum", vars())
618        else:
619            print "Please specify a valid (Beam/IF/Pol)"
620        return
621
622    def plot(self, what='spectrum',col='Pol', panel=None):
623        """
624        Plot the spectra contained in the scan. Alternatively you can also
625        Plot Tsys vs Time
626        Parameters:
627            what:     a choice of 'spectrum' (default) or 'tsys'
628            col:      which out of Beams/IFs/Pols should be colour stacked
629            panel:    set up multiple panels, currently not working.
630        """
631        print "Warning! Not fully functional. Use plotter.plot() instead"
632       
633        validcol = {'Beam':self.nbeam(),'IF':self.nif(),'Pol':self.npol()}
634
635        validyax = ['spectrum','tsys']
636        from asap.asaplot import ASAPlot
637        if not self._p:
638            self._p = ASAPlot()
639            #print "Plotting not enabled"
640            #return
641        if self._p.is_dead:
642            del self._p
643            self._p = ASAPlot()
644        npan = 1
645        x = None
646        if what == 'tsys':
647            n = self.nrow()
648            if n < 2:
649                print "Only one integration. Can't plot."
650                return
651        self._p.hold()
652        self._p.clear()
653        if panel == 'Time':
654            npan = self.nrow()
655            self._p.set_panels(rows=npan)
656        xlab,ylab,tlab = None,None,None
657        self._vb = False
658        sel = self.get_cursor()       
659        for i in range(npan):
660            if npan > 1:
661                self._p.subplot(i)
662            for j in range(validcol[col]):
663                x = None
664                y = None
665                m = None
666                tlab = self._getsourcename(i)
667                import re
668                tlab = re.sub('_S','',tlab)
669                if col == 'Beam':
670                    self.setbeam(j)
671                elif col == 'IF':
672                    self.setif(j)
673                elif col == 'Pol':
674                    self.setpol(j)
675                if what == 'tsys':
676                    x = range(self.nrow())
677                    xlab = 'Time [pixel]'
678                    m = list(ones(len(x)))
679                    y = []
680                    ylab = r'$T_{sys}$'
681                    for k in range(len(x)):
682                        y.append(self._gettsys(k))
683                else:
684                    x,xlab = self.get_abcissa(i)
685                    y = self._getspectrum(i)
686                    ylab = r'Flux'
687                    m = self._getmask(i)
688                llab = col+' '+str(j)
689                self._p.set_line(label=llab)
690                self._p.plot(x,y,m)
691            self._p.set_axes('xlabel',xlab)
692            self._p.set_axes('ylabel',ylab)
693            self._p.set_axes('title',tlab)
694        self._p.release()
695        self.set_cursor(sel[0],sel[1],sel[2])
696        self._vb = rcParams['verbose']
697        return
698
699        print out
700
701    def _print_values(self, dat, label='', timestamps=[]):
702        d = dat['data']
703        a = dat['axes']
704        shp = d.getshape()
705        out = ''
706        for i in range(shp[3]):
707            out += '%s [%s]:\n' % (a[3],timestamps[i])
708            t = d[:,:,:,i]
709            for j in range(shp[0]):
710                if shp[0] > 1: out +=  ' %s[%d] ' % (a[0],j)
711                for k in range(shp[1]):
712                    if shp[1] > 1: out +=  ' %s[%d] ' % (a[1],k)
713                    for l in range(shp[2]):
714                        if shp[2] > 1: out +=  ' %s[%d] ' % (a[2],l)
715                        out += '= %3.3f\n' % (t[j,k,l])
716            out += "-"*80
717            out += "\n"
718        print "-"*80
719        print " ", label
720        print "-"*80
721        print out
722
723    def history(self):
724        hist = list(self._gethistory())
725        print "-"*80
726        for h in hist:
727            if h.startswith("---"):
728                print h
729            else:
730                items = h.split("##")
731                date = items[0]
732                func = items[1]
733                items = items[2:]
734                print date           
735                print "Function: %s\n  Parameters:" % (func)
736                for i in items:
737                    s = i.split("=")
738                    print "   %s = %s" % (s[0],s[1])
739                print "-"*80
740        return
741
742    #
743    # Maths business
744    #
745
746    def average_time(self, mask=None, scanav=True, weight=None):
747        """
748        Return the (time) average of a scan, or apply it 'insitu'.
749        Note:
750            in channels only
751            The cursor of the output scan is set to 0.
752        Parameters:
753            one scan or comma separated  scans
754            mask:     an optional mask (only used for 'var' and 'tsys'
755                      weighting)
756            scanav:   True (default) averages each scan separately
757                      False averages all scans together,
758            weight:   Weighting scheme. 'none' (default), 'var' (variance
759                      weighted), 'tsys'
760
761        Example:
762            # time average the scantable without using a mask
763            newscan = scan.average_time()           
764        """
765        varlist = vars()
766        if weight is None: weight = 'none'
767        if mask is None: mask = ()
768        from asap._asap import average as _av       
769        s = scantable(_av((self,), mask, scanav, weight))
770        s._add_history("average_time",varlist)
771        return s
772       
773    def convert_flux(self, jyperk=None, eta=None, d=None, insitu=None,
774                     allaxes=None):
775        """
776        Return a scan where all spectra are converted to either
777        Jansky or Kelvin depending upon the flux units of the scan table.
778        By default the function tries to look the values up internally.
779        If it can't find them (or if you want to over-ride), you must
780        specify EITHER jyperk OR eta (and D which it will try to look up
781        also if you don't set it). jyperk takes precedence if you set both.
782        Parameters:
783            jyperk:      the Jy / K conversion factor
784            eta:         the aperture efficiency
785            d:           the geomtric diameter (metres)
786            insitu:      if False a new scantable is returned.
787                         Otherwise, the scaling is done in-situ
788                         The default is taken from .asaprc (False)
789            allaxes:         if True apply to all spectra. Otherwise
790                         apply only to the selected (beam/pol/if)spectra only
791                         The default is taken from .asaprc (True if none)
792        """
793        if allaxes is None: allaxes = rcParams['scantable.allaxes']
794        if insitu is None: insitu = rcParams['insitu']
795        varlist = vars()
796        if jyperk is None: jyperk = -1.0
797        if d is None: d = -1.0
798        if eta is None: eta = -1.0
799        if not insitu:
800            from asap._asap import convertflux as _convert
801            s = scantable(_convert(self, d, eta, jyperk, allaxes))
802            s._add_history("convert_flux", varlist)
803            return s
804        else:
805            from asap._asap import convertflux_insitu as _convert
806            _convert(self, d, eta, jyperk, allaxes)
807            self._add_history("convert_flux", varlist)
808            return
809
810    def gain_el(self, poly=None, filename="", method="linear",
811                insitu=None, allaxes=None):
812        """
813        Return a scan after applying a gain-elevation correction.
814        The correction can be made via either a polynomial or a
815        table-based interpolation (and extrapolation if necessary).
816        You specify polynomial coefficients, an ascii table or neither.
817        If you specify neither, then a polynomial correction will be made
818        with built in coefficients known for certain telescopes (an error
819        will occur if the instrument is not known).
820        The data and Tsys are *divided* by the scaling factors.
821        Parameters:
822            poly:        Polynomial coefficients (default None) to compute a
823                         gain-elevation correction as a function of
824                         elevation (in degrees).
825            filename:    The name of an ascii file holding correction factors.
826                         The first row of the ascii file must give the column
827                         names and these MUST include columns
828                         "ELEVATION" (degrees) and "FACTOR" (multiply data
829                         by this) somewhere.
830                         The second row must give the data type of the
831                         column. Use 'R' for Real and 'I' for Integer.
832                         An example file would be
833                         (actual factors are arbitrary) :
834
835                         TIME ELEVATION FACTOR
836                         R R R
837                         0.1 0 0.8
838                         0.2 20 0.85
839                         0.3 40 0.9
840                         0.4 60 0.85
841                         0.5 80 0.8
842                         0.6 90 0.75
843            method:      Interpolation method when correcting from a table.
844                         Values are  "nearest", "linear" (default), "cubic"
845                         and "spline"
846            insitu:      if False a new scantable is returned.
847                         Otherwise, the scaling is done in-situ
848                         The default is taken from .asaprc (False)
849            allaxes:     If True apply to all spectra. Otherwise
850                         apply only to the selected (beam/pol/if) spectra only
851                         The default is taken from .asaprc (True if none)
852        """
853
854        if allaxes is None: allaxes = rcParams['scantable.allaxes']
855        if insitu is None: insitu = rcParams['insitu']
856        varlist = vars()
857        if poly is None:
858           poly = ()
859        from os.path import expandvars
860        filename = expandvars(filename)
861        if not insitu:
862            from asap._asap import gainel as _gainEl
863            s = scantable(_gainEl(self, poly, filename, method, allaxes))
864            s._add_history("gain_el", varlist)
865            return s
866        else:
867            from asap._asap import gainel_insitu as _gainEl
868            _gainEl(self, poly, filename, method, allaxes)
869            self._add_history("gain_el", varlist)
870            return
871   
872    def freq_align(self, reftime=None, method='cubic', perif=False,
873                   insitu=None):
874        """
875        Return a scan where all rows have been aligned in frequency/velocity.
876        The alignment frequency frame (e.g. LSRK) is that set by function
877        set_freqframe.
878        Parameters:
879            reftime:     reference time to align at. By default, the time of
880                         the first row of data is used.
881            method:      Interpolation method for regridding the spectra.
882                         Choose from "nearest", "linear", "cubic" (default)
883                         and "spline"
884            perif:       Generate aligners per freqID (no doppler tracking) or
885                         per IF (scan-based doppler tracking)
886            insitu:      if False a new scantable is returned.
887                         Otherwise, the scaling is done in-situ
888                         The default is taken from .asaprc (False)
889        """
890        if insitu is None: insitu = rcParams['insitu']
891        varlist = vars()
892        if reftime is None: reftime = ''
893        perfreqid = not perif
894        if not insitu:
895            from asap._asap import freq_align as _align
896            s = scantable(_align(self, reftime, method, perfreqid))
897            s._add_history("freq_align", varlist)
898            return s
899        else:
900            from asap._asap import freq_align_insitu as _align
901            _align(self, reftime, method, perfreqid)
902            self._add_history("freq_align", varlist)
903            return
904
905    def opacity(self, tau, insitu=None, allaxes=None):
906        """
907        Apply an opacity correction. The data
908        and Tsys are multiplied by the correction factor.
909        Parameters:
910            tau:         Opacity from which the correction factor is
911                         exp(tau*ZD)
912                         where ZD is the zenith-distance
913            insitu:      if False a new scantable is returned.
914                         Otherwise, the scaling is done in-situ
915                         The default is taken from .asaprc (False)
916            allaxes:     if True apply to all spectra. Otherwise
917                         apply only to the selected (beam/pol/if)spectra only
918                         The default is taken from .asaprc (True if none)
919        """
920        if allaxes is None: allaxes = rcParams['scantable.allaxes']
921        if insitu is None: insitu = rcParams['insitu']
922        varlist = vars()
923        if not insitu:
924            from asap._asap import opacity as _opacity
925            s = scantable(_opacity(self, tau, allaxes))
926            s._add_history("opacity", varlist)
927            return s
928        else:
929            from asap._asap import opacity_insitu as _opacity
930            _opacity(self, tau, allaxes)
931            self._add_history("opacity", varlist)
932            return
933
934    def bin(self, width=5, insitu=None):
935        """
936        Return a scan where all spectra have been binned up.
937            width:       The bin width (default=5) in pixels
938            insitu:      if False a new scantable is returned.
939                         Otherwise, the scaling is done in-situ
940                         The default is taken from .asaprc (False)
941        """
942        if insitu is None: insitu = rcParams['insitu']
943        varlist = vars()
944        if not insitu:
945            from asap._asap import bin as _bin
946            s = scantable(_bin(self, width))
947            s._add_history("bin",varlist)
948            return s
949        else:
950            from asap._asap import bin_insitu as _bin
951            _bin(self, width)
952            self._add_history("bin",varlist)
953            return
954
955   
956    def resample(self, width=5, method='cubic', insitu=None):
957        """
958        Return a scan where all spectra have been binned up
959            width:       The bin width (default=5) in pixels
960            method:      Interpolation method when correcting from a table.
961                         Values are  "nearest", "linear", "cubic" (default)
962                         and "spline"
963            insitu:      if False a new scantable is returned.
964                         Otherwise, the scaling is done in-situ
965                         The default is taken from .asaprc (False)
966        """
967        if insitu is None: insitu = rcParams['insitu']
968        varlist = vars()
969        if not insitu:
970            from asap._asap import resample as _resample
971            s = scantable(_resample(self, method, width))
972            s._add_history("resample",varlist)
973            return s
974        else:
975            from asap._asap import resample_insitu as _resample
976            _resample(self, method, width)
977            self._add_history("resample",varlist)
978            return
979
980    def average_pol(self, mask=None, weight='none', insitu=None):
981        """
982        Average the Polarisations together.
983        The polarisation cursor of the output scan is set to 0
984        Parameters:
985            scan:        The scantable
986            mask:        An optional mask defining the region, where the
987                         averaging will be applied. The output will have all
988                         specified points masked.
989            weight:      Weighting scheme. 'none' (default), or
990                         'var' (variance weighted)
991            insitu:      if False a new scantable is returned.
992                         Otherwise, the scaling is done in-situ
993                         The default is taken from .asaprc (False)
994        Example:
995            polav = average_pols(myscan)
996        """
997        if insitu is None: insitu = rcParams['insitu']
998        varlist = vars()
999
1000        if mask is None:
1001            mask = ()
1002        if not insitu:
1003            from asap._asap import averagepol as _avpol
1004            s = scantable(_avpol(self, mask, weight))
1005            s._add_history("average_pol",varlist)
1006            return s
1007        else:
1008            from asap._asap import averagepol_insitu as _avpol
1009            _avpol(self, mask, weight)
1010            self._add_history("average_pol",varlist)
1011            return
1012
1013    def smooth(self, kernel="hanning", width=5.0, insitu=None, allaxes=None):
1014        """
1015        Smooth the spectrum by the specified kernel (conserving flux).
1016        Parameters:
1017            scan:       The input scan
1018            kernel:     The type of smoothing kernel. Select from
1019                        'hanning' (default), 'gaussian' and 'boxcar'.
1020                        The first three characters are sufficient.
1021            width:      The width of the kernel in pixels. For hanning this is
1022                        ignored otherwise it defauls to 5 pixels.
1023                        For 'gaussian' it is the Full Width Half
1024                        Maximum. For 'boxcar' it is the full width.
1025            insitu:     if False a new scantable is returned.
1026                        Otherwise, the scaling is done in-situ
1027                        The default is taken from .asaprc (False)
1028            allaxes:    If True (default) apply to all spectra. Otherwise
1029                        apply only to the selected (beam/pol/if)spectra only
1030                        The default is taken from .asaprc (True if none)
1031        Example:
1032             none
1033        """
1034        if allaxes is None: allaxes = rcParams['scantable.allaxes']
1035        if insitu is None: insitu = rcParams['insitu']
1036        varlist = vars()
1037        if not insitu:
1038            from asap._asap import smooth as _smooth
1039            s = scantable(_smooth(self,kernel,width,allaxes))
1040            s._add_history("smooth", varlist)
1041            return s
1042        else:
1043            from asap._asap import smooth_insitu as _smooth
1044            _smooth(self,kernel,width,allaxes)
1045            self._add_history("smooth", varlist)
1046            return
1047
1048    def poly_baseline(self, mask=None, order=0, insitu=None):
1049        """
1050        Return a scan which has been baselined (all rows) by a polynomial.
1051        Parameters:
1052            scan:    a scantable
1053            mask:    an optional mask
1054            order:   the order of the polynomial (default is 0)
1055            insitu:      if False a new scantable is returned.
1056                         Otherwise, the scaling is done in-situ
1057                         The default is taken from .asaprc (False)
1058        Example:
1059            # return a scan baselined by a third order polynomial,
1060            # not using a mask
1061            bscan = scan.poly_baseline(order=3)
1062        """
1063        if insitu is None: insitu = rcParams['insitu']
1064        varlist = vars()
1065        if mask is None:
1066            from numarray import ones
1067            mask = list(ones(scan.nchan()))
1068        from asap.asapfitter import fitter
1069        f = fitter()
1070        f._verbose(True)
1071        f.set_scan(self, mask)
1072        f.set_function(poly=order)
1073        sf = f.auto_fit(insitu)
1074        if insitu:
1075            self._add_history("poly_baseline", varlist)
1076            return
1077        else:
1078            sf._add_history("poly_baseline", varlist)
1079            return sf
1080
1081    def auto_poly_baseline(self, mask=None, edge=(0,0), order=0,
1082                           threshold=3,insitu=None):
1083        """
1084        Return a scan which has been baselined (all rows) by a polynomial.
1085        Spectral lines are detected first using linefinder and masked out
1086        to avoid them affecting the baseline solution.
1087
1088        Parameters:
1089            scan:    a scantable
1090            mask:       an optional mask retreived from scantable
1091            edge:       an optional number of channel to drop at
1092                        the edge of spectrum. If only one value is
1093                        specified, the same number will be dropped from
1094                        both sides of the spectrum. Default is to keep
1095                        all channels
1096            order:      the order of the polynomial (default is 0)
1097            threshold:  the threshold used by line finder. It is better to
1098                        keep it large as only strong lines affect the
1099                        baseline solution.
1100            insitu:     if False a new scantable is returned.
1101                        Otherwise, the scaling is done in-situ
1102                        The default is taken from .asaprc (False)
1103
1104        Example:
1105            sc2=auto_poly_baseline(order=7)
1106        """
1107        if insitu is None: insitu = rcParams['insitu']
1108        varlist = vars()
1109        from asap.asapfitter import fitter
1110        from asap.asaplinefind import linefinder
1111        from asap import _is_sequence_or_number as _is_valid
1112       
1113        if not _is_valid(edge, int):
1114            print "Parameter 'edge' as to be an integer or a pair of integers"
1115            return
1116       
1117        # setup fitter
1118        f = fitter()
1119        f._verbose(True)
1120        f.set_function(poly=order)
1121
1122        # setup line finder
1123        fl=linefinder()
1124        fl.set_options(threshold=threshold)
1125
1126        if not insitu:
1127            workscan=self.copy()
1128        else:
1129            workscan=self
1130
1131        vb=workscan._vb
1132        # remember the verbose parameter and selection
1133        workscan._vb=False
1134        sel=workscan.get_cursor()
1135        rows=range(workscan.nrow())
1136        for i in range(workscan.nbeam()):
1137            workscan.setbeam(i)
1138            for j in range(workscan.nif()):
1139                workscan.setif(j)
1140                for k in range(workscan.npol()):
1141                    workscan.setpol(k)
1142                    if f._vb:
1143                       print "Processing:"
1144                       print 'Beam[%d], IF[%d], Pol[%d]' % (i,j,k)
1145                    for iRow in rows:
1146                       fl.set_scan(workscan,mask,edge)
1147                       fl.find_lines(iRow)
1148                       f.set_scan(workscan, fl.get_mask())
1149                       f.x=workscan._getabcissa(iRow)
1150                       f.y=workscan._getspectrum(iRow)
1151                       f.data=None
1152                       f.fit()
1153                       x=f.get_parameters()
1154                       workscan._setspectrum(f.fitter.getresidual(),iRow)
1155        workscan.set_cursor(sel[0],sel[1],sel[2])
1156        workscan._vb = vb
1157        if not insitu:
1158           return workscan
1159
1160    def rotate_linpolphase(self, angle, allaxes=None):
1161        """
1162        Rotate the phase of the complex polarization O=Q+iU correlation. 
1163        This is always done in situ in the raw data.  So if you call this
1164        function more than once then each call rotates the phase further.
1165        Parameters:
1166            angle:   The angle (degrees) to rotate (add) by.
1167            allaxes: If True apply to all spectra. Otherwise
1168                     apply only to the selected (beam/pol/if)spectra only.
1169                     The default is taken from .asaprc (True if none)
1170        Examples:
1171            scan.rotate_linpolphase(2.3)
1172    """
1173        if allaxes is None: allaxes = rcParams['scantable.allaxes']
1174        varlist = vars()
1175        from asap._asap import _rotate_linpolphase as _rotate
1176        _rotate(self, angle, allaxes)
1177        self._add_history("rotate_linpolphase", varlist)
1178        return
1179   
1180
1181    def rotate_xyphase(self, angle, allaxes=None):
1182        """
1183        Rotate the phase of the XY correlation.  This is always done in situ
1184        in the data.  So if you call this function more than once
1185        then each call rotates the phase further.
1186        Parameters:
1187            angle:   The angle (degrees) to rotate (add) by.
1188            allaxes: If True apply to all spectra. Otherwise
1189                     apply only to the selected (beam/pol/if)spectra only.
1190                     The default is taken from .asaprc (True if none)
1191        Examples:
1192            scan.rotate_xyphase(2.3)
1193        """
1194        if allaxes is None: allaxes = rcParams['scantable.allaxes']
1195        varlist = vars()
1196        from asap._asap import rotate_xyphase as _rotate_xyphase
1197        _rotate_xyphase(self, angle, allaxes)
1198        self._add_history("rotate_xyphase", varlist)
1199        return
1200
1201
1202    def add(self, offset, insitu=None, allaxes=None):
1203        """
1204        Return a scan where all spectra have the offset added
1205        Parameters:
1206            offset:      the offset
1207            insitu:      if False a new scantable is returned.
1208                         Otherwise, the scaling is done in-situ
1209                         The default is taken from .asaprc (False)
1210            allaxes:     if True apply to all spectra. Otherwise
1211                         apply only to the selected (beam/pol/if)spectra only
1212                         The default is taken from .asaprc (True if none)
1213        """
1214        if allaxes is None: allaxes = rcParams['scantable.allaxes']
1215        if insitu is None: insitu = rcParams['insitu']
1216        varlist = vars()
1217        if not insitu:
1218            from asap._asap import add as _add
1219            s = scantable(_add(self, offset, allaxes))
1220            s._add_history("add",varlist)
1221            return s
1222        else:
1223            from asap._asap import add_insitu as _add
1224            _add(self, offset, allaxes)
1225            self._add_history("add",varlist)
1226            return
1227
1228    def scale(self, factor, insitu=None, allaxes=None, tsys=True):
1229        """
1230        Return a scan where all spectra are scaled by the give 'factor'
1231        Parameters:
1232            factor:      the scaling factor
1233            insitu:      if False a new scantable is returned.
1234                         Otherwise, the scaling is done in-situ
1235                         The default is taken from .asaprc (False)
1236            allaxes:     if True apply to all spectra. Otherwise
1237                         apply only to the selected (beam/pol/if)spectra only.
1238                         The default is taken from .asaprc (True if none)
1239            tsys:        if True (default) then apply the operation to Tsys
1240                         as well as the data
1241        """
1242        if allaxes is None: allaxes = rcParams['scantable.allaxes']
1243        if insitu is None: insitu = rcParams['insitu']
1244        varlist = vars()
1245        if not insitu:
1246            from asap._asap import scale as _scale
1247            s = scantable(_scale(self, factor, allaxes, tsys))
1248            s._add_history("scale",varlist)
1249            return s
1250        else:
1251            from asap._asap import scale_insitu as _scale
1252            _scale(self, factor, allaxes, tsys)
1253            self._add_history("scale",varlist)
1254            return
1255
1256
1257    def quotient(self, other, isreference=True, preserve=True):
1258        """
1259        Return the quotient of a 'source' (signal) scan and a 'reference' scan.
1260        The reference can have just one row, even if the signal has many.
1261        Otherwise they must have the same number of rows.
1262        The cursor of the output scan is set to 0
1263        Parameters:
1264            source:         the 'other' scan
1265            isreference:    if the 'other' scan is the reference (default)
1266                            or source
1267            preserve:       you can preserve (default) the continuum or
1268                            remove it.  The equations used are
1269                            preserve - Output = Tref * (sig/ref) - Tref
1270                            remove   - Output = Tref * (sig/ref) - Tsig
1271        Example:
1272            # src is a scantable for an 'on' scan, ref for an 'off' scantable
1273            q1 = src.quotient(ref)
1274            q2 = ref.quotient(src, isreference=False)
1275            # gives the same result
1276        """
1277        from asap._asap import quotient as _quot
1278        if isreference:
1279            return scantable(_quot(self, other, preserve))
1280        else:
1281            return scantable(_quot(other, self, preserve))
1282
1283    def __add__(self, other):
1284        varlist = vars()
1285        s = None
1286        if isinstance(other, scantable):
1287            from asap._asap import b_operate as _bop           
1288            s = scantable(_bop(self, other, 'add', True))
1289        elif isinstance(other, float):
1290            from asap._asap import add as _add
1291            s = scantable(_add(self, other, True))
1292        else:
1293            print "Other input is not a scantable or float value"
1294            return
1295        s._add_history("operator +", varlist)
1296        return s
1297
1298
1299    def __sub__(self, other):
1300        """
1301        implicit on all axes and on Tsys
1302        """
1303        varlist = vars()
1304        s = None
1305        if isinstance(other, scantable):
1306            from asap._asap import b_operate as _bop           
1307            s = scantable(_bop(self, other, 'sub', True))
1308        elif isinstance(other, float):
1309            from asap._asap import add as _add
1310            s = scantable(_add(self, -other, True))
1311        else:
1312            print "Other input is not a scantable or float value"
1313            return
1314        s._add_history("operator -", varlist)
1315        return s
1316   
1317    def __mul__(self, other):
1318        """
1319        implicit on all axes and on Tsys
1320        """
1321        varlist = vars()
1322        s = None
1323        if isinstance(other, scantable):
1324            from asap._asap import b_operate as _bop           
1325            s = scantable(_bop(self, other, 'mul', True))
1326        elif isinstance(other, float):
1327            if other == 0.0:
1328                print "Multiplying by zero is not recommended."
1329                return
1330            from asap._asap import scale as _sca
1331            s = scantable(_sca(self, other, True))
1332        else:
1333            print "Other input is not a scantable or float value"
1334            return
1335        s._add_history("operator *", varlist)
1336        return s
1337   
1338
1339    def __div__(self, other):
1340        """
1341        implicit on all axes and on Tsys
1342        """
1343        varlist = vars()
1344        s = None
1345        if isinstance(other, scantable):
1346            from asap._asap import b_operate as _bop           
1347            s = scantable(_bop(self, other, 'div', True))
1348        elif isinstance(other, float):
1349            if other == 0.0:
1350                print "Dividing by zero is not recommended"
1351                return
1352            from asap._asap import scale as _sca
1353            s = scantable(_sca(self, 1.0/other, True))
1354        else:
1355            print "Other input is not a scantable or float value"
1356            return
1357        s._add_history("operator /", varlist)
1358        return s
1359
1360    def _add_history(self, funcname, parameters):
1361        # create date
1362        sep = "##"
1363        from datetime import datetime
1364        dstr = datetime.now().strftime('%Y/%m/%d %H:%M:%S')
1365        hist = dstr+sep
1366        hist += funcname+sep#cdate+sep
1367        if parameters.has_key('self'): del parameters['self']
1368        for k,v in parameters.iteritems():
1369            if type(v) is dict:
1370                for k2,v2 in v.iteritems():
1371                    hist += k2
1372                    hist += "="
1373                    if isinstance(v2,scantable):
1374                        hist += 'scantable'
1375                    elif k2 == 'mask':
1376                        if isinstance(v2,list) or isinstance(v2,tuple):
1377                            hist += str(self._zip_mask(v2))
1378                        else:
1379                            hist += str(v2)
1380                    else:
1381                        hist += str(v2)
1382            else:
1383                hist += k
1384                hist += "="
1385                if isinstance(v,scantable):
1386                    hist += 'scantable'
1387                elif k == 'mask':
1388                    if isinstance(v,list) or isinstance(v,tuple):
1389                        hist += str(self._zip_mask(v))
1390                    else:
1391                        hist += str(v)
1392                else:
1393                    hist += str(v)
1394            hist += sep
1395        hist = hist[:-2] # remove trailing '##'
1396        self._addhistory(hist)
1397       
1398
1399    def _zip_mask(self, mask):
1400        mask = list(mask)
1401        i = 0
1402        segments = []
1403        while mask[i:].count(1):
1404            i += mask[i:].index(1)
1405            if mask[i:].count(0):
1406                j = i + mask[i:].index(0)
1407            else:
1408                j = len(mask)               
1409            segments.append([i,j])
1410            i = j       
1411        return segments
Note: See TracBrowser for help on using the repository browser.