source: tags/asap1/python/scantable.py

Last change on this file was 794, checked in by mar637, 18 years ago

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