Changeset 1603 for branches/alma


Ignore:
Timestamp:
07/18/09 06:35:47 (15 years ago)
Author:
TakTsutsumi
Message:

New Development: No, merge with asap2.3.1

JIRA Issue: Yes CAS-1450

Ready to Release: Yes/No?

Interface Changes: Yes/No?

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes

Module(s): single dish

Description: Upgrade of alma branch based on ASAP2.2.0

(rev.1562) to ASAP2.3.1 (rev.1561)


Location:
branches/alma
Files:
30 edited

Legend:

Unmodified
Added
Removed

  • branches/alma/python/__init__.py

    r1494 r1603  
    3232    os.mkdir(userdir)
    3333    #shutil.copyfile(asapdata+"/data/ipythonrc-asap", userdir+"/ipythonrc-asap")
     34    # commented out by TT on 2009.06.23 for casapy use
     35    ##shutil.copyfile(asapdata+"/data/ipy_user_conf.py",
     36    ##                userdir+"/ipy_user_conf.py")
    3437    f = file(userdir+"/asapuserfuncs.py", "w")
    3538    f.close()
    3639    f = file(userdir+"/ipythonrc", "w")
    3740    f.close()
     41# commented out by TT on 2009.06.23 for casapy use
     42##else:
     43    # upgrade to support later ipython versions
     44    ##if not os.path.exists(userdir+"/ipy_user_conf.py"):
     45    ##    shutil.copyfile(asapdata+"/data/ipy_user_conf.py",
     46    ##                    userdir+"/ipy_user_conf.py")
     47
    3848# remove from namespace
    3949del asapdata, userdir, shutil, platform
    … …  
    99109    'plotter.histogram'  : [False, _validate_bool],
    100110    'plotter.papertype'  : ['A4', str],
     111    'plotter.xaxisformatting' : ['asap', str],
    101112
    102113    # scantable
    … …  
    105116    'scantable.freqframe' : ['LSRK', str],  #default frequency frame
    106117    'scantable.verbosesummary'   : [False, _validate_bool],
    107     'scantable.storage'   : ['memory', str]
     118    'scantable.storage'   : ['memory', str],
     119    'scantable.history'   : [True, _validate_bool],
     120    'scantable.reference'      : ['.*(e|w|_R)$', str]
    108121    # fitter
    109122    }
    … …  
    136149plotter.ganged             : True
    137150
    138 # decimate the number of points plotted bya afactor of
     151# decimate the number of points plotted by a factor of
    139152# nchan/1024
    140153plotter.decimate           : False
    … …  
    150163plotter.papertype          : A4
    151164
     165# The formatting style of the xaxis
     166plotter.xaxisformatting    : 'asap' or 'mpl'
     167
    152168# scantable
    153169
    154170# default storage of scantable ('memory'/'disk')
    155171scantable.storage          : memory
     172
     173# write history of each call to scantable
     174scantable.history          : True
     175
    156176# default ouput format when saving
    157177scantable.save             : ASAP
     178
    158179# auto averaging on read
    159180scantable.autoaverage      : True
    … …  
    166187scantable.verbosesummary   : False
    167188
     189# Control the identification of reference (off) scans
     190# This is has to be a regular expression
     191scantable.reference         : .*(e|w|_R)$
    168192# Fitter
    169193"""
    … …  
    242266    for k,v in kwargs.items():
    243267        name = aliases.get(k) or k
    244         key = '%s.%s' % (group, name)
     268        if len(group):
     269            key = '%s.%s' % (group, name)
     270        else:
     271            key = name
    245272        if not rcParams.has_key(key):
    246273            raise KeyError('Unrecognized key "%s" for group "%s" and name "%s"' % (key, group, name))
    … …  
    352379if rcParams['useplotter']:
    353380    try:
    354         from  asapplotter import asapplotter
    355         gui = os.environ.has_key('DISPLAY') and rcParams['plotter.gui']
    356         if gui:
    357             import pylab as xyplotter
    358             plotter = asapplotter(gui)
    359             del gui
     381        from asapplotter import asapplotter
     382        gui = os.environ.has_key('DISPLAY') and rcParams['plotter.gui']
     383        if gui:
     384            import matplotlib
     385            matplotlib.use("TkAgg")
     386        import pylab
     387        xyplotter = pylab
     388        plotter = asapplotter(gui)
     389        del gui
    360390    except ImportError:
    361         print "Matplotlib not installed. No plotting available"
     391        print "Matplotlib not installed. No plotting available"
    362392
    363393__date__ = '$Date$'.split()[1]
    364 __version__  = '2.2.0 alma'
     394__version__  = '2.3.1 alma'
    365395# nrao casapy specific, get revision number
    366396#__revision__ = ' unknown '
    … …  
    369399if os.path.isfile(revinfo):
    370400    f = file(revinfo)
    371     f.readline()
     401    #f.readline()
    372402    revsionno=f.readline()
    373403    f.close()
    … …  
    427457            get_fluxunit    - get the brightness flux unit
    428458            set_fluxunit    - set the brightness flux unit
     459            set_sourcetype  - set the type of the source - source or reference
    429460            create_mask     - return an mask in the current unit
    430461                              for the given region. The specified regions
    … …  
    432463            get_restfreqs   - get the current list of rest frequencies
    433464            set_restfreqs   - set a list of rest frequencies
    434             shift_refpix    - shift the reference pixel of the IFs
     465            shift_refpix    - shift the reference pixel of the IFs
     466            set_spectrum    - overwrite the spectrum for a given row
     467            get_spectrum    - retrieve the spectrum for a given
     468            get_mask        - retrieve the mask for a given
    435469            flag            - flag selected channels in the data
    436470            lag_flag        - flag specified frequency in the data
    … …  
    439473            nbeam,nif,nchan,npol - the number of beams/IFs/Pols/Chans
    440474            nscan           - the number of scans in the scantable
    441             nrow            - te number of spectra in the scantable
     475            nrow            - the number of spectra in the scantable
    442476            history         - print the history of the scantable
    443477            get_fit         - get a fit which has been stored witnh the data

  • branches/alma/python/asapfitter.py

    r1461 r1603  
    33from asap import print_log
    44from asap import _n_bools
     5from asap import mask_and
    56
    67class fitter:
    … …  
    5354        Parameters:
    5455            thescan:     a scantable
    55             mask:        a msk retireved from the scantable
     56            mask:        a msk retrieved from the scantable
    5657        """
    5758        if not thescan:
    … …  
    142143                self.x = self.data._getabcissa(row)
    143144                self.y = self.data._getspectrum(row)
     145                self.mask = mask_and(self.mask, self.data._getmask(row))
    144146                from asap import asaplog
    145147                asaplog.push("Fitting:")
    146148                i = row
    147                 out = "Scan[%d] Beam[%d] IF[%d] Pol[%d] Cycle[%d]" % (self.data.getscan(i),self.data.getbeam(i),self.data.getif(i),self.data.getpol(i), self.data.getcycle(i))
     149                out = "Scan[%d] Beam[%d] IF[%d] Pol[%d] Cycle[%d]" % (self.data.getscan(i),
     150                                                                      self.data.getbeam(i),
     151                                                                      self.data.getif(i),
     152                                                                      self.data.getpol(i),
     153                                                                      self.data.getcycle(i))
    148154                asaplog.push(out,False)
    149155        self.fitter.setdata(self.x, self.y, self.mask)
    … …  
    245251
    246252    def set_gauss_parameters(self, peak, centre, fwhm,
    247                              peakfixed=0, centerfixed=0,
     253                             peakfixed=0, centrefixed=0,
    248254                             fwhmfixed=0,
    249255                             component=0):
    … …  
    253259            peak, centre, fwhm:  The gaussian parameters
    254260            peakfixed,
    255             centerfixed,
     261            centrefixed,
    256262            fwhmfixed:           Optional parameters to indicate if
    257263                                 the paramters should be held fixed during
    … …  
    270276        if 0 <= component < len(self.components):
    271277            d = {'params':[peak, centre, fwhm],
    272                  'fixed':[peakfixed, centerfixed, fwhmfixed]}
     278                 'fixed':[peakfixed, centrefixed, fwhmfixed]}
    273279            self.set_parameters(d, component)
    274280        else:
    … …  
    604610        asaplog.push("Fitting:")
    605611        for r in rows:
    606             out = " Scan[%d] Beam[%d] IF[%d] Pol[%d] Cycle[%d]" % (scan.getscan(r),scan.getbeam(r),scan.getif(r),scan.getpol(r), scan.getcycle(r))
     612            out = " Scan[%d] Beam[%d] IF[%d] Pol[%d] Cycle[%d]" % (scan.getscan(r),
     613                                                                   scan.getbeam(r),
     614                                                                   scan.getif(r),
     615                                                                   scan.getpol(r),
     616                                                                   scan.getcycle(r))
    607617            asaplog.push(out, False)
    608618            self.x = scan._getabcissa(r)
    609619            self.y = scan._getspectrum(r)
     620            self.mask = mask_and(self.mask, scan._getmask(r))
    610621            self.data = None
    611622            self.fit()
     623            x = self.get_parameters()
    612624            fpar = self.get_parameters()
    613625            if plot:

  • branches/alma/python/asaplotbase.py

    r1446 r1603  
    1515from matplotlib.ticker import OldScalarFormatter
    1616from matplotlib.ticker import NullLocator
    17 from matplotlib.transforms import blend_xy_sep_transform
     17
     18# API change in mpl >= 0.98
     19try:
     20    from matplotlib.transforms import blended_transform_factory
     21except ImportError:
     22    from matplotlib.transforms import blend_xy_sep_transform  as blended_transform_factory
    1823
    1924if int(matplotlib.__version__.split(".")[1]) < 87:
    … …  
    161166        y2 = range(l2)
    162167        m2 = range(l2)
    163         #ymsk = y.raw_mask()
    164         #ydat = y.raw_data()
    165         ymsk = y.mask
    166         ydat = y.data
     168        ymsk = None
     169        ydat = None
     170        if hasattr(y, "raw_mask"):
     171            # numpy < 1.1
     172            ymsk = y.raw_mask()
     173            ydat = y.raw_data()
     174        else:
     175            ymsk = y.mask
     176            ydat = y.data
    167177        for i in range(l2):
    168178            x2[i] = x[i/2]
    … …  
    410420                if fname[-3:].lower() == ".ps":
    411421                    from matplotlib import __version__ as mv
    412                     w = self.figure.figwidth.get()
    413                     h = self.figure.figheight.get()
     422                    w = self.figure.get_figwidth()
     423                    h = self.figure.get_figheight()
    414424
    415425                    if orientation is None:
    … …  
    428438                    ow = ds * w
    429439                    oh = ds * h
    430                     self.figure.set_figsize_inches((ow, oh))
     440                    self.figure.set_size_inches((ow, oh))
    431441                    self.figure.savefig(fname, orientation=orientation,
    432442                                        papertype=papertype.lower())
    433                     self.figure.set_figsize_inches((w, h))
     443                    self.figure.set_size_inches((w, h))
    434444                    print 'Written file %s' % (fname)
    435445                else:
    … …  
    617627                    self.subplots[i]['axes'] = self.figure.add_subplot(rows,
    618628                                                cols, i+1)
    619                     self.subplots[i]['axes'].xaxis.set_major_formatter(OldScalarFormatter())
     629                    if asaprcParams['plotter.xaxisformatting'] == 'mpl':
     630                        self.subplots[i]['axes'].xaxis.set_major_formatter(OldScalarFormatter())
    620631                else:
    621632                    if i == 0:
    622633                        self.subplots[i]['axes'] = self.figure.add_subplot(rows,
    623634                                                cols, i+1)
    624                         self.subplots[i]['axes'].xaxis.set_major_formatter(OldScalarFormatter())
     635                        if asaprcParams['plotter.xaxisformatting'] != 'mpl':
     636                           
     637                            self.subplots[i]['axes'].xaxis.set_major_formatter(OldScalarFormatter())
    625638                    else:
    626639                        self.subplots[i]['axes'] = self.figure.add_subplot(rows,
    … …  
    709722            for sp in self.subplots:
    710723                ax = sp['axes']
    711                 s = rcParams['axes.titlesize']
    712                 tsize = s-(self.cols+self.rows-1)
    713                 ax.title.set_size(max(tsize,9))
     724                s = ax.title.get_size()
     725                tsize = s-(self.cols+self.rows)
     726                ax.title.set_size(tsize)
    714727                fp = FP(size=rcParams['axes.labelsize'])
    715728                setp(ax.get_xticklabels(), fontsize=xts)
    … …  
    770783        if rotate > 0.0: lbloffset = 0.03*len(label)
    771784        peakoffset = 0.01
    772         xy0 = ax.transData.xy_tup((x,y))
    773         # get relative coords
    774         xy = ax.transAxes.inverse_xy_tup(xy0)
     785        xy = None
     786        xy0 = None
     787        # matplotlib api change 0.98 is using transform now
     788        if hasattr(ax.transData, "inverse_xy_tup"):
     789            # get relative coords
     790            xy0 = ax.transData.xy_tup((x,y))
     791            xy = ax.transAxes.inverse_xy_tup(xy0)
     792        else:
     793            xy0 = ax.transData.transform((x,y))
     794            # get relative coords
     795            xy = ax.transAxes.inverted().transform(xy0)
    775796        if location.lower() == 'top':
    776797            ymax = 1.0-lbloffset
    … …  
    783804            valign = 'top'
    784805            ylbl = ymin-0.01
    785         trans = blend_xy_sep_transform(ax.transData, ax.transAxes)
     806        trans = blended_transform_factory(ax.transData, ax.transAxes)
    786807        l = ax.axvline(x, ymin, ymax, color='black', **kwargs)
    787808        t = ax.text(x, ylbl ,label, verticalalignment=valign,

  • branches/alma/python/asaplotgui.py

    r1153 r1603  
    55from asap.asaplotbase import *
    66import Tkinter as Tk
     7import matplotlib
    78from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg, \
    89        FigureManagerTkAgg
    910# Force use of the newfangled toolbar.
    10 import matplotlib
    11 matplotlib.use("TkAgg")
    1211matplotlib.rcParams['toolbar'] = 'toolbar2'
    1312

  • branches/alma/python/linecatalog.py

    r1259 r1603  
    11"""
    2 A representation of a spectra line catalog.
     2A representation of a spectral line catalog.
    33
    44Author: Malte Marquarding
    … …  
    3434            else:
    3535                raise IOError(msg)
     36
     37    def __repr__(self):
     38        return lcbase.summary(self, -1)
    3639
    3740    def summary(self):

  • branches/alma/python/scantable.py

    r1530 r1603  
    4545            Scantable.__init__(self, filename)
    4646        else:
    47             if isinstance(filename, str):
     47            if isinstance(filename, str):# or \
     48#                (isinstance(filename, list) or isinstance(filename, tuple)) \
     49#                  and isinstance(filename[-1], str):
    4850                import os.path
    4951                filename = os.path.expandvars(filename)
    … …  
    9698                                       'MS2' (saves as an aips++
    9799                                              MeasurementSet V2)
     100                                       'FITS' (save as image FITS - not
     101                                               readable by class)
     102                                       'CLASS' (save as FITS readable by CLASS)
    98103            overwrite:   If the file should be overwritten if it exists.
    99104                         The default False is to return with warning
    … …  
    279284            return info
    280285
     286    def get_spectrum(self, rowno):
     287        """Return the spectrum for the current row in the scantable as a list.
     288        Parameters:
     289             rowno:   the row number to retrieve the spectrum from       
     290        """
     291        return self._getspectrum(rowno)
     292
     293    def get_mask(self, rowno):
     294        """Return the mask for the current row in the scantable as a list.
     295        Parameters:
     296             rowno:   the row number to retrieve the mask from       
     297        """
     298        return self._getmask(rowno)
     299
     300    def set_spectrum(self, spec, rowno):
     301        """Return the spectrum for the current row in the scantable as a list.
     302        Parameters:
     303             spec:   the spectrum
     304             rowno:    the row number to set the spectrum for       
     305        """
     306        assert(len(spec) == self.nchan())
     307        return self._setspectrum(spec, rowno)
    281308
    282309    def get_selection(self):
    … …  
    518545        times = self._get_column(self._gettime, row)
    519546        if not asdatetime:
    520             return times 
     547            return times
    521548        format = "%Y/%m/%d/%H:%M:%S"
    522549        if isinstance(times, list):
    … …  
    755782                          to remove
    756783            width:        the width of the frequency to remove, to remove a
    757                           range of frequencies aroung the centre.
     784                          range of frequencies around the centre.
    758785            unit:         the frequency unit (default "GHz")
    759786        Notes:
    … …  
    964991        IF 1 gets restfreq 2e9.
    965992        ********NEED TO BE UPDATED end************
    966         You can also specify the frequencies via a linecatalog/
     993        You can also specify the frequencies via a linecatalog.
    967994
    968995        Parameters:
    … …  
    16871714            return s
    16881715
     1716    def set_sourcetype(self, match, matchtype="pattern",
     1717                       sourcetype="reference"):
     1718        """
     1719        Set the type of the source to be an source or reference scan
     1720        using the provided pattern:
     1721        Parameters:
     1722            match:          a Unix style pattern, regular expression or selector
     1723            matchtype:      'pattern' (default) UNIX style pattern or
     1724                            'regex' regular expression
     1725            sourcetype:     the type of the source to use (source/reference)
     1726        """
     1727        varlist = vars()
     1728        basesel = self.get_selection()
     1729        stype = -1
     1730        if sourcetype.lower().startswith("r"):
     1731            stype = 1
     1732        elif sourcetype.lower().startswith("s"):
     1733            stype = 0
     1734        else:
     1735            raise ValueError("Illegal sourcetype use s(ource) or r(eference)")
     1736        if matchtype.lower().startswith("p"):
     1737            matchtype = "pattern"
     1738        elif matchtype.lower().startswith("r"):
     1739            matchtype = "regex"
     1740        else:
     1741            raise ValueError("Illegal matchtype, use p(attern) or r(egex)")
     1742        sel = selector()
     1743        if isinstance(match, selector):
     1744            sel = match
     1745        else:
     1746            sel.set_query("SRCNAME == %s('%s')" % (matchtype, match))
     1747        self.set_selection(basesel+sel)
     1748        self._setsourcetype(stype)
     1749        self.set_selection(basesel)
     1750        s._add_history("set_sourcetype", varlist)
     1751
    16891752    def auto_quotient(self, preserve=True, mode='paired'):
    16901753        """
    … …  
    17031766                            '_e'/'_w' (Tid) and matches
    17041767                            on/off pairs from the observing pattern
    1705                 'time'
    1706                    finds the closest off in time
     1768                            'time'
     1769                            finds the closest off in time
    17071770
    17081771        """
    … …  
    18641927            return fit.as_dict()
    18651928
     1929    def flag_nans(self):
     1930        """
     1931        Utility function to flag NaN values in the scantable.
     1932        """
     1933        import numpy
     1934        basesel = self.get_selection()
     1935        for i in range(self.nrow()):
     1936            sel = selector()+basesel
     1937            sel.set_scans(self.getscan(i))
     1938            sel.set_beams(self.getbeam(i))
     1939            sel.set_ifs(self.getif(i))
     1940            sel.set_polarisations(self.getpol(i))
     1941            self.set_selection(sel)
     1942            nans = numpy.isnan(self._getspectrum(0))
     1943        if numpy.any(nans):
     1944            bnans = [ bool(v) for v in nans]
     1945            self.flag(bnans)
     1946        self.set_selection(basesel)
     1947       
     1948
    18661949    def _add_history(self, funcname, parameters):
     1950        if not rcParams['scantable.history']:
     1951            return
    18671952        # create date
    18681953        sep = "##"
    … …  
    19542039            tbl = Scantable(stype)
    19552040            r = stfiller(tbl)
     2041            rx = rcParams['scantable.reference']
     2042            r._setreferenceexpr(rx)
    19562043            msg = "Importing %s..." % (name)
    19572044            asaplog.push(msg, False)
    … …  
    19592046            r._open(name, -1, -1, getpt)
    19602047            r._read()
    1961             #tbl = r._getdata()
    19622048            if average:
    19632049                tbl = self._math._average((tbl, ), (), 'NONE', 'SCAN')
    1964                 #tbl = tbl2
    19652050            if not first:
    19662051                tbl = self._math._merge([self, tbl])
    1967                 #tbl = tbl2
    19682052            Scantable.__init__(self, tbl)
    19692053            r._close()
    … …  
    19742058        #self.set_freqframe(rcParams['scantable.freqframe'])
    19752059
     2060    def __getitem__(self, key):
     2061        if key < 0:
     2062            key += self.nrow()
     2063        if key >= self.nrow():
     2064            raise IndexError("Row index out of range.")
     2065        return self._getspectrum(key)
     2066
     2067    def __setitem__(self, key, value):
     2068        if key < 0:
     2069            key += self.nrow()
     2070        if key >= self.nrow():
     2071            raise IndexError("Row index out of range.")
     2072        if not hasattr(value, "__len__") or \
     2073                len(value) > self.nchan(self.getif(key)):
     2074            raise ValueError("Spectrum length doesn't match.")
     2075        return self._setspectrum(value, key)
     2076
     2077    def __len__(self):
     2078        return self.nrow()
     2079
     2080    def __iter__(self):
     2081        for i in range(len(self)):
     2082            yield self[i]

  • branches/alma/python/selector.py

    r1349 r1603  
    4747   
    4848    # for the americans
    49     set_polarization = set_polarisations
     49    set_polarizations = set_polarisations
     50    # for the lazy
     51    set_pols = set_polarisations
    5052
    5153    def set_ifs(self, ifs=[]):
    … …  
    163165        prefix = "SELECT FROM $1 WHERE "
    164166        return self._gettaql().replace(prefix, "")
     167
    165168    def get_name(self):
    166169        print "NYI"

  • branches/alma/src/LineCatalog.cpp

    r1259 r1603  
    6262void LineCatalog::setPattern(const std::string& name, const std::string& stype)
    6363{
    64   //std::string mode = stype+"('"+name+"')";
    65   std::string taql = "SELECT FROM $1 WHERE Column1 == pattern('"+name+"')";
    66   cerr << taql << endl;
     64  std::string mode = stype+"('"+name+"')";
     65  std::string taql = "SELECT FROM $1 WHERE Column1 == " + mode;
    6766  Table tmp = tableCommand(taql, table_);
    6867  if (tmp.nrow() > 0) table_ = tmp.sort("Column2");
    … …  
    110109      }
    111110  }
    112   /// @todo implement me
    113111  return String(oss);
    114112}
    115113
    116 /*!
    117     \fn asap::LineCatalog::getName(int row)
    118  */
     114
    119115std::string LineCatalog::getName(uint row) const
    120116{
    … …  
    123119}
    124120
    125 double asap::LineCatalog::getFrequency(uint row) const
     121double LineCatalog::getFrequency(uint row) const
    126122{
    127   ROScalarColumn<Double> col(table_, "Column2");
    128   return col(row);
     123  return getDouble("Column2", row);
    129124}
    130125
    131 double asap::LineCatalog::getStrength(uint row) const
     126double LineCatalog::getStrength(uint row) const
    132127{
    133   ROScalarColumn<Double> col(table_, "Column3");
    134   return col(row);
     128  return getDouble("Column4", row);
    135129}
    136130
     131double LineCatalog::getDouble(const std::string& colname, uint row) const {
     132  DataType dtype = table_.tableDesc().columnDesc(colname).dataType();
     133  if (dtype == TpDouble) {
     134    ROScalarColumn<Double> col(table_, colname);
     135    return col(row);
     136  } else if (dtype == TpInt) {
     137    ROScalarColumn<Int> col(table_, colname);
     138    return Double(col(row));
     139  } else {
     140    throw AipsError("Column " + colname + "doesn't contain numerical values." );
     141  }
     142}
    137143
    138144} // namespace

  • branches/alma/src/LineCatalog.h

    r1353 r1603  
    111111  casa::Table setLimits(double lmin, double lmax, const std::string& colname);
    112112
     113  double getDouble(const std::string& colname, uint row) const;
     114
    113115  // the table with seelection
    114116  casa::Table table_;

  • branches/alma/src/Makefile

    r1520 r1603  
    4343RPFITSLIB := -lrpfits
    4444
     45#wcs
     46WCSLIB := -lwcs
     47
    4548G2CROOT := /usr
    46 G2CARCH := $(G2CROOT)/lib/gcc/i386-apple-darwin8.7.1/4.2.0/libgcc.a
    47 G2CARCH := $(G2CROOT)/lib/gcc/powerpc-apple-darwin8.7.0/4.2.0/libgcc.a
    48 G2CARCH := $(G2CROOT)/lib/gcc/i386-redhat-linux/4.1.0/libgcc.a
     49#G2CARCH := $(G2CROOT)/lib/gcc/i386-apple-darwin8.7.1/4.2.0/libgcc.a
     50#G2CARCH := $(G2CROOT)/lib/gcc/powerpc-apple-darwin8.7.0/4.2.0/libgcc.a
     51G2CARCH := $(G2CROOT)/lib/gcc/i386-redhat-linux/4.1.2/libgcc.a
    4952#G2CLIB := $(G2CROOT)/lib/libgfortran.a
    50 G2CLIB := -lg2c
     53G2CLIB := $(G2CARCH)
     54#G2CLIB := -lg2c
    5155
    5256# This assumes all casa libs are static only (*.a)
    … …  
    5660             -lcasa_lattices -lcasa_fits -lcasa_measures -lcasa_measures_f \
    5761             -lcasa_tables -lcasa_scimath -lcasa_scimath_f -lcasa_casa  \
    58              $(CASALIB)/libwcs.a \
     62             $(WCSLIB) \
    5963             $(RPFITSLIB) $(CFITSIOLIB) $(G2CLIB) -lstdc++
    6064
    … …  
    127131             STWriter.o \
    128132             STAsciiWriter.o \
     133             STFITSImageWriter.o \
    129134             Scantable.o \
    130135             Templates.o
    … …  
    172177             STPolLinear.h \
    173178             STWriter.h \
    174              STAsciiWriter.h
     179             STAsciiWriter.h \
     180             STFITSImageWriter.h
    175181
    176182STATICCCLIB := libasap.a

  • branches/alma/src/MathUtils.cpp

    r1530 r1603  
    3535#include <casa/Arrays/MaskedArray.h>
    3636#include <casa/Arrays/MaskArrMath.h>
     37#include <casa/Arrays/VectorSTLIterator.h>
    3738#include <casa/BasicSL/String.h>
    3839#include <scimath/Mathematics/MedianSlider.h>
    … …  
    106107{
    107108  std::vector<std::string> out;
    108   Vector<String>::const_iterator it = in.begin();
    109   for (uInt i=0; it != in.end(); ++it,++i) {
     109  out.reserve(in.nelements());
     110  for (Array<String>::const_iterator it = in.begin(); it != in.end(); ++it) {
    110111    out.push_back(*it);
    111112  }
    … …  
    116117{
    117118  Vector<String> out(in.size());
    118   uInt i=0;
    119   std::vector<std::string>::const_iterator it;
    120   for (it=in.begin();it != in.end();++it,++i) {
    121     out[i] = casa::String(*it);
     119  Array<String>::iterator oit = out.begin();
     120  for (std::vector<std::string>::const_iterator it=in.begin() ;
     121       it != in.end(); ++it,++oit) {
     122    *oit = *it;
    122123  }
    123124  return out;

  • branches/alma/src/MathUtils.h

    r1514 r1603  
    8787
    8888/**
    89  * Convert a std::vector of std::string
    90  * to a casa::Vector casa::String
    91  * @param in
    92  * @return
     89 * Convert casa implementations to stl
     90 * @param in casa string
     91 * @return a std vector of std strings
    9392 */
    9493std::vector<std::string> tovectorstring(const casa::Vector<casa::String>& in);
    9594
    9695/**
    97  * Convert a casa::Vector of casa::String
    98  * to a stl std::vector of stl std::string
     96 * convert stl implementations to casa versions
    9997 * @param in
    10098 * @return

  • branches/alma/src/RowAccumulator.cpp

    r1446 r1603  
    4444    Vector<Bool> dummymsk(m.nelements(), True);
    4545    spectrum_.setData(dummy, dummymsk);
    46     n_.setData(Vector<Float>(v.nelements(), 0.0), dummymsk);
     46    n_.setData(Vector<uInt>(v.nelements(), 0), dummymsk);
    4747    weightSum_.setData(Vector<Float>(v.nelements(), 0.0), dummymsk);
    4848    tsysSum_.resize(tsys.nelements()); tsysSum_=0.0;
    … …  
    5050  // add spectrum related weights, so far it is variance only.
    5151  Float totalweight = 1.0;
    52   totalweight *= addTsys(tsys);
     52
    5353  // only add these if not everything masked
    5454  if ( !allEQ(m, False) ) {
     55    totalweight *= addTsys(tsys);
    5556    totalweight *= addInterval(interval);
    5657    addTime(time);
    … …  
    7980  weightSum_ += wadd;
    8081  spectrum_ += data;
    81   const MaskedArray<Float> inc(Vector<Float>(m.nelements(),1.0), m);
     82  const MaskedArray<uInt> inc(Vector<uInt>(m.nelements(),1), m);
    8283  n_ += inc;
    8384}
    … …  
    125126casa::Double asap::RowAccumulator::getTime( ) const
    126127{
    127   Float n = max(n_);
    128   if (n < 1.0) n = 1.0;
    129   return timeSum_/n;
     128  uInt n = max(n_);
     129  return timeSum_/Float(n);
    130130}
    131131
    … …  
    138138{
    139139  // Return the "total" mask - False where no points have been accumulated.
    140   return (n_.getArray() > Float(0.0));
     140  return (n_.getArray() > uInt(0));
    141141}
    142142
    … …  
    144144{
    145145  // @fixme this assumes tsys.nelements() == 1
    146   return tsysSum_/max(n_);
     146  return tsysSum_/Float(max(n_));
    147147}
    148148
    … …  
    158158  return initialized_;
    159159}
     160

  • branches/alma/src/RowAccumulator.h

    r1446 r1603  
    103103  //these are Vectors
    104104  casa::MaskedArray<casa::Float> spectrum_;
    105   casa::MaskedArray<casa::Float> n_, weightSum_;
     105  casa::MaskedArray<casa::Float> weightSum_;
     106  casa::MaskedArray<casa::uInt> n_;
    106107
    107108  casa::Vector<casa::Bool> userMask_;

  • branches/alma/src/STFiller.cpp

    r1533 r1603  
    55//
    66//
    7 // Author: Malte Marquarding <asap@atnf.csiro.au>, (C) 2006
     7// Author: Malte Marquarding <asap@atnf.csiro.au>, (C) 2006-2007
    88//
    99// Copyright: See COPYING file that comes with this distribution
    … …  
    2828#include <measures/Measures/MeasConvert.h>
    2929
     30#include <atnf/PKSIO/PKSrecord.h>
    3031#include <atnf/PKSIO/PKSreader.h>
     32#ifdef HAS_ALMA
     33 #include <casa/System/ProgressMeter.h>
     34#endif
    3135#include <casa/System/ProgressMeter.h>
    3236#include <atnf/PKSIO/NROReader.h>
    3337
    3438#include <time.h>
     39
    3540
    3641#include "STDefs.h"
    … …  
    4853  header_(0),
    4954  table_(0),
     55  refRx_(".*(e|w|_R)$"),
    5056  nreader_(0)
    5157{
    … …  
    5662  header_(0),
    5763  table_(stbl),
     64  refRx_(".*(e|w|_R)$"),
    5865  nreader_(0)
    5966{
    … …  
    6471  header_(0),
    6572  table_(0),
     73  refRx_(".*(e|w|_R)$"),
    6674  nreader_(0)
    6775{
    … …  
    141149  Int status = reader_->getHeader(header_->observer, header_->project,
    142150                                  header_->antennaname, header_->antennaposition,
    143                                   header_->obstype,header_->equinox,
     151                                  header_->obstype,
     152                                  header_->fluxunit,
     153                                  header_->equinox,
    144154                                  header_->freqref,
    145155                                  header_->utc, header_->reffreq,
    146                                   header_->bandwidth,
    147                                   header_->fluxunit);
     156                                  header_->bandwidth);
    148157
    149158  if (status) {
    … …  
    166175  Bool throwIt = False;
    167176  Instrument inst = STAttr::convertInstrument(header_->antennaname, throwIt);
    168   //header_->fluxunit = "Jy";
     177
    169178  if (inst==ATMOPRA || inst==TIDBINBILLA) {
    170      header_->fluxunit = "K";
     179    header_->fluxunit = "K";
     180  } else {
     181    // downcase for use with Quanta
     182    if (header_->fluxunit == "JY") {
     183      header_->fluxunit = "Jy";
     184    }
    171185  }
    172186  STAttr stattr;
    … …  
    275289  //
    276290
     291/**
    277292  Int    beamNo, IFno, refBeam, scanNo, cycleNo;
    278293  Float  azimuth, elevation, focusAxi, focusRot, focusTan,
    … …  
    287302  Complex         xCalFctr;
    288303  Vector<Complex> xPol;
     304**/
     305
    289306  Double min = 0.0;
    290307  Double max = nInDataRow;
     308#ifdef HAS_ALMA
    291309  ProgressMeter fillpm(min, max, "Data importing progress");
     310#endif
     311  PKSrecord pksrec;
    292312  int n = 0;
    293313  while ( status == 0 ) {
    294     status = reader_->read(scanNo, cycleNo, mjd, interval, fieldName,
    295                           srcName, srcDir, srcPM, srcVel, obsType, IFno,
    296                           refFreq, bandwidth, freqInc, restFreq, tcal, tcalTime,
    297                           azimuth, elevation, parAngle, focusAxi,
    298                           focusTan, focusRot, temperature, pressure,
    299                           humidity, windSpeed, windAz, refBeam,
    300                           beamNo, direction, scanRate,
    301                           tsys, sigma, calFctr, baseLin, baseSub,
    302                           spectra, flagtra, xCalFctr, xPol);
     314    status = reader_->read(pksrec);
    303315    if ( status != 0 ) break;
    304316    n += 1;
    305    
     317
    306318    Regex filterrx(".*[SL|PA]$");
    307319    Regex obsrx("^AT.+");
    308     if ( header_->antennaname.matches(obsrx) && obsType.matches(filterrx)) {
     320    if ( header_->antennaname.matches(obsrx) &&
     321         pksrec.obsType.matches(filterrx)) {
    309322        //cerr << "ignoring paddle scan" << endl;
    310323        continue;
    … …  
    314327    // fields that don't get used and are just passed through asap
    315328    RecordFieldPtr<Array<Double> > srateCol(rec, "SCANRATE");
    316     *srateCol = scanRate;
     329    // MRC changed type from double to float
     330    Vector<Double> sratedbl(pksrec.scanRate.nelements());
     331    convertArray(sratedbl, pksrec.scanRate);
     332    *srateCol = sratedbl;
    317333    RecordFieldPtr<Array<Double> > spmCol(rec, "SRCPROPERMOTION");
    318     *spmCol = srcPM;
     334    *spmCol = pksrec.srcPM;
    319335    RecordFieldPtr<Array<Double> > sdirCol(rec, "SRCDIRECTION");
    320     *sdirCol = srcDir;
     336    *sdirCol = pksrec.srcDir;
    321337    RecordFieldPtr<Double> svelCol(rec, "SRCVELOCITY");
    322     *svelCol = srcVel;
     338    *svelCol = pksrec.srcVel;
    323339    // the real stuff
    324340    RecordFieldPtr<Int> fitCol(rec, "FIT_ID");
    325341    *fitCol = -1;
    326342    RecordFieldPtr<uInt> scanoCol(rec, "SCANNO");
    327     *scanoCol = scanNo-1;
     343    *scanoCol = pksrec.scanNo-1;
    328344    RecordFieldPtr<uInt> cyclenoCol(rec, "CYCLENO");
    329     *cyclenoCol = cycleNo-1;
     345    *cyclenoCol = pksrec.cycleNo-1;
    330346    RecordFieldPtr<Double> mjdCol(rec, "TIME");
    331     *mjdCol = mjd;
     347    *mjdCol = pksrec.mjd;
    332348    RecordFieldPtr<Double> intCol(rec, "INTERVAL");
    333     *intCol = interval;
     349    *intCol = pksrec.interval;
    334350    RecordFieldPtr<String> srcnCol(rec, "SRCNAME");
    335351    RecordFieldPtr<Int> srctCol(rec, "SRCTYPE");
    336352    RecordFieldPtr<String> fieldnCol(rec, "FIELDNAME");
    337     *fieldnCol = fieldName;
     353    *fieldnCol = pksrec.fieldName;
    338354    // try to auto-identify if it is on or off.
    339     Regex rx(".*[e|w|_R]$");
     355    Regex rx(refRx_);
    340356    Regex rx2("_S$");
    341     Int match = srcName.matches(rx);
     357    Int match = pksrec.srcName.matches(rx);
    342358    if (match) {
    343       *srcnCol = srcName;
     359      *srcnCol = pksrec.srcName;
    344360    } else {
    345       *srcnCol = srcName.before(rx2);
    346     }
    347     //*srcnCol = srcName;//.before(rx2);
     361      *srcnCol = pksrec.srcName.before(rx2);
     362    }
     363    //*srcnCol = pksrec.srcName;//.before(rx2);
    348364    *srctCol = match;
    349365    RecordFieldPtr<uInt> beamCol(rec, "BEAMNO");
    350     *beamCol = beamNo-beamOffset_-1;
     366    *beamCol = pksrec.beamNo-beamOffset_-1;
    351367    RecordFieldPtr<Int> rbCol(rec, "REFBEAMNO");
    352368    Int rb = -1;
    353     if (nBeam_ > 1 ) rb = refBeam-1;
     369    if (nBeam_ > 1 ) rb = pksrec.refBeam-1;
    354370    *rbCol = rb;
    355371    RecordFieldPtr<uInt> ifCol(rec, "IFNO");
    356     *ifCol = IFno-ifOffset_- 1;
     372    *ifCol = pksrec.IFno-ifOffset_- 1;
    357373    uInt id;
    358374    /// @todo this has to change when nchan isn't global anymore
    359375    id = table_->frequencies().addEntry(Double(header_->nchan/2),
    360                                             refFreq, freqInc);
     376                                            pksrec.refFreq, pksrec.freqInc);
    361377    RecordFieldPtr<uInt> mfreqidCol(rec, "FREQ_ID");
    362378    *mfreqidCol = id;
    363379
    364     id = table_->molecules().addEntry(restFreq);
     380    id = table_->molecules().addEntry(pksrec.restFreq);
    365381    RecordFieldPtr<uInt> molidCol(rec, "MOLECULE_ID");
    366382    *molidCol = id;
    367383
    368     id = table_->tcal().addEntry(tcalTime, tcal);
     384    id = table_->tcal().addEntry(pksrec.tcalTime, pksrec.tcal);
    369385    RecordFieldPtr<uInt> mcalidCol(rec, "TCAL_ID");
    370386    *mcalidCol = id;
    371     id = table_->weather().addEntry(temperature, pressure, humidity,
    372                                     windSpeed, windAz);
     387    id = table_->weather().addEntry(pksrec.temperature, pksrec.pressure,
     388                                    pksrec.humidity, pksrec.windSpeed,
     389                                    pksrec.windAz);
    373390    RecordFieldPtr<uInt> mweatheridCol(rec, "WEATHER_ID");
    374391    *mweatheridCol = id;
    375392    RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID");
    376     id = table_->focus().addEntry(focusAxi, focusTan, focusRot);
     393    id = table_->focus().addEntry(pksrec.focusAxi, pksrec.focusTan,
     394                                  pksrec.focusRot);
    377395    *mfocusidCol = id;
    378396    RecordFieldPtr<Array<Double> > dirCol(rec, "DIRECTION");
    379     *dirCol = direction;
     397    *dirCol = pksrec.direction;
    380398    RecordFieldPtr<Float> azCol(rec, "AZIMUTH");
    381     *azCol = azimuth;
     399    *azCol = pksrec.azimuth;
    382400    RecordFieldPtr<Float> elCol(rec, "ELEVATION");
    383     *elCol = elevation;
     401    *elCol = pksrec.elevation;
    384402
    385403    RecordFieldPtr<Float> parCol(rec, "PARANGLE");
    386     *parCol = parAngle;
     404    *parCol = pksrec.parAngle;
    387405
    388406    RecordFieldPtr< Array<Float> > specCol(rec, "SPECTRA");
    … …  
    395413    Vector<Float> tsysvec(1);
    396414    // Why is spectra.ncolumn() == 3 for haveXPol_ == True
    397     uInt npol = (spectra.ncolumn()==1 ? 1: 2);
     415    uInt npol = (pksrec.spectra.ncolumn()==1 ? 1: 2);
    398416    for ( uInt i=0; i< npol; ++i ) {
    399       tsysvec = tsys(i);
     417      tsysvec = pksrec.tsys(i);
    400418      *tsysCol = tsysvec;
    401419      *polnoCol = i;
    402420
    403       *specCol = spectra.column(i);
    404       *flagCol = flagtra.column(i);
     421      *specCol = pksrec.spectra.column(i);
     422      *flagCol = pksrec.flagged.column(i);
    405423      table_->table().addRow();
    406424      row.put(table_->table().nrow()-1, rec);
    … …  
    408426    if ( haveXPol_[0] ) {
    409427      // no tsys given for xpol, so emulate it
    410       tsysvec = sqrt(tsys[0]*tsys[1]);
     428      tsysvec = sqrt(pksrec.tsys[0]*pksrec.tsys[1]);
    411429      *tsysCol = tsysvec;
    412430      // add real part of cross pol
    413431      *polnoCol = 2;
    414       Vector<Float> r(real(xPol));
     432      Vector<Float> r(real(pksrec.xPol));
    415433      *specCol = r;
    416434      // make up flags from linears
    417435      /// @fixme this has to be a bitwise or of both pols
    418       *flagCol = flagtra.column(0);// | flagtra.column(1);
     436      *flagCol = pksrec.flagged.column(0);// | pksrec.flagged.column(1);
    419437      table_->table().addRow();
    420438      row.put(table_->table().nrow()-1, rec);
    421439      // ad imaginary part of cross pol
    422440      *polnoCol = 3;
    423       Vector<Float> im(imag(xPol));
     441      Vector<Float> im(imag(pksrec.xPol));
    424442      *specCol = im;
    425443      table_->table().addRow();
    426444      row.put(table_->table().nrow()-1, rec);
    427445    }
     446#ifdef HAS_ALMA
    428447    fillpm._update(n);
     448#endif
    429449  }
    430450  if (status > 0) {
    … …  
    432452    throw(AipsError("Reading error occured, data possibly corrupted."));
    433453  }
     454#ifdef HAS_ALMA
    434455  fillpm.done();
     456#endif
    435457  return status;
    436458}

  • branches/alma/src/STFiller.h

    r1495 r1603  
    111111  casa::Bool fileCheck() ;
    112112
     113  void setReferenceExpr(const std::string& rx) { refRx_ = rx; }
     114
    113115private:
    114116
    … …  
    120122  casa::uInt ifOffset_, beamOffset_;
    121123  casa::Vector<casa::Bool> haveXPol_;
     124  casa::String refRx_;
    122125  NROReader *nreader_ ;
    123126  casa::Bool isNRO_ ;

  • branches/alma/src/STFrequencies.cpp

    r1446 r1603  
    157157  // get first row - there should only be one matching id
    158158  const TableRecord& rec = row.get(0);
    159 
    160159  return SpectralCoordinate( getFrame(true), rec.asDouble("REFVAL"),
    161160                             rec.asDouble("INCREMENT"),
    … …  
    165164/**
    166165SpectralCoordinate
    167   asap::STFrequencies::getSpectralCoordinate( const MDirection& md,
    168                                               const MPosition& mp,
    169                                               const MEpoch& me,
    170                                               Double restfreq, uInt id ) const
     166  STFrequencies::getSpectralCoordinate( const MDirection& md,
     167                                        const MPosition& mp,
     168                                        const MEpoch& me,
     169                                        Double restfreq, uInt id ) const
    171170**/
    172171SpectralCoordinate

  • branches/alma/src/STHeader.cpp

    r901 r1603  
    5353  conforms = (this->antennaname == other.antennaname
    5454              && this->equinox == other.equinox
    55               && this->obstype == other.obstype
    5655              && this->fluxunit == other.fluxunit
    5756              );
    5857  return conforms;
     58}
     59
     60String STHeader::diff( const STHeader& other )
     61{
     62  ostringstream thediff;
     63  if ( this->equinox != other.equinox ) {
     64    thediff  << "Equinox: "  << this->equinox << " <-> "
     65             << other.equinox << endl;
     66  }
     67  if ( this->obstype != other.obstype ) {
     68    thediff << "Obs. Type: " << this->obstype << " <-> "
     69            << other.obstype << endl;
     70  }
     71  if ( this->fluxunit != other.fluxunit ) {
     72    thediff << "Flux unit: " << this->fluxunit << " <-> "
     73            << other.fluxunit << endl;
     74  }
     75  return String(thediff);
    5976}
    6077

  • branches/alma/src/STHeader.h

    r1386 r1603  
    5050
    5151  bool conformant(const STHeader& other);
     52  casa::String diff( const STHeader& other );
     53
    5254
    5355  casa::Int nchan;

  • branches/alma/src/STLineFinder.cpp

    r1315 r1603  
    870870       // iterator through lines
    871871       std::list<std::pair<int,int> >::const_iterator cli=lines.begin();
    872        for (int ch=0;ch<int(res_mask.size());++ch)
     872       for (int ch=0;ch<int(res_mask.size());++ch) {
    873873            if (ch<edge.first || ch>=edge.second) res_mask[ch]=false;
    874874            else if (!mask[ch]) res_mask[ch]=false;
    875875            else {
    876876                    res_mask[ch]=!invert; // no line by default
    877                     if (cli==lines.end()) continue;
    878                     if (ch>=cli->first && ch<cli->second)
    879                         res_mask[ch]=invert; // this is a line
    880                     if (ch>=cli->second)
    881                         ++cli; // next line in the list
    882                  }
    883 
     877                    if (cli!=lines.end())
     878                        if (ch>=cli->first && ch<cli->second)
     879                             res_mask[ch]=invert; // this is a line
     880            }
     881            if (cli!=lines.end())
     882                if (ch>=cli->second) {
     883                    ++cli; // next line in the list
     884                }
     885       }
    884886       return res_mask;
    885887  }

  • branches/alma/src/STMath.cpp

    r1516 r1603  
    297297                                             bool droprows)
    298298{
    299   if (insitu_) return in;
     299  if (insitu_) {
     300    return in;
     301  }
    300302  else {
    301303    // clone
    302     Scantable* tabp = new Scantable(*in, Bool(droprows));
    303     return CountedPtr<Scantable>(tabp);
     304    return CountedPtr<Scantable>(new Scantable(*in, Bool(droprows)));
    304305  }
    305306}
    … …  
    15131514  ArrayColumn<Float> specCol(tab, "SPECTRA");
    15141515  ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
     1516  ArrayColumn<Float> tsysCol(tab, "TSYS");
    15151517  for ( uInt i=0; i<tab.nrow(); ++i) {
    15161518    Float zdist = Float(C::pi_2) - elev(i);
    … …  
    15201522    specCol.put(i, ma.getArray());
    15211523    flagCol.put(i, flagsFromMA(ma));
     1524    Vector<Float> tsys;
     1525    tsysCol.get(i, tsys);
     1526    tsys *= factor;
     1527    tsysCol.put(i, tsys);
    15221528  }
    15231529  return out;
    … …  
    16141620  while ( it != in.end() ){
    16151621    if ( ! (*it)->conformant(*out) ) {
    1616       // log message: "ignoring scantable i, as it isn't
    1617       // conformant with the other(s)"
    1618       cerr << "oh oh" << endl;
    1619       ++it;
    1620       continue;
     1622      // non conformant.
     1623      pushLog(String("Warning: Can't merge scantables as header info differs."));
    16211624    }
    16221625    out->appendToHistoryTable((*it)->history());
    … …  
    16281631        Table thetab = freqit.table();
    16291632        uInt nrow = tout.nrow();
    1630         //tout.addRow(thetab.nrow());
     1633        tout.addRow(thetab.nrow());
    16311634        TableCopy::copyRows(tout, thetab, nrow, 0, thetab.nrow());
    16321635        ROTableRow row(thetab);
    … …  
    18261829
    18271830  InterpolateArray1D<Double,Float>::InterpolationMethod interp = stringToIMethod(method);
     1831  /*
     1832  // Comment from MV.
     1833  // the following code has been commented out because different FREQ_IDs have to be aligned together even
     1834  // if the frame doesn't change. So far, lack of this check didn't cause any problems.
    18281835  // test if user frame is different to base frame
    18291836  if ( in->frequencies().getFrameString(true)
    … …  
    18321839                    " (use set_freqframe) or it is aligned already."));
    18331840  }
     1841  */
    18341842  MFrequency::Types system = in->frequencies().getFrame();
    18351843  MVTime mvt(refEpoch.getValue());
    … …  
    18571865
    18581866    ROArrayColumn<Float> sCol(t, "SPECTRA");
    1859     MDirection direction = dirCol(0);
    1860     uInt nchan = sCol(0).nelements();
     1867    const MDirection direction = dirCol(0);
     1868    const uInt nchan = sCol(0).nelements();
     1869
     1870    // skip operations if there is nothing to align
     1871    if (fiter.pastEnd()) {
     1872        continue;
     1873    }
     1874
     1875    Table ftab = fiter.table();
     1876    // align all frequency ids with respect to the first encountered id
     1877    ScalarColumn<uInt> freqidCol(ftab, "FREQ_ID");
     1878    // get the SpectralCoordinate for the freqid, which we are iterating over
     1879    SpectralCoordinate sC = in->frequencies().getSpectralCoordinate(freqidCol(0));
     1880    FrequencyAligner<Float> fa( sC, nchan, refEpoch,
     1881                                direction, refPos, system );
     1882    // realign the SpectralCoordinate and put into the output Scantable
     1883    Vector<String> units(1);
     1884    units = String("Hz");
     1885    Bool linear=True;
     1886    SpectralCoordinate sc2 = fa.alignedSpectralCoordinate(linear);
     1887    sc2.setWorldAxisUnits(units);
     1888    const uInt id = out->frequencies().addEntry(sc2.referencePixel()[0],
     1889                                                sc2.referenceValue()[0],
     1890                                                sc2.increment()[0]);
    18611891    while ( !fiter.pastEnd() ) {
    1862       Table ftab = fiter.table();
    1863       ScalarColumn<uInt> freqidCol(ftab, "FREQ_ID");
    1864       // get the SpectralCoordinate for the freqid, which we are iterating over
    1865       SpectralCoordinate sC = in->frequencies().getSpectralCoordinate(freqidCol(0));
    1866       FrequencyAligner<Float> fa( sC, nchan, refEpoch,
    1867                                   direction, refPos, system );
    1868       // realign the SpectralCoordinate and put into the output Scantable
    1869       Vector<String> units(1);
    1870       units = String("Hz");
    1871       Bool linear=True;
    1872       SpectralCoordinate sc2 = fa.alignedSpectralCoordinate(linear);
    1873       sc2.setWorldAxisUnits(units);
    1874       uInt id = out->frequencies().addEntry(sc2.referencePixel()[0],
    1875                                             sc2.referenceValue()[0],
    1876                                             sc2.increment()[0]);
    1877       TableVector<uInt> tvec(ftab, "FREQ_ID");
    1878       tvec = id;
     1892      ftab = fiter.table();
     1893      // spectral coordinate for the current FREQ_ID
     1894      ScalarColumn<uInt> freqidCol2(ftab, "FREQ_ID");
     1895      sC = in->frequencies().getSpectralCoordinate(freqidCol2(0));
    18791896      // create the "global" abcissa for alignment with same FREQ_ID
    18801897      Vector<Double> abc(nchan);
    1881       Double w;
    18821898      for (uInt i=0; i<nchan; i++) {
    1883         sC.toWorld(w,Double(i));
    1884         abc[i] = w;
    1885       }
     1899           Double w;
     1900           sC.toWorld(w,Double(i));
     1901           abc[i] = w;
     1902      }
     1903      TableVector<uInt> tvec(ftab, "FREQ_ID");
     1904      // assign new frequency id to all rows
     1905      tvec = id;
    18861906      // cache abcissa for same time stamps, so iterate over those
    18871907      TableIterator timeiter(ftab, "TIME");
    … …  
    18921912        MEpoch::ROScalarColumn timeCol(tab, "TIME");
    18931913        // use align abcissa cache after the first row
     1914        // these rows should be just be POLNO
    18941915        bool first = true;
    1895         // these rows should be just be POLNO
    18961916        for (int i=0; i<int(tab.nrow()); ++i) {
    18971917          // input values

  • branches/alma/src/STMath.h

    r1516 r1603  
    6565    * average a vector of Scantables
    6666    * @param in the vector of Scantables to average
    67     * @param an optional mask to apply on specific weights
     67    * @param mask an optional mask to apply on specific weights
    6868    * @param weight weighting scheme
    6969    * @param avmode the mode ov averaging. Per "SCAN" or "ALL".

  • branches/alma/src/STSelector.cpp

    r1387 r1603  
    143143  if ( taql_.size() > 0 ) {
    144144    Table tmpt = tab;
     145    std::string pytaql = "USING STYLE PYTHON " + taql_;
    145146
    146147    if ( !query.isNull() ) { // taql and selection
    147       tmpt = tableCommand(taql_, tab(query));
     148      tmpt = tableCommand(pytaql, tab(query));
    148149    } else { // taql only
    149       tmpt = tableCommand(taql_, tab);
     150      tmpt = tableCommand(pytaql, tab);
    150151    }
    151152    return sort(tmpt);

  • branches/alma/src/STWeather.cpp

    r856 r1603  
    108108void STWeather::getEntry( Float& temperature, Float& pressure,
    109109                          Float& humidity, Float& windspeed, Float& windaz,
    110                           uInt id )
     110                          uInt id ) const
    111111{
    112112  Table t = table_(table_.col("ID") == Int(id) );

  • branches/alma/src/STWeather.h

    r1353 r1603  
    4444                       casa::Float& humidity,
    4545                       casa::Float& windspeed, casa::Float& windaz,
    46                        casa::uInt id);
     46                       casa::uInt id) const;
    4747
    4848  const casa::String& name() const { return name_; }

  • branches/alma/src/STWriter.cpp

    r1446 r1603  
    3939#include <casa/Utilities/Assert.h>
    4040
     41#include <atnf/PKSIO/PKSrecord.h>
    4142#include <atnf/PKSIO/PKSMS2writer.h>
    4243#include <atnf/PKSIO/PKSSDwriter.h>
    … …  
    4748#include <tables/Tables/ArrayColumn.h>
    4849
    49 //#include "SDFITSImageWriter.h"
     50#include "STFITSImageWriter.h"
    5051#include "STAsciiWriter.h"
    5152#include "STHeader.h"
    … …  
    5859STWriter::STWriter(const std::string &format)
    5960{
     61  format_ = format;
     62  String t(format_);
     63  t.upcase();
     64  if (t == "MS2") {
     65    writer_ = new PKSMS2writer();
     66  } else if (t == "SDFITS") {
     67    writer_ = new PKSSDwriter();
     68  } else if (t == "ASCII" || t == "FITS" || t == "CLASS") {
     69    writer_ = 0;
     70  } else {
     71    throw (AipsError("Unrecognized export format"));
     72  }
     73}
     74
     75STWriter::~STWriter()
     76{
     77   if (writer_) {
     78     delete writer_;
     79   }
     80}
     81
     82Int STWriter::setFormat(const std::string &format)
     83{
     84  if (format != format_) {
     85    if (writer_) delete writer_;
     86  }
     87
    6088  format_ = format;
    6189  String t(format_);
    … …  
    6593  } else if (t== "SDFITS") {
    6694    writer_ = new PKSSDwriter();
    67   } else if (t== "ASCII") {
    68     writer_ = 0;
    69   } else {
    70     throw (AipsError("Unrecognized export format"));
    71   }
    72 }
    73 
    74 STWriter::~STWriter()
    75 {
    76    if (writer_) {
    77      delete writer_;
    78    }
    79 }
    80 
    81 Int STWriter::setFormat(const std::string &format)
    82 {
    83   if (format != format_) {
    84     if (writer_) delete writer_;
    85   }
    86 
    87   format_ = format;
    88   String t(format_);
    89   t.upcase();
    90   if (t== "MS2") {
    91     writer_ = new PKSMS2writer();
    92   } else if (t== "SDFITS") {
    93     writer_ = new PKSSDwriter();
    94   } else if (t== "ASCII") {
     95  } else if (t == "ASCII" || t == "FITS" || t == "CLASS") {
    9596    writer_ = 0;
    9697  } else {
    … …  
    110111    } else {
    111112      return 1;
     113    }
     114  } else if ( format_ == "FITS" || format_ == "CLASS") {
     115    STFITSImageWriter iw;
     116    if (format_ == "CLASS") {
     117      iw.setClass(True);
     118    }
     119    if (iw.write(*in, filename)) {
     120      return 0;
    112121    }
    113122  }
    … …  
    139148  Int status;
    140149  status = writer_->create(String(filename), hdr.observer, hdr.project,
    141                                hdr.antennaname, hdr.antennaposition,
    142                                hdr.obstype, hdr.equinox, hdr.freqref,
    143                                nChan, nPol, havexpol, False, fluxUnit);
     150                           hdr.antennaname, hdr.antennaposition,
     151                           hdr.obstype, hdr.fluxunit,
     152                           hdr.equinox, hdr.freqref,
     153                           nChan, nPol, havexpol, False);
    144154  if ( status ) {
    145155    throw(AipsError("Failed to create output file"));
    146156  }
    147157
    148   Double          srcVel;
    149 
    150   String          fieldName, srcName, tcalTime;
    151   Vector<Float>   calFctr, sigma, tcal, tsys;
    152   Vector<Double>  direction(2), scanRate(2), srcDir(2), srcPM(2);
    153   Matrix<Float>   spectra;
    154   Matrix<uChar>   flagtra;
    155   Complex         xCalFctr;
     158
    156159  Int count = 0;
    157   Int scanno = 1;
     160  PKSrecord pksrec;
     161  pksrec.scanNo = 1;
    158162  // use spearate iterators to ensure renumbering of all numbers
    159163  TableIterator scanit(table, "SCANNO");
    … …  
    161165    Table stable = scanit.table();
    162166    TableIterator beamit(stable, "BEAMNO");
    163     Int beamno = 1;
     167    pksrec.beamNo = 1;
    164168    while (!beamit.pastEnd() ) {
    165169      Table btable = beamit.table();
    166       // position only varies by beam
    167       // No, we want to pointing data which varies by cycle!
    168170      MDirection::ScalarColumn dirCol(btable, "DIRECTION");
    169       Vector<Double> direction = dirCol(0).getAngle("rad").getValue();
     171      pksrec.direction = dirCol(0).getAngle("rad").getValue();
    170172      TableIterator cycit(btable, "CYCLENO");
    171173      ROArrayColumn<Double> srateCol(btable, "SCANRATE");
    172       srateCol.get(0, scanRate);
     174      Vector<Double> sratedbl;
     175      srateCol.get(0, sratedbl);
     176      Vector<Float> srateflt(sratedbl.nelements());
     177      convertArray(srateflt, sratedbl);
     178      //pksrec.scanRate = srateflt;
     179      pksrec.scanRate = sratedbl;
    173180      ROArrayColumn<Double> spmCol(btable, "SRCPROPERMOTION");
    174       spmCol.get(0, srcPM);
     181      spmCol.get(0, pksrec.srcPM);
    175182      ROArrayColumn <Double> sdirCol(btable, "SRCDIRECTION");
    176       sdirCol.get(0, srcDir);
     183      sdirCol.get(0, pksrec.srcDir);
    177184      ROScalarColumn<Double> svelCol(btable, "SRCVELOCITY");
    178       svelCol.get(0, srcVel);
     185      svelCol.get(0, pksrec.srcVel);
    179186      ROScalarColumn<uInt> bCol(btable, "BEAMNO");
    180       beamno = bCol(0)+1;
    181       Int cycno = 1;
     187      pksrec.beamNo = bCol(0)+1;
     188      pksrec.cycleNo = 1;
    182189      while (!cycit.pastEnd() ) {
    183190        Table ctable = cycit.table();
     191        TableIterator ifit(ctable, "IFNO");
    184192        //MDirection::ScalarColumn dirCol(ctable, "DIRECTION");
    185         //Vector<Double> direction = dirCol(0).getAngle("rad").getValue();
    186         TableIterator ifit(ctable, "IFNO");
    187         Int ifno = 1;
     193        //pksrec.direction = dirCol(0).getAngle("rad").getValue();
     194        pksrec.IFno = 1;
    188195        while (!ifit.pastEnd() ) {
    189196          Table itable = ifit.table();
    … …  
    192199          const TableRecord& rec = row.get(0);
    193200          ROArrayColumn<Float> specCol(itable, "SPECTRA");
    194           ifno = rec.asuInt("IFNO")+1;
     201          pksrec.IFno = rec.asuInt("IFNO")+1;
    195202          uInt nchan = specCol(0).nelements();
    196           //Double cdelt,crval,crpix, restfreq;
    197           Double cdelt,crval,crpix;
    198           Vector<Double> restfreq;
    199           Float focusAxi, focusTan, focusRot,
    200                 temperature, pressure, humidity, windSpeed, windAz;
     203          Double crval,crpix;
     204          //Vector<Double> restfreq;
    201205          Float tmp0,tmp1,tmp2,tmp3,tmp4;
    202           Vector<Float> tcalval;
    203           //String stmp0,stmp1, tcalt;
    204206          String tcalt;
    205207          Vector<String> stmp0, stmp1;
    206           in->frequencies().getEntry(crpix,crval,cdelt, rec.asuInt("FREQ_ID"));
    207           in->focus().getEntry(focusAxi, focusTan, focusRot,
    208                                tmp0,tmp1,tmp2,tmp3,tmp4,
     208          in->frequencies().getEntry(crpix,crval, pksrec.freqInc,
     209                                     rec.asuInt("FREQ_ID"));
     210          in->focus().getEntry(pksrec.focusAxi, pksrec.focusTan,
     211                               pksrec.focusRot, tmp0,tmp1,tmp2,tmp3,tmp4,
    209212                               rec.asuInt("FOCUS_ID"));
    210           in->molecules().getEntry(restfreq,stmp0,stmp1,rec.asuInt("MOLECULE_ID"));
    211           in->tcal().getEntry(tcalt,tcalval,rec.asuInt("TCAL_ID"));
    212           in->weather().getEntry(temperature, pressure, humidity,
    213                                  windSpeed, windAz,
    214                                  rec.asuInt("WEATHER_ID"));
     213          in->molecules().getEntry(pksrec.restFreq,stmp0,stmp1,
     214                                   rec.asuInt("MOLECULE_ID"));
     215          in->tcal().getEntry(pksrec.tcalTime, pksrec.tcal,
     216                              rec.asuInt("TCAL_ID"));
     217          in->weather().getEntry(pksrec.temperature, pksrec.pressure,
     218                                 pksrec.humidity, pksrec.windSpeed,
     219                                 pksrec.windAz, rec.asuInt("WEATHER_ID"));
    215220          Double pixel = Double(nchan/2);
    216           Double refFreqNew = (pixel-crpix)*cdelt + crval;
     221          pksrec.refFreq = (pixel-crpix)*pksrec.freqInc + crval;
    217222          // ok, now we have nrows for the n polarizations in this table
    218           Matrix<Float> specs;
    219           Matrix<uChar> flags;
    220           Vector<Complex> xpol;
    221           polConversion(specs, flags, xpol, itable);
    222           Vector<Float> tsys = tsysFromTable(itable);
     223          polConversion(pksrec.spectra, pksrec.flagged, pksrec.xPol, itable);
     224          pksrec.tsys = tsysFromTable(itable);
    223225          // dummy data
    224           //uInt npol;
    225           //if ( hdr.antennaname == "GBT" ) {
    226           //  npol = nPolUsed;
    227           //}
    228           //else {   
    229           //  npol = specs.ncolumn();
    230           //}
    231           uInt npol = specs.ncolumn();
    232 
    233           Matrix<Float>   baseLin(npol,2, 0.0f);
    234           Matrix<Float>   baseSub(npol,9, 0.0f);
    235           Complex         xCalFctr;
    236           Vector<Double>  scanRate(2, 0.0);
    237           Vector<Float>   sigma(npol, 0.0f);
    238           Vector<Float>   calFctr(npol, 0.0f);
    239           status = writer_->write(scanno, cycno, rec.asDouble("TIME"),
    240                                       rec.asDouble("INTERVAL"),
    241                                       rec.asString("FIELDNAME"),
    242                                       rec.asString("SRCNAME"),
    243                                       srcDir, srcPM, srcVel,hdr.obstype,
    244                                       ifno,
    245                                       refFreqNew, nchan*abs(cdelt), cdelt,
    246                                       restfreq,
    247                                       tcalval,
    248                                       tcalt,
    249                                       rec.asFloat("AZIMUTH"),
    250                                       rec.asFloat("ELEVATION"),
    251                                       rec.asFloat("PARANGLE"),
    252                                       focusAxi, focusTan, focusRot,
    253                                       temperature,
    254                                       pressure, humidity, windSpeed, windAz,
    255                                       rec.asInt("REFBEAMNO")+1, beamno,
    256                                       direction,
    257                                       scanRate,
    258                                       tsys,
    259                                       sigma, calFctr,// not in scantable
    260                                       baseLin, baseSub,// not in scantable
    261                                       specs, flags,
    262                                       xCalFctr,//
    263                                       xpol);
     226          uInt npol = pksrec.spectra.ncolumn();
     227
     228          pksrec.mjd       = rec.asDouble("TIME");
     229          pksrec.interval  = rec.asDouble("INTERVAL");
     230          pksrec.fieldName = rec.asString("FIELDNAME");
     231          pksrec.srcName   = rec.asString("SRCNAME");
     232          pksrec.obsType   = hdr.obstype;
     233          pksrec.bandwidth = nchan * abs(pksrec.freqInc);
     234          pksrec.azimuth   = rec.asFloat("AZIMUTH");
     235          pksrec.elevation = rec.asFloat("ELEVATION");
     236          pksrec.parAngle  = rec.asFloat("PARANGLE");
     237          pksrec.refBeam   = rec.asInt("REFBEAMNO") + 1;
     238          pksrec.sigma.resize(npol);
     239          pksrec.sigma     = 0.0f;
     240          pksrec.calFctr.resize(npol);
     241          pksrec.calFctr   = 0.0f;
     242          pksrec.baseLin.resize(npol,2);
     243          pksrec.baseLin   = 0.0f;
     244          pksrec.baseSub.resize(npol,9);
     245          pksrec.baseSub   = 0.0f;
     246          pksrec.xCalFctr  = 0.0;
     247
     248          status = writer_->write(pksrec);
    264249          if ( status ) {
    265250            writer_->close();
    … …  
    267252          }
    268253          ++count;
    269           ++ifno;
     254          //++pksrec.IFno;
    270255          ++ifit;
    271256        }
    272         ++cycno;
     257        ++pksrec.cycleNo;
    273258        ++cycit;
    274259      }
    275       ++beamno;
     260      //++pksrec.beamNo;
    276261      ++beamit;
    277262    }
    278     ++scanno;
     263    ++pksrec.scanNo;
    279264    ++scanit;
    280265  }
    … …  
    286271  if ( format_ == "MS2" ) {
    287272    replacePtTab(table, filename);
    288   } 
     273  }
    289274  return 0;
    290275}

  • branches/alma/src/Scantable.cpp

    r1451 r1603  
    213213  td.addColumn(ScalarColumnDesc<String>("SRCNAME"));
    214214  // Type of source (on=0, off=1, other=-1)
    215   td.addColumn(ScalarColumnDesc<Int>("SRCTYPE", Int(-1)));
     215  ScalarColumnDesc<Int> stypeColumn("SRCTYPE");
     216  stypeColumn.setDefault(Int(-1));
     217  td.addColumn(stypeColumn);
    216218  td.addColumn(ScalarColumnDesc<String>("FIELDNAME"));
    217219
    … …  
    535537    return int(n);
    536538  } else {
    537     // take the first SCANNO,POLNO,BEAMNO,CYCLENO as nbeam shouldn't vary with these
     539    // take the first SCANNO,POLNO,BEAMNO,CYCLENO as nbeam shouldn't
     540    // vary with these
    538541    Table t = table_(table_.col("IFNO") == ifno);
    539542    if ( t.nrow() == 0 ) return 0;
    … …  
    885888}
    886889
     890MEpoch Scantable::getEpoch(int whichrow) const
     891{
     892  if (whichrow > -1) {
     893    return timeCol_(uInt(whichrow));
     894  } else {
     895    Double tm;
     896    table_.keywordSet().get("UTC",tm);
     897    return MEpoch(MVEpoch(tm)); 
     898  }
     899}
     900
    887901std::string Scantable::getDirectionString(int whichrow) const
    888902{
    … …  
    892906std::vector< double > Scantable::getAbcissa( int whichrow ) const
    893907{
    894   if ( whichrow > int(table_.nrow()) ) throw(AipsError("Illegal ro number"));
     908  if ( whichrow > int(table_.nrow()) ) throw(AipsError("Illegal row number"));
    895909  std::vector<double> stlout;
    896910  int nchan = specCol_(whichrow).nelements();

  • branches/alma/src/Scantable.h

    r1446 r1603  
    158158  casa::MEpoch::Types getTimeReference() const;
    159159
     160
     161  casa::MEpoch getEpoch(int whichrow) const;
     162
    160163  /**
    161164   * Get global antenna position
    … …  
    164167  casa::MPosition getAntennaPosition() const;
    165168
    166         /**
    167         * the @ref casa::MDirection for a specific row
    168         * @param[in] whichrow the row number
    169         * return casa::MDirection
    170         */
     169  /**
     170  * the @ref casa::MDirection for a specific row
     171  * @param[in] whichrow the row number
     172  * return casa::MDirection
     173  */
    171174  casa::MDirection getDirection( int whichrow ) const;
    172 
    173         /**
    174          * get the direction as a string
    175         * @param[in] whichrow the row number
    176         * return the direction string
    177         */
     175 
     176  /**
     177   * get the direction type as a string, e.g. "J2000"
     178  * @param[in] whichrow the row number
     179  * return the direction string
     180  */
    178181  std::string getDirectionString( int whichrow ) const;
    179182
    180         /**
    181         * set the direction type as a string, e.g. "J2000"
    182         * @param[in] refstr the direction type
    183         */
     183  /**
     184  * set the direction type as a string, e.g. "J2000"
     185  * @param[in] refstr the direction type
     186  */
    184187  void setDirectionRefString(const std::string& refstr="");
     188
    185189  /**
    186190   * get the direction reference string
    187191   * @return a string describing the direction reference
    188192   */
    189   std::string getDirectionRefString() const;    /**
    190          * get the direction type as a string, e.g. "J2000"
    191          * param[in] whichrow the row number
    192          * return the direction string
    193          */
    194 
     193  std::string getDirectionRefString() const;   
    195194
    196195  /**

  • branches/alma/src/python_STFiller.cpp

    r1386 r1603  
    4747        .def("_close", &STFillerWrapper::close)
    4848        .def("_getdata", &STFillerWrapper::getScantable)
     49        .def("_setreferenceexpr", &STFillerWrapper::setReferenceExpr)
    4950      ;
    5051    };

  • branches/alma/src/python_asap.cpp

    r1126 r1603  
    3939#include "ScantableWrapper.h"
    4040
     41#ifndef HAVE_LIBPYRAP
     42  #include "pyconversions.h"
     43#else
     44  #include <pyrap/Converters/PycExcp.h>
     45  #include <pyrap/Converters/PycBasicData.h>
     46#endif
    4147
    42 #include "pyconversions.h"
    4348#include "python_asap.h"
    4449
     50#ifndef HAVE_LIBPYRAP
    4551namespace asap {
    4652  namespace python {
     53
    4754void translate_ex(const casa::AipsError& e)
    4855{
    … …  
    5360  }
    5461}
     62#endif
     63
    5564using namespace boost::python;
    5665
    … …  
    6574  asap::python::python_STWriter();
    6675  asap::python::python_LineCatalog();
     76  asap::python::python_Logger();
    6777
    68   asap::python::python_Logger();
     78#ifndef HAVE_LIBPYRAP
     79  // Use built-in pyconversions.h
    6980  register_exception_translator<casa::AipsError>(&asap::python::translate_ex);
    70 
    71   //std_vector_to_tuple <  > ();
    7281  from_python_sequence < std::vector< asap::ScantableWrapper >,
    7382    variable_capacity_policy > ();
    74 
    7583  std_vector_to_tuple < int > ();
    7684  from_python_sequence < std::vector < int >,
    … …  
    9199  from_python_sequence < std::vector < bool >,
    92100    variable_capacity_policy > ();
     101#else
     102  casa::pyrap::register_convert_excp();
     103  casa::pyrap::register_convert_basicdata();
     104  casa::pyrap::register_convert_std_vector<asap::ScantableWrapper>();
     105#endif
    93106}
Note: See TracChangeset for help on using the changeset viewer.