Changeset 1446


Ignore:
Timestamp:
11/12/08 17:04:01 (16 years ago)
Author:
TakTsutsumi
Message:

Merged recent updates (since 2007) from nrao-asap

Location:
branches/alma
Files:
2 added
25 edited

Legend:

Unmodified
Added
Removed
  • branches/alma/python/asapfitter.py

    r1389 r1446  
    586586            scan = self.data
    587587        rows = xrange(scan.nrow())
     588        # Save parameters of baseline fits as a class attribute.
     589        # NOTICE: This does not reflect changes in scantable!
     590        if len(rows) > 0: self.blpars=[]
    588591        from asap import asaplog
    589592        asaplog.push("Fitting:")
     
    595598            self.data = None
    596599            self.fit()
    597             x = self.get_parameters()
     600            fpar = self.get_parameters()
    598601            if plot:
    599602                self.plot(residual=True)
    600603                x = raw_input("Accept fit ([y]/n): ")
    601604                if x.upper() == 'N':
     605                    self.blpars.append(None)
    602606                    continue
    603607            scan._setspectrum(self.fitter.getresidual(), r)
     608            self.blpars.append(fpar)
    604609        if plot:
    605610            self._p.unmap()
  • branches/alma/python/asaplotbase.py

    r1259 r1446  
    161161        y2 = range(l2)
    162162        m2 = range(l2)
    163         ymsk = y.raw_mask()
    164         ydat = y.raw_data()
     163        #ymsk = y.raw_mask()
     164        #ydat = y.raw_data()
     165        ymsk = y.mask
     166        ydat = y.data
    165167        for i in range(l2):
    166168            x2[i] = x[i/2]
     
    707709            for sp in self.subplots:
    708710                ax = sp['axes']
    709                 s = ax.title.get_size()
    710                 tsize = s-(self.cols+self.rows)
    711                 ax.title.set_size(tsize)
     711                s = rcParams['axes.titlesize']
     712                tsize = s-(self.cols+self.rows-1)
     713                ax.title.set_size(max(tsize,9))
    712714                fp = FP(size=rcParams['axes.labelsize'])
    713715                setp(ax.get_xticklabels(), fontsize=xts)
  • branches/alma/python/asapmath.py

    r1389 r1446  
    4545    if kwargs.has_key('align'):
    4646        align = kwargs.get('align')
     47    compel = False
     48    if kwargs.has_key('compel'):
     49        compel = kwargs.get('compel')
    4750    varlist = vars()
    4851    if isinstance(args[0],list):
     
    8689        del merged
    8790    else:
    88         s = scantable(stm._average(alignedlst, mask, weight.upper(), scanav))
     91        #s = scantable(stm._average(alignedlst, mask, weight.upper(), scanav))
     92        s = scantable(stm._new_average(alignedlst, compel, mask, weight.upper(), scanav))
    8993    s._add_history("average_time",varlist)
    9094    print_log()
  • branches/alma/python/asapplotter.py

    r1389 r1446  
    605605        if n > 1:
    606606            ganged = rcParams['plotter.ganged']
     607            ###Start Mod: 2008.09.22 kana ###
     608            if self._panelling == 'i':
     609                ganged = False
     610            ###End Mod#######################
    607611            if self._rows and self._cols:
    608612                n = min(n,self._rows*self._cols)
     
    688692                if stackcount == nstack:
    689693                    allxlim.sort()
    690                     self._plotter.axes.set_xlim([allxlim[0],allxlim[-1]])
     694                    ###Start Mod: 2008.09.22 kana ###
     695                    #self._plotter.axes.set_xlim([allxlim[0],allxlim[-1]])
     696                    self._plotter.subplots[panelcount-1]['axes'].set_xlim([allxlim[0],allxlim[-1]])
     697                    ###End Mod#######################
    691698                    # clear
    692699                    allxlim =[]
     
    741748        return userlabel or d[mode]
    742749
    743     def plotazel(self, scan=None):
     750    def plotazel(self, scan=None, outfile=None):
    744751        """
    745752        plot azimuth and elevation  versus time of a scantable
     
    750757        from matplotlib.numerix import array, pi
    751758        self._data = scan
     759        self._outfile = outfile
    752760        dates = self._data.get_time()
    753761        t = PL.date2num(dates)
    754762        tz = timezone('UTC')
    755763        PL.cla()
    756         PL.ioff()
     764        #PL.ioff()
    757765        PL.clf()
    758766        tdel = max(t) - min(t)
     
    772780            minloc = MinuteLocator(20)
    773781        PL.title(dstr)
    774         PL.plot_date(t,el,'b,', tz=tz)
    775         #ax.grid(True)
     782
     783        if tdel == 0.0:
     784            th = (t - PL.floor(t))*24.0
     785            PL.plot(th,el,'o',markersize=2, markerfacecolor='b', markeredgecolor='b')
     786        else:
     787            PL.plot_date(t,el,'o', markersize=2, markerfacecolor='b', markeredgecolor='b',tz=tz)
     788            #ax.grid(True)
     789            ax.xaxis.set_major_formatter(timefmt)
     790            ax.xaxis.set_major_locator(majloc)
     791            ax.xaxis.set_minor_locator(minloc)
    776792        ax.yaxis.grid(True)
    777793        yloc = MultipleLocator(30)
    778794        ax.set_ylim(0,90)
    779         ax.xaxis.set_major_formatter(timefmt)
    780         ax.xaxis.set_major_locator(majloc)
    781         ax.xaxis.set_minor_locator(minloc)
    782795        ax.yaxis.set_major_locator(yloc)
    783796        if tdel > 1.0:
     
    785798        #    PL.setp(labels, fontsize=10, rotation=45)
    786799            PL.setp(labels, fontsize=10)
     800
    787801        # Az plot
    788802        az = array(self._data.get_azimuth())*180./pi
     
    792806
    793807        ax = PL.subplot(2,1,2)
    794         PL.xlabel('Time (UT)')
     808        #PL.xlabel('Time (UT [hour])')
    795809        PL.ylabel('Az [deg.]')
    796         PL.plot_date(t,az,'b,', tz=tz)
     810        if tdel == 0.0:
     811            PL.plot(th,az,'o',markersize=2, markeredgecolor='b',markerfacecolor='b')
     812        else:
     813            PL.plot_date(t,az,'o', markersize=2,markeredgecolor='b',markerfacecolor='b',tz=tz)
     814            ax.xaxis.set_major_formatter(timefmt)
     815            ax.xaxis.set_major_locator(majloc)
     816            ax.xaxis.set_minor_locator(minloc)
     817        #ax.grid(True)
    797818        ax.set_ylim(0,360)
    798         #ax.grid(True)
    799819        ax.yaxis.grid(True)
    800820        #hfmt = DateFormatter('%H')
    801821        #hloc = HourLocator()
    802822        yloc = MultipleLocator(60)
    803         ax.xaxis.set_major_formatter(timefmt)
    804         ax.xaxis.set_major_locator(majloc)
    805         ax.xaxis.set_minor_locator(minloc)
    806823        ax.yaxis.set_major_locator(yloc)
    807824        if tdel > 1.0:
    808825            labels = ax.get_xticklabels()
    809826            PL.setp(labels, fontsize=10)
    810         PL.ion()
     827            PL.xlabel('Time (UT [day])')
     828        else:
     829            PL.xlabel('Time (UT [hour])')
     830
     831        #PL.ion()
    811832        PL.draw()
    812 
    813     def plotpointing(self, scan=None):
     833        if (self._outfile is not None):
     834           PL.savefig(self._outfile)
     835
     836    def plotpointing(self, scan=None, outfile=None):
    814837        """
    815838        plot telescope pointings
     
    820843        from matplotlib.numerix import array, pi, zeros
    821844        self._data = scan
     845        self._outfile = outfile
    822846        dir = array(self._data.get_directionval()).transpose()
    823847        ra = dir[0]*180./pi
    824848        dec = dir[1]*180./pi
    825849        PL.cla()
    826         PL.ioff()
     850        #PL.ioff()
    827851        PL.clf()
    828852        ax = PL.axes([0.1,0.1,0.8,0.8])
     
    835859        [xmin,xmax,ymin,ymax] = PL.axis()
    836860        PL.axis([xmax,xmin,ymin,ymax])
    837         PL.ion()
     861        #PL.ion()
    838862        PL.draw()
    839 
     863        if (self._outfile is not None):
     864           PL.savefig(self._outfile)
     865
     866    # plot total power data
     867    # plotting in time is not yet implemented..
     868    def plottp(self, scan=None, outfile=None):
     869        if self._plotter.is_dead:
     870            self._plotter = self._newplotter()
     871        self._plotter.hold()
     872        self._plotter.clear()
     873        from asap import scantable
     874        if not self._data and not scan:
     875            msg = "Input is not a scantable"
     876            if rcParams['verbose']:
     877                print msg
     878                return
     879            raise TypeError(msg)
     880        if isinstance(scan, scantable):
     881            if self._data is not None:
     882                if scan != self._data:
     883                    self._data = scan
     884                    # reset
     885                    self._reset()
     886            else:
     887                self._data = scan
     888                self._reset()
     889        # ranges become invalid when abcissa changes?
     890        #if self._abcunit and self._abcunit != self._data.get_unit():
     891        #    self._minmaxx = None
     892        #    self._minmaxy = None
     893        #    self._abcunit = self._data.get_unit()
     894        #    self._datamask = None
     895        self._plottp(self._data)
     896        if self._minmaxy is not None:
     897            self._plotter.set_limits(ylim=self._minmaxy)
     898        self._plotter.release()
     899        self._plotter.tidy()
     900        self._plotter.show(hardrefresh=False)
     901        print_log()
     902        return
     903
     904    def _plottp(self,scan):
     905        """
     906        private method for plotting total power data
     907        """
     908        from matplotlib.numerix import ma, array, arange, logical_not
     909        r=0
     910        nr = scan.nrow()
     911        a0,b0 = -1,-1
     912        allxlim = []
     913        allylim = []
     914        y=[]
     915        self._plotter.set_panels()
     916        self._plotter.palette(0)
     917        #title
     918        #xlab = self._abcissa and self._abcissa[panelcount] \
     919        #       or scan._getabcissalabel()
     920        #ylab = self._ordinate and self._ordinate[panelcount] \
     921        #       or scan._get_ordinate_label()
     922        xlab = self._abcissa or 'row number' #or Time
     923        ylab = self._ordinate or scan._get_ordinate_label()
     924        self._plotter.set_axes('xlabel',xlab)
     925        self._plotter.set_axes('ylabel',ylab)
     926        lbl = self._get_label(scan, r, 's', self._title)
     927        if isinstance(lbl, list) or isinstance(lbl, tuple):
     928        #    if 0 <= panelcount < len(lbl):
     929        #        lbl = lbl[panelcount]
     930        #    else:
     931                # get default label
     932             lbl = self._get_label(scan, r, self._panelling, None)
     933        self._plotter.set_axes('title',lbl)
     934        y=array(scan._get_column(scan._getspectrum,-1))
     935        m = array(scan._get_column(scan._getmask,-1))
     936        y = ma.masked_array(y,mask=logical_not(array(m,copy=False)))
     937        x = arange(len(y))
     938        # try to handle spectral data somewhat...
     939        l,m = y.shape
     940        if m > 1:
     941            y=y.mean(axis=1)
     942        plotit = self._plotter.plot
     943        llbl = self._get_label(scan, r, self._stacking, None)
     944        self._plotter.set_line(label=llbl)
     945        if len(x) > 0:
     946            plotit(x,y)
     947
     948
     949    # forwards to matplotlib.Figure.text
     950    def figtext(self, *args, **kwargs):
     951        """
     952        Add text to figure at location x,y (relative 0-1 coords).
     953        This method forwards *args and **kwargs to a Matplotlib method,
     954        matplotlib.Figure.text.
     955        See the method help for detailed information.
     956        """
     957        self._plotter.text(*args, **kwargs)
     958    # end matplotlib.Figure.text forwarding function
     959
  • branches/alma/python/scantable.py

    r1401 r1446  
    300300        self._setselection(selection)
    301301
     302    def get_row(self, row=0, insitu=None):
     303        """
     304        Select a row in the scantable.
     305        Return a scantable with single row.
     306        Parameters:
     307            row: row no of integration, default is 0.
     308            insitu: if False a new scantable is returned.
     309                    Otherwise, the scaling is done in-situ
     310                    The default is taken from .asaprc (False)
     311        """
     312        if insitu is None: insitu = rcParams['insitu']
     313        if not insitu:
     314            workscan = self.copy()
     315        else:
     316            workscan = self
     317        # Select a row
     318        sel=selector()
     319        sel.set_scans([workscan.getscan(row)])
     320        sel.set_cycles([workscan.getcycle(row)])
     321        sel.set_beams([workscan.getbeam(row)])
     322        sel.set_ifs([workscan.getif(row)])
     323        sel.set_polarisations([workscan.getpol(row)])
     324        sel.set_name(workscan._getsourcename(row))
     325        workscan.set_selection(sel)
     326        if not workscan.nrow() == 1:
     327            msg = "Cloud not identify single row. %d rows selected."%(workscan.nrow())
     328            raise RuntimeError(msg)
     329        del sel
     330        if insitu:
     331            self._assign(workscan)
     332        else:
     333            return workscan
     334
    302335    def stats(self, stat='stddev', mask=None):
    303336        """
     
    528561        return self._get_column(self._getdirectionvec, row)
    529562
    530 
    531563    def set_unit(self, unit='channel'):
    532564        """
     
    768800        return msk
    769801
    770     def get_restfreqs(self):
     802    def get_masklist(self, mask=None, row=0):
     803        """
     804        Compute and return a list of mask windows, [min, max].
     805        Parameters:
     806            mask:       channel mask, created with create_mask.
     807            row:        calcutate the masklist using the specified row
     808                        for unit conversions, default is row=0
     809                        only necessary if frequency varies over rows.
     810        Returns:
     811            [min, max], [min2, max2], ...
     812                Pairs of start/end points (inclusive)specifying
     813                the masked regions
     814        """
     815        if not (isinstance(mask,list) or isinstance(mask, tuple)):
     816            raise TypeError("The mask should be list or tuple.")
     817        if len(mask) < 2:
     818            raise TypeError("The mask elements should be > 1")
     819        if self.nchan() != len(mask):
     820            msg = "Number of channels in scantable != number of mask elements"
     821            raise TypeError(msg)
     822        data = self._getabcissa(row)
     823        u = self._getcoordinfo()[0]
     824        if rcParams['verbose']:
     825            if u == "": u = "channel"
     826            msg = "The current mask window unit is %s" % u
     827            i = self._check_ifs()
     828            if not i:
     829                msg += "\nThis mask is only valid for IF=%d" % (self.getif(i))
     830            asaplog.push(msg)
     831        masklist=[]
     832        ist, ien = None, None
     833        ist, ien=self.get_mask_indices(mask)
     834        if ist is not None and ien is not None:
     835            for i in xrange(len(ist)):
     836                range=[data[ist[i]],data[ien[i]]]
     837                range.sort()
     838                masklist.append([range[0],range[1]])
     839        return masklist
     840
     841    def get_mask_indices(self, mask=None):
     842        """
     843        Compute and Return lists of mask start indices and mask end indices.
     844         Parameters:
     845            mask:       channel mask, created with create_mask.
     846        Returns:
     847            List of mask start indices and that of mask end indices,
     848            i.e., [istart1,istart2,....], [iend1,iend2,....].
     849        """
     850        if not (isinstance(mask,list) or isinstance(mask, tuple)):
     851            raise TypeError("The mask should be list or tuple.")
     852        if len(mask) < 2:
     853            raise TypeError("The mask elements should be > 1")
     854        istart=[]
     855        iend=[]
     856        if mask[0]: istart.append(0)
     857        for i in range(len(mask)-1):
     858            if not mask[i] and mask[i+1]:
     859                istart.append(i+1)
     860            elif mask[i] and not mask[i+1]:
     861                iend.append(i)
     862        if mask[len(mask)-1]: iend.append(len(mask)-1)
     863        if len(istart) != len(iend):
     864            raise RuntimeError("Numbers of mask start != mask end.")
     865        for i in range(len(istart)):
     866            if istart[i] > iend[i]:
     867                raise RuntimeError("Mask start index > mask end index")
     868                break
     869        return istart,iend
     870
     871#    def get_restfreqs(self):
     872#        """
     873#        Get the restfrequency(s) stored in this scantable.
     874#        The return value(s) are always of unit 'Hz'
     875#        Parameters:
     876#            none
     877#        Returns:
     878#            a list of doubles
     879#        """
     880#        return list(self._getrestfreqs())
     881
     882    def get_restfreqs(self, ids=None):
    771883        """
    772884        Get the restfrequency(s) stored in this scantable.
    773885        The return value(s) are always of unit 'Hz'
    774886        Parameters:
    775             none
     887            ids: (optional) a list of MOLECULE_ID for that restfrequency(s) to
     888                 be retrieved
    776889        Returns:
    777             a list of doubles
    778         """
    779         return list(self._getrestfreqs())
    780 
     890            dictionary containing ids and a list of doubles for each id
     891        """
     892        if ids is None:
     893            rfreqs={}
     894            idlist = self.getmolnos()
     895            for i in idlist:
     896                rfreqs[i]=list(self._getrestfreqs(i))
     897            return rfreqs
     898        else:
     899            if type(ids)==list or type(ids)==tuple:
     900                rfreqs={}
     901                for i in ids:
     902                    rfreqs[i]=list(self._getrestfreqs(i))
     903                return rfreqs
     904            else:
     905                return list(self._getrestfreqs(ids))
     906            #return list(self._getrestfreqs(ids))
    781907
    782908    def set_restfreqs(self, freqs=None, unit='Hz'):
    783909        """
     910        ********NEED TO BE UPDATED begin************
    784911        Set or replace the restfrequency specified and
    785912        If the 'freqs' argument holds a scalar,
     
    793920        E.g. 'freqs=[1e9, 2e9]'  would mean IF 0 gets restfreq 1e9 and
    794921        IF 1 gets restfreq 2e9.
     922        ********NEED TO BE UPDATED end************
    795923        You can also specify the frequencies via a linecatalog/
    796924
     
    800928
    801929        Example:
    802             # set the given restfrequency for the whole table
     930            # set the given restfrequency for the all currently selected IFs
    803931            scan.set_restfreqs(freqs=1.4e9)
    804             # If thee number of IFs in the data is >= 2 IF0 gets the first
    805             # value IF1 the second...
    806             scan.set_restfreqs(freqs=[1.4e9, 1.67e9])
     932            # set multiple restfrequencies to all the selected data
     933            scan.set_restfreqs(freqs=[1.4e9, 1.41e9, 1.42e9])
     934            # If the number of IFs in the data is >= 2 the IF0 gets the first
     935            # value IF1 the second... NOTE that freqs needs to be
     936            # specified in list of list (e.g. [[],[],...] ).
     937            scan.set_restfreqs(freqs=[[1.4e9],[1.67e9]])
    807938            #set the given restfrequency for the whole table (by name)
    808939            scan.set_restfreqs(freqs="OH1667")
     
    824955        # simple  value
    825956        if isinstance(freqs, int) or isinstance(freqs, float):
    826             self._setrestfreqs(freqs, "",unit)
     957            # TT mod
     958            #self._setrestfreqs(freqs, "",unit)
     959            self._setrestfreqs([freqs], [""],unit)
    827960        # list of values
    828961        elif isinstance(freqs, list) or isinstance(freqs, tuple):
    829962            # list values are scalars
    830963            if isinstance(freqs[-1], int) or isinstance(freqs[-1], float):
     964                self._setrestfreqs(freqs, [""],unit)
     965            # list values are tuples, (value, name)
     966            elif isinstance(freqs[-1], dict):
     967                #sel = selector()
     968                #savesel = self._getselection()
     969                #iflist = self.getifnos()
     970                #for i in xrange(len(freqs)):
     971                #    sel.set_ifs(iflist[i])
     972                #    self._setselection(sel)
     973                #    self._setrestfreqs(freqs[i], "",unit)
     974                #self._setselection(savesel)
     975                self._setrestfreqs(freqs["value"],
     976                                   freqs["name"], "MHz")
     977            elif isinstance(freqs[-1], list) or isinstance(freqs[-1], tuple):
    831978                sel = selector()
    832979                savesel = self._getselection()
    833980                iflist = self.getifnos()
    834                 for i in xrange(len(freqs)):
    835                     sel.set_ifs(iflist[i])
    836                     self._setselection(sel)
    837                     self._setrestfreqs(freqs[i], "",unit)
    838                 self._setselection(savesel)
    839             # list values are tuples, (value, name)
    840             elif isinstance(freqs[-1], dict):
    841                 sel = selector()
    842                 savesel = self._getselection()
    843                 iflist = self.getifnos()
     981                if len(freqs)>len(iflist):
     982                    raise ValueError("number of elements in list of list exeeds the current IF selections")
    844983                for i in xrange(len(freqs)):
    845984                    sel.set_ifs(iflist[i])
     
    852991            sel = selector()
    853992            savesel = self._getselection()
    854             iflist = self.getifnos()
    855993            for i in xrange(freqs.nrow()):
    856994                sel.set_ifs(iflist[i])
     
    12541392            f = fitter()
    12551393            f.set_scan(self, mask)
    1256             #f.set_function(poly=order)
    12571394            if uselin:
    12581395                f.set_function(lpoly=order)
     
    12601397                f.set_function(poly=order)
    12611398            s = f.auto_fit(insitu, plot=plot)
     1399            # Save parameters of baseline fits as a class attribute.
     1400            # NOTICE: It does not reflect changes in scantable!
     1401            self.blpars = f.blpars
    12621402            s._add_history("poly_baseline", varlist)
    12631403            print_log()
     
    13551495
    13561496        rows = range(workscan.nrow())
     1497        # Save parameters of baseline fits & masklists as a class attribute.
     1498        # NOTICE: It does not reflect changes in scantable!
     1499        if len(rows) > 0:
     1500            self.blpars=[]
     1501            self.masklists=[]
    13571502        asaplog.push("Processing:")
    13581503        for r in rows:
     
    13711516            # setup line finder
    13721517            fl.find_lines(r, mask, curedge)
     1518            outmask=fl.get_mask()
    13731519            f.set_scan(workscan, fl.get_mask())
    13741520            f.x = workscan._getabcissa(r)
     
    13761522            f.data = None
    13771523            f.fit()
    1378             x = f.get_parameters()
     1524           
     1525            # Show mask list
     1526            masklist=workscan.get_masklist(fl.get_mask(),row=r)
     1527            msg = "mask range: "+str(masklist)
     1528            asaplog.push(msg, False)
     1529
     1530            fpar = f.get_parameters()
    13791531            if plot:
    13801532                f.plot(residual=True)
    13811533                x = raw_input("Accept fit ( [y]/n ): ")
    13821534                if x.upper() == 'N':
     1535                    self.blpars.append(None)
     1536                    self.masklists.append(None)
    13831537                    continue
    13841538            workscan._setspectrum(f.fitter.getresidual(), r)
     1539            self.blpars.append(fpar)
     1540            self.masklists.append(masklist)
    13851541        if plot:
    13861542            f._p.unmap()
     
    17731929        if unit is not None:
    17741930            self.set_fluxunit(unit)
    1775         self.set_freqframe(rcParams['scantable.freqframe'])
    1776 
     1931        #self.set_freqframe(rcParams['scantable.freqframe'])
     1932
  • branches/alma/src/MathUtils.cpp

    r1400 r1446  
    3737#include <casa/BasicSL/String.h>
    3838#include <scimath/Mathematics/MedianSlider.h>
     39#include <casa/Exceptions/Error.h>
    3940
    4041#include "MathUtils.h"
     
    4748   String str(which);
    4849   str.upcase();
    49    if (str.contains(String("MIN"))) {
     50   if (str.matches(String("MIN"))) {
    5051      return min(data);
    51    } else if (str.contains(String("MAX"))) {
     52   } else if (str.matches(String("MAX"))) {
    5253      return max(data);
    53    } else if (str.contains(String("SUMSQ"))) {
     54   } else if (str.matches(String("SUMSQ"))) {
    5455      return sumsquares(data);
    55    } else if (str.contains(String("SUM"))) {
     56   } else if (str.matches(String("SUM"))) {
    5657      return sum(data);
    57    } else if (str.contains(String("MEAN"))) {
     58   } else if (str.matches(String("MEAN"))) {
    5859      return mean(data);
    59    } else if (str.contains(String("VAR"))) {
     60   } else if (str.matches(String("VAR"))) {
    6061      return variance(data);
    61    } else if (str.contains(String("STDDEV"))) {
     62      } else if (str.matches(String("STDDEV"))) {
    6263      return stddev(data);
    63    } else if (str.contains(String("AVDEV"))) {
     64   } else if (str.matches(String("AVDEV"))) {
    6465      return avdev(data);
    65    } else if (str.contains(String("RMS"))) {
     66   } else if (str.matches(String("RMS"))) {
    6667      uInt n = data.nelementsValid();
    6768      return sqrt(sumsquares(data)/n);
    68    } else if (str.contains(String("MED"))) {
     69   } else if (str.matches(String("MEDIAN"))) {
    6970      return median(data);
    70    }
     71   } else {
     72      String msg = str + " is not a valid type of statistics";
     73      throw(AipsError(msg));
     74   }
    7175   return 0.0;
    7276}
  • branches/alma/src/MathUtils.h

    r1373 r1446  
    8282
    8383/**
    84  * Convert casa implementations to stl
    85  * @param in casa string
    86  * @return a std vector of std strings
     84 * Convert a std::vector of std::string
     85 * to a casa::Vector casa::String
     86 * @param in
     87 * @return
    8788 */
    8889std::vector<std::string> tovectorstring(const casa::Vector<casa::String>& in);
    8990
    9091/**
    91  * convert stl implementations to casa versions
     92 * Convert a casa::Vector of casa::String
     93 * to a stl std::vector of stl std::string
    9294 * @param in
    9395 * @return
  • branches/alma/src/RowAccumulator.cpp

    r1352 r1446  
    152152  userMask_ = m;
    153153}
     154
     155// Added by TT  check the state of RowAccumulator
     156casa::Bool RowAccumulator::state() const
     157{
     158  return initialized_;
     159}
  • branches/alma/src/RowAccumulator.h

    r1353 r1446  
    8585    */
    8686  void reset();
     87  /**
     88    * check the initialization state
     89    */
     90  casa::Bool state() const;
    8791
    8892private:
  • branches/alma/src/STAsciiWriter.cpp

    r1375 r1446  
    122122    String wcs = stable.frequencies().print(rec.asuInt("FREQ_ID"), True);
    123123    addLine(of, "WCS", wcs);
    124     addLine(of, "Rest Freq.",
    125             stable.molecules().getRestFrequency(rec.asuInt("MOLECULE_ID") ));
     124    std::vector<double> restfreqs= stable.molecules().getRestFrequency(rec.asuInt("MOLECULE_ID"));
     125    int nf = restfreqs.size();
     126    //addLine(of, "Rest Freq.",
     127    //        stable.molecules().getRestFrequency(rec.asuInt("MOLECULE_ID") ));
     128    addLine(of, "Rest Freq.", restfreqs[0]);
     129    for ( unsigned int i=1; i<nf; ++i) {
     130        addLine(of, " ", restfreqs[i]);
     131    }
    126132    of << setfill('#') << setw(70) << "" << setfill(' ') << endl;
    127133
  • branches/alma/src/STFiller.cpp

    r1387 r1446  
    205205    }
    206206  }
     207  String freqFrame = header_->freqref;
     208  //translate frequency reference frame back to
     209  //MS style (as PKSMS2reader converts the original frame
     210  //in FITS standard style)
     211  if (freqFrame == "TOPOCENT") {
     212    freqFrame = "TOPO";
     213  } else if (freqFrame == "GEOCENER") {
     214    freqFrame = "GEO";
     215  } else if (freqFrame == "BARYCENT") {
     216    freqFrame = "BARY";
     217  } else if (freqFrame == "GALACTOC") {
     218    freqFrame = "GALACTO";
     219  } else if (freqFrame == "LOCALGRP") {
     220    freqFrame = "LGROUP";
     221  } else if (freqFrame == "CMBDIPOL") {
     222    freqFrame = "CMB";
     223  } else if (freqFrame == "SOURCE") {
     224    freqFrame = "REST";
     225  }
     226  table_->frequencies().setFrame(freqFrame);
    207227     
    208228}
     
    222242  Float  azimuth, elevation, focusAxi, focusRot, focusTan,
    223243    humidity, parAngle, pressure, temperature, windAz, windSpeed;
    224   Double bandwidth, freqInc, interval, mjd, refFreq, restFreq, srcVel;
     244  Double bandwidth, freqInc, interval, mjd, refFreq, srcVel;
    225245  String          fieldName, srcName, tcalTime, obsType;
    226246  Vector<Float>   calFctr, sigma, tcal, tsys;
    227247  Matrix<Float>   baseLin, baseSub;
    228   Vector<Double>  direction(2), scanRate(2), srcDir(2), srcPM(2);
     248  Vector<Double>  direction(2), scanRate(2), srcDir(2), srcPM(2), restFreq;
    229249  Matrix<Float>   spectra;
    230250  Matrix<uChar>   flagtra;
  • branches/alma/src/STFrequencies.cpp

    r1375 r1446  
    128128}
    129129
     130void STFrequencies::setEntry( Double refpix, Double refval, Double inc, uInt id )
     131{
     132  Table t = table_(table_.col("ID") == Int(id) );
     133  if (t.nrow() == 0 ) {
     134    throw(AipsError("STFrequencies::getEntry - freqID out of range"));
     135  }
     136  for ( uInt i = 0 ; i < table_.nrow() ; i++ ) {
     137    uInt fid ;
     138    idCol_.get( i, fid ) ;
     139    if ( fid == id ) {
     140      refpixCol_.put( i, refpix ) ;
     141      refvalCol_.put( i, refval ) ;
     142      incrCol_.put( i, inc ) ;
     143    }
     144  }
     145}
     146
    130147SpectralCoordinate STFrequencies::getSpectralCoordinate( uInt id ) const
    131148{
     
    140157  // get first row - there should only be one matching id
    141158  const TableRecord& rec = row.get(0);
     159
    142160  return SpectralCoordinate( getFrame(true), rec.asDouble("REFVAL"),
    143161                             rec.asDouble("INCREMENT"),
     
    145163}
    146164
     165/**
    147166SpectralCoordinate
    148   STFrequencies::getSpectralCoordinate( const MDirection& md,
    149                                         const MPosition& mp,
    150                                         const MEpoch& me,
    151                                         Double restfreq, uInt id ) const
     167  asap::STFrequencies::getSpectralCoordinate( const MDirection& md,
     168                                              const MPosition& mp,
     169                                              const MEpoch& me,
     170                                              Double restfreq, uInt id ) const
     171**/
     172SpectralCoordinate
     173  asap::STFrequencies::getSpectralCoordinate( const MDirection& md,
     174                                              const MPosition& mp,
     175                                              const MEpoch& me,
     176                                              Vector<Double> restfreq, uInt id ) const
    152177{
    153178  SpectralCoordinate spc = getSpectralCoordinate(id);
    154   spc.setRestFrequency(restfreq, True);
     179  //spc.setRestFrequency(restfreq, True);
     180  // for now just use the first rest frequency
     181  if (restfreq.nelements()==0 ) {
     182    restfreq.resize(1);
     183    restfreq[0] = 0;
     184  }
     185  spc.setRestFrequency(restfreq[0], True);
    155186  if ( !spc.setReferenceConversion(getFrame(), me, mp, md) ) {
    156187    throw(AipsError("Couldn't convert frequency frame."));
  • branches/alma/src/STFrequencies.h

    r1375 r1446  
    5959                 casa::Double& inc, casa::uInt id );
    6060
     61  /***
     62   * Set the frequency values for a specific id via references
     63   * @param refpix the reference pixel
     64   * @param refval the reference value
     65   * @param inc    the increment
     66   * @param id     the identifier
     67   *
     68   * 17/09/2008 Takeshi Nakazato
     69   ***/
     70  void setEntry( casa::Double refpix, casa::Double refval,
     71                 casa::Double inc, casa::uInt id ) ;
     72
    6173
    6274  bool conformant(const STFrequencies& other) const;
     
    6981  casa::SpectralCoordinate getSpectralCoordinate( casa::uInt freqID ) const;
    7082
     83  /**
    7184  casa::SpectralCoordinate getSpectralCoordinate( const casa::MDirection& md,
    7285                                                  const casa::MPosition& mp,
     
    7588                                                  casa::uInt freqID
    7689                                                  ) const;
     90  **/
     91  casa::SpectralCoordinate getSpectralCoordinate( const casa::MDirection& md,
     92                                                  const casa::MPosition& mp,
     93                                                  const casa::MEpoch& me,
     94                                                  casa::Vector<casa::Double> restfreq,
     95                                                  casa::uInt freqID
     96                                                  ) const;
    7797
    7898  /**
  • branches/alma/src/STMath.cpp

    r1387 r1446  
    191191    }
    192192    //write out
    193     Vector<uChar> flg(msk.shape());
    194     convertArray(flg, !msk);
    195     flagColOut.put(i, flg);
    196     specColOut.put(i, acc.getSpectrum());
    197     tsysColOut.put(i, acc.getTsys());
    198     intColOut.put(i, acc.getInterval());
    199     mjdColOut.put(i, acc.getTime());
    200     // we should only have one cycle now -> reset it to be 0
    201     // frequency switched data has different CYCLENO for different IFNO
    202     // which requires resetting this value
    203     cycColOut.put(i, uInt(0));
     193    if (acc.state()) {
     194      Vector<uChar> flg(msk.shape());
     195      convertArray(flg, !msk);
     196      flagColOut.put(i, flg);
     197      specColOut.put(i, acc.getSpectrum());
     198      tsysColOut.put(i, acc.getTsys());
     199      intColOut.put(i, acc.getInterval());
     200      mjdColOut.put(i, acc.getTime());
     201      // we should only have one cycle now -> reset it to be 0
     202      // frequency switched data has different CYCLENO for different IFNO
     203      // which requires resetting this value
     204      cycColOut.put(i, uInt(0));
     205    } else {
     206      ostringstream oss;
     207      oss << "For output row="<<i<<", all input rows of data are flagged. no averaging" << endl;
     208      pushLog(String(oss));
     209    }
    204210    acc.reset();
    205211  }
     
    947953{
    948954
    949 
     955 
    950956  STSelector sel;
    951957  CountedPtr< Scantable > ws = getScantable(s, false);
    952958  CountedPtr< Scantable > sig, sigwcal, ref, refwcal;
    953   CountedPtr< Scantable > calsig, calref, out;
     959  CountedPtr< Scantable > calsig, calref, out, out1, out2;
     960  Bool nofold=False;
    954961
    955962  //split the data
     
    973980  calref = dototalpower(refwcal, ref, tcal=tcal);
    974981
    975   out=dosigref(calsig,calref,smoothref,tsysv,tau);
    976 
     982  out1=dosigref(calsig,calref,smoothref,tsysv,tau);
     983  out2=dosigref(calref,calsig,smoothref,tsysv,tau);
     984
     985  Table& tabout1=out1->table();
     986  Table& tabout2=out2->table();
     987  ROScalarColumn<uInt> freqidCol1(tabout1, "FREQ_ID");
     988  ROScalarColumn<uInt> freqidCol2(tabout2, "FREQ_ID");
     989  ROArrayColumn<Float> specCol(tabout2, "SPECTRA");
     990  Vector<Float> spec; specCol.get(0, spec);
     991  uInt nchan = spec.nelements();
     992  uInt freqid1; freqidCol1.get(0,freqid1);
     993  uInt freqid2; freqidCol2.get(0,freqid2);
     994  Double rp1, rp2, rv1, rv2, inc1, inc2;
     995  out1->frequencies().getEntry(rp1, rv1, inc1, freqid1);
     996  out2->frequencies().getEntry(rp2, rv2, inc2, freqid2);
     997  if (rp1==rp2) {
     998    Double foffset = rv1 - rv2;
     999    uInt choffset = static_cast<uInt>(foffset/abs(inc2));
     1000    if (choffset >= nchan) {
     1001      cerr<<"out-band frequency switching, no folding"<<endl;
     1002      nofold = True;
     1003    }
     1004  }
     1005 
     1006  if (nofold) {
     1007    std::vector< CountedPtr< Scantable > > tabs;
     1008    tabs.push_back(out1);
     1009    tabs.push_back(out2);
     1010    out = merge(tabs);
     1011  }
     1012  else { //folding is not implemented yet
     1013    out = out1;
     1014  }
     1015   
    9771016  return out;
    9781017}
    979 
    9801018
    9811019CountedPtr< Scantable > STMath::freqSwitch( const CountedPtr< Scantable >& in )
     
    15651603          id = out->frequencies().addEntry(rp, rv, inc);
    15661604          freqidcol.put(k,id);
    1567           String name,fname;Double rf;
     1605          //String name,fname;Double rf;
     1606          Vector<String> name,fname;Vector<Double> rf;
    15681607          (*it)->molecules().getEntry(rf, name, fname, rec.asuInt("MOLECULE_ID"));
    15691608          id = out->molecules().addEntry(rf, name, fname);
     
    19812020  return out;
    19822021}
     2022
     2023// Averaging spectra with different channel/resolution
     2024CountedPtr<Scantable>
     2025STMath::new_average( const std::vector<CountedPtr<Scantable> >& in,
     2026                     const bool& compel,
     2027                     const std::vector<bool>& mask,
     2028                     const std::string& weight,
     2029                     const std::string& avmode )
     2030  throw ( casa::AipsError )
     2031{
     2032  if ( avmode == "SCAN" && in.size() != 1 )
     2033    throw(AipsError("Can't perform 'SCAN' averaging on multiple tables.\n"
     2034                    "Use merge first."));
     2035 
     2036  CountedPtr<Scantable> out ;     // processed result
     2037  if ( compel ) {
     2038    std::vector< CountedPtr<Scantable> > newin ; // input for average process
     2039    uInt insize = in.size() ;    // number of scantables
     2040
     2041    // TEST: do normal average in each table before IF grouping
     2042    cout << "Do preliminary averaging" << endl ;
     2043    vector< CountedPtr<Scantable> > tmpin( insize ) ;
     2044    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
     2045      vector< CountedPtr<Scantable> > v( 1, in[itable] ) ;
     2046      tmpin[itable] = average( v, mask, weight, avmode ) ;
     2047    }
     2048
     2049    // warning
     2050    cout << "Average spectra with different spectral resolution" << endl ;
     2051    cout << endl ;
     2052
     2053    // temporarily set coordinfo
     2054    vector<string> oldinfo( insize ) ;
     2055    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
     2056      vector<string> coordinfo = in[itable]->getCoordInfo() ;
     2057      oldinfo[itable] = coordinfo[0] ;
     2058      coordinfo[0] = "Hz" ;
     2059      tmpin[itable]->setCoordInfo( coordinfo ) ;
     2060    }
     2061
     2062    // columns
     2063    ScalarColumn<uInt> freqIDCol ;
     2064    ScalarColumn<uInt> ifnoCol ;
     2065    ScalarColumn<uInt> scannoCol ;
     2066
     2067
     2068    // check IF frequency coverage
     2069    //cout << "Check IF settings in each table" << endl ;
     2070    vector< vector<uInt> > freqid( insize );
     2071    vector< vector<double> > iffreq( insize ) ;
     2072    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
     2073      uInt rows = tmpin[itable]->nrow() ;
     2074      uInt freqnrows = tmpin[itable]->frequencies().table().nrow() ;
     2075      for ( uInt irow = 0 ; irow < rows ; irow++ ) {
     2076        if ( freqid[itable].size() == freqnrows ) {
     2077          break ;
     2078        }
     2079        else {
     2080          freqIDCol.attach( tmpin[itable]->table(), "FREQ_ID" ) ;
     2081          ifnoCol.attach( tmpin[itable]->table(), "IFNO" ) ;
     2082          uInt id = freqIDCol( irow ) ;
     2083          if ( freqid[itable].size() == 0 || count( freqid[itable].begin(), freqid[itable].end(), id ) == 0 ) {
     2084            //cout << "itable = " << itable << ": IF " << id << " is included in the list" << endl ;
     2085            vector<double> abcissa = tmpin[itable]->getAbcissa( irow ) ;
     2086            freqid[itable].push_back( id ) ;
     2087            iffreq[itable].push_back( abcissa[0] - 0.5 * ( abcissa[1] - abcissa[0] ) ) ;
     2088            iffreq[itable].push_back( abcissa[abcissa.size()-1] + 0.5 * ( abcissa[1] - abcissa[0] ) ) ;
     2089          }
     2090        }
     2091      }
     2092    }
     2093
     2094    // debug
     2095    //cout << "IF settings summary:" << endl ;
     2096    //for ( uInt i = 0 ; i < freqid.size() ; i++ ) {
     2097    //cout << "   Table" << i << endl ;
     2098    //for ( uInt j = 0 ; j < freqid[i].size() ; j++ ) {
     2099    //cout << "      id = " << freqid[i][j] << " (min,max) = (" << iffreq[i][2*j] << "," << iffreq[i][2*j+1] << ")" << endl ;
     2100    //}
     2101    //}
     2102    //cout << endl ;
     2103
     2104    // IF grouping based on their frequency coverage
     2105    //cout << "IF grouping based on their frequency coverage" << endl ;
     2106
     2107    // IF group
     2108    // ifgrp[numgrp][nummember*2]
     2109    // ifgrp = [table00, freqrow00, table01, freqrow01, ...],
     2110    //         [table11, freqrow11, table11, freqrow11, ...],
     2111    //         ...
     2112    // ifgfreq[numgrp*2]
     2113    // ifgfreq = [min0, max0, min1, max1, ...]
     2114    vector< vector<uInt> > ifgrp ;
     2115    vector<double> ifgfreq ;
     2116
     2117    // parameter for IF grouping
     2118    // groupmode = OR    retrieve all region
     2119    //             AND   only retrieve overlaped region
     2120    //string groupmode = "AND" ;
     2121    string groupmode = "OR" ;
     2122    uInt sizecr = 0 ;
     2123    if ( groupmode == "AND" )
     2124      sizecr = 2 ;
     2125    else if ( groupmode == "OR" )
     2126      sizecr = 0 ;
     2127
     2128    vector<double> sortedfreq ;
     2129    for ( uInt i = 0 ; i < iffreq.size() ; i++ ) {
     2130      for ( uInt j = 0 ; j < iffreq[i].size() ; j++ ) {
     2131        if ( count( sortedfreq.begin(), sortedfreq.end(), iffreq[i][j] ) == 0 )
     2132          sortedfreq.push_back( iffreq[i][j] ) ;
     2133      }
     2134    }
     2135    sort( sortedfreq.begin(), sortedfreq.end() ) ;
     2136    for ( vector<double>::iterator i = sortedfreq.begin() ; i != sortedfreq.end()-1 ; i++ ) {
     2137      ifgfreq.push_back( *i ) ;
     2138      ifgfreq.push_back( *(i+1) ) ;
     2139    }
     2140    ifgrp.resize( ifgfreq.size()/2 ) ;
     2141    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
     2142      for ( uInt iif = 0 ; iif < freqid[itable].size() ; iif++ ) {
     2143        double range0 = iffreq[itable][2*iif] ;
     2144        double range1 = iffreq[itable][2*iif+1] ;
     2145        for ( uInt j = 0 ; j < ifgrp.size() ; j++ ) {
     2146          double fmin = max( range0, ifgfreq[2*j] ) ;
     2147          double fmax = min( range1, ifgfreq[2*j+1] ) ;
     2148          if ( fmin < fmax ) {
     2149            ifgrp[j].push_back( itable ) ;
     2150            ifgrp[j].push_back( freqid[itable][iif] ) ;
     2151          }
     2152        }
     2153      }
     2154    }
     2155    vector< vector<uInt> >::iterator fiter = ifgrp.begin() ;
     2156    vector<double>::iterator giter = ifgfreq.begin() ;
     2157    while( fiter != ifgrp.end() ) {
     2158      if ( fiter->size() <= sizecr ) {
     2159        fiter = ifgrp.erase( fiter ) ;
     2160        giter = ifgfreq.erase( giter ) ;
     2161        giter = ifgfreq.erase( giter ) ;
     2162      }
     2163      else {
     2164        fiter++ ;
     2165        advance( giter, 2 ) ;
     2166      }
     2167    }
     2168
     2169    // freqgrp[numgrp][nummember]
     2170    // freqgrp = [ifgrp00, ifgrp01, ifgrp02, ...],
     2171    //           [ifgrp10, ifgrp11, ifgrp12, ...],
     2172    //           ...
     2173    // freqrange[numgrp*2]
     2174    // freqrange = [min0, max0, min1, max1, ...]
     2175    vector< vector<uInt> > freqgrp ;
     2176    double freqrange = 0.0 ;
     2177    uInt grpnum = 0 ;
     2178    for ( uInt i = 0 ; i < ifgrp.size() ; i++ ) {
     2179      // Assumed that ifgfreq was sorted
     2180      if ( grpnum != 0 && freqrange == ifgfreq[2*i] ) {
     2181        freqgrp[grpnum-1].push_back( i ) ;
     2182      }
     2183      else {
     2184        vector<uInt> grp0( 1, i ) ;
     2185        freqgrp.push_back( grp0 ) ;
     2186        grpnum++ ;
     2187      }
     2188      freqrange = ifgfreq[2*i+1] ;
     2189    }
     2190       
     2191
     2192    // print IF groups
     2193    cout << "IF Group summary: " << endl ;
     2194    cout << "   GROUP_ID [FREQ_MIN, FREQ_MAX]: (TABLE_ID, FREQ_ID)" << endl ;
     2195    for ( uInt i = 0 ; i < ifgrp.size() ; i++ ) {
     2196      cout << "   GROUP " << setw( 2 ) << i << " [" << ifgfreq[2*i] << "," << ifgfreq[2*i+1] << "]: " ;
     2197      for ( uInt j = 0 ; j < ifgrp[i].size()/2 ; j++ ) {
     2198        cout << "(" << ifgrp[i][2*j] << "," << ifgrp[i][2*j+1] << ") " ;
     2199      }
     2200      cout << endl ;
     2201    }
     2202    cout << endl ;
     2203   
     2204    // print frequency group
     2205    cout << "Frequency Group summary: " << endl ;
     2206    cout << "   GROUP_ID [FREQ_MIN, FREQ_MAX]: IF_GROUP_ID" << endl ;
     2207    for ( uInt i = 0 ; i < freqgrp.size() ; i++ ) {
     2208      cout << "   GROUP " << setw( 2 ) << i << " [" << ifgfreq[2*freqgrp[i][0]] << "," << ifgfreq[2*freqgrp[i][freqgrp[i].size()-1]+1] << "]: " ;
     2209      for ( uInt j = 0 ; j < freqgrp[i].size() ; j++ ) {
     2210        cout << freqgrp[i][j] << " " ;
     2211      }
     2212      cout << endl ;
     2213    }
     2214    cout << endl ;
     2215
     2216    // membership check
     2217    // groups[numtable][numIF][nummembership]
     2218    vector< vector< vector<uInt> > > groups( insize ) ;
     2219    for ( uInt i = 0 ; i < insize ; i++ ) {
     2220      groups[i].resize( freqid[i].size() ) ;
     2221    }
     2222    for ( uInt igrp = 0 ; igrp < ifgrp.size() ; igrp++ ) {
     2223      for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
     2224        uInt tableid = ifgrp[igrp][2*imem] ;
     2225        vector<uInt>::iterator iter = find( freqid[tableid].begin(), freqid[tableid].end(), ifgrp[igrp][2*imem+1] ) ;
     2226        if ( iter != freqid[tableid].end() ) {
     2227          uInt rowid = distance( freqid[tableid].begin(), iter ) ;
     2228          groups[tableid][rowid].push_back( igrp ) ;
     2229        }
     2230      }
     2231    }
     2232
     2233    // print membership
     2234    //for ( uInt i = 0 ; i < insize ; i++ ) {
     2235    //cout << "Table " << i << endl ;
     2236    //for ( uInt j = 0 ; j < groups[i].size() ; j++ ) {
     2237    //cout << "   FREQ_ID " <<  setw( 2 ) << freqid[i][j] << ": " ;
     2238    //for ( uInt k = 0 ; k < groups[i][j].size() ; k++ ) {
     2239    //cout << setw( 2 ) << groups[i][j][k] << " " ;
     2240    //}
     2241    //cout << endl ;
     2242    //}
     2243    //}
     2244
     2245    // set back coordinfo
     2246    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
     2247      vector<string> coordinfo = tmpin[itable]->getCoordInfo() ;
     2248      coordinfo[0] = oldinfo[itable] ;
     2249      tmpin[itable]->setCoordInfo( coordinfo ) ;
     2250    }
     2251
     2252    // Create additional table if needed
     2253    bool oldInsitu = insitu_ ;
     2254    setInsitu( false ) ;
     2255    vector< vector<uInt> > addrow( insize ) ;
     2256    vector<uInt> addtable( insize, 0 ) ;
     2257    vector<uInt> newtableids( insize ) ;
     2258    vector<uInt> newifids( insize, 0 ) ;
     2259    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
     2260      //cout << "Table " << setw(2) << itable << ": " ;
     2261      for ( uInt ifrow = 0 ; ifrow < groups[itable].size() ; ifrow++ ) {
     2262        addrow[itable].push_back( groups[itable][ifrow].size()-1 ) ;
     2263        //cout << addrow[itable][ifrow] << " " ;
     2264      }
     2265      addtable[itable] = *max_element( addrow[itable].begin(), addrow[itable].end() ) ;
     2266      //cout << "(" << addtable[itable] << ")" << endl ;
     2267    }
     2268    newin.resize( insize ) ;
     2269    copy( tmpin.begin(), tmpin.end(), newin.begin() ) ;
     2270    for ( uInt i = 0 ; i < insize ; i++ ) {
     2271      newtableids[i] = i ;
     2272    }
     2273    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
     2274      for ( uInt iadd = 0 ; iadd < addtable[itable] ; iadd++ ) {
     2275        CountedPtr<Scantable> add = getScantable( newin[itable], false ) ;
     2276        vector<int> freqidlist ;
     2277        for ( uInt i = 0 ; i < groups[itable].size() ; i++ ) {
     2278          if ( groups[itable][i].size() > iadd + 1 ) {
     2279            freqidlist.push_back( freqid[itable][i] ) ;
     2280          }
     2281        }
     2282        stringstream taqlstream ;
     2283        taqlstream << "SELECT FROM $1 WHERE FREQ_ID IN [" ;
     2284        for ( uInt i = 0 ; i < freqidlist.size() ; i++ ) {
     2285          taqlstream << i ;
     2286          if ( i < freqidlist.size() - 1 )
     2287            taqlstream << "," ;
     2288          else
     2289            taqlstream << "]" ;
     2290        }
     2291        string taql = taqlstream.str() ;
     2292        //cout << "taql = " << taql << endl ;
     2293        STSelector selector = STSelector() ;
     2294        selector.setTaQL( taql ) ;
     2295        add->setSelection( selector ) ;
     2296        newin.push_back( add ) ;
     2297        newtableids.push_back( itable ) ;
     2298        newifids.push_back( iadd + 1 ) ;
     2299      }
     2300    }
     2301
     2302    // udpate ifgrp
     2303    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
     2304      for ( uInt iadd = 0 ; iadd < addtable[itable] ; iadd++ ) {
     2305        for ( uInt ifrow = 0 ; ifrow < groups[itable].size() ; ifrow++ ) {
     2306          if ( groups[itable][ifrow].size() > iadd + 1 ) {
     2307            uInt igrp = groups[itable][ifrow][iadd+1] ;
     2308            for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
     2309              if ( ifgrp[igrp][2*imem] == newtableids[iadd+insize] && ifgrp[igrp][2*imem+1] == freqid[newtableids[iadd+insize]][ifrow] ) {
     2310                ifgrp[igrp][2*imem] = insize + iadd ;
     2311              }
     2312            }
     2313          }
     2314        }
     2315      }
     2316    }
     2317
     2318    // print IF groups again for debug
     2319    //cout << "IF Group summary: " << endl ;
     2320    //cout << "   GROUP_ID [FREQ_MIN, FREQ_MAX]: (TABLE_ID, FREQ_ID)" << endl ;
     2321    //for ( uInt i = 0 ; i < ifgrp.size() ; i++ ) {
     2322    //cout << "   GROUP " << setw( 2 ) << i << " [" << ifgfreq[2*i] << "," << ifgfreq[2*i+1] << "]: " ;
     2323    //for ( uInt j = 0 ; j < ifgrp[i].size()/2 ; j++ ) {
     2324    //cout << "(" << ifgrp[i][2*j] << "," << ifgrp[i][2*j+1] << ") " ;
     2325    //}
     2326    //cout << endl ;
     2327    //}
     2328    //cout << endl ;
     2329
     2330    // reset SCANNO and IF: IF number is reset by the result of sortation
     2331    cout << "All scan number is set to 0" << endl ;
     2332    //cout << "All IF number is set to IF group index" << endl ;
     2333    cout << endl ;
     2334    insize = newin.size() ;
     2335    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
     2336      uInt rows = newin[itable]->nrow() ;
     2337      Table &tmpt = newin[itable]->table() ;
     2338      freqIDCol.attach( tmpt, "FREQ_ID" ) ;
     2339      scannoCol.attach( tmpt, "SCANNO" ) ;
     2340      ifnoCol.attach( tmpt, "IFNO" ) ;
     2341      for ( uInt irow=0 ; irow < rows ; irow++ ) {
     2342        scannoCol.put( irow, 0 ) ;
     2343        uInt freqID = freqIDCol( irow ) ;
     2344        vector<uInt>::iterator iter = find( freqid[newtableids[itable]].begin(), freqid[newtableids[itable]].end(), freqID ) ;
     2345        if ( iter != freqid[newtableids[itable]].end() ) {
     2346          uInt index = distance( freqid[newtableids[itable]].begin(), iter ) ;
     2347          ifnoCol.put( irow, groups[newtableids[itable]][index][newifids[itable]] ) ;
     2348        }
     2349        else {
     2350          throw(AipsError("IF grouping was wrong in additional tables.")) ;
     2351        }
     2352      }
     2353    }
     2354    oldinfo.resize( insize ) ;
     2355    setInsitu( oldInsitu ) ;
     2356
     2357    // temporarily set coordinfo
     2358    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
     2359      vector<string> coordinfo = newin[itable]->getCoordInfo() ;
     2360      oldinfo[itable] = coordinfo[0] ;
     2361      coordinfo[0] = "Hz" ;
     2362      newin[itable]->setCoordInfo( coordinfo ) ;
     2363    }
     2364
     2365    // save column values in the vector
     2366    vector< vector<uInt> > freqTableIdVec( insize ) ;
     2367    vector< vector<uInt> > freqIdVec( insize ) ;
     2368    vector< vector<uInt> > ifNoVec( insize ) ;
     2369    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
     2370      ScalarColumn<uInt> freqIDs ;
     2371      freqIDs.attach( newin[itable]->frequencies().table(), "ID" ) ;
     2372      ifnoCol.attach( newin[itable]->table(), "IFNO" ) ;
     2373      freqIDCol.attach( newin[itable]->table(), "FREQ_ID" ) ;
     2374      for ( uInt irow = 0 ; irow < newin[itable]->frequencies().table().nrow() ; irow++ ) {
     2375        freqTableIdVec[itable].push_back( freqIDs( irow ) ) ;
     2376      }
     2377      for ( uInt irow = 0 ; irow < newin[itable]->table().nrow() ; irow++ ) {
     2378        freqIdVec[itable].push_back( freqIDCol( irow ) ) ;
     2379        ifNoVec[itable].push_back( ifnoCol( irow ) ) ;
     2380      }
     2381    }
     2382
     2383    // reset spectra and flagtra: pick up common part of frequency coverage
     2384    //cout << "Pick common frequency range and align resolution" << endl ;
     2385    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
     2386      uInt rows = newin[itable]->nrow() ;
     2387      int nminchan = -1 ;
     2388      int nmaxchan = -1 ;
     2389      vector<uInt> freqIdUpdate ;
     2390      for ( uInt irow = 0 ; irow < rows ; irow++ ) {
     2391        uInt ifno = ifNoVec[itable][irow] ;  // IFNO is reset by group index
     2392        double minfreq = ifgfreq[2*ifno] ;
     2393        double maxfreq = ifgfreq[2*ifno+1] ;
     2394        //cout << "frequency range: [" << minfreq << "," << maxfreq << "]" << endl ;
     2395        vector<double> abcissa = newin[itable]->getAbcissa( irow ) ;
     2396        int nchan = abcissa.size() ;
     2397        double resol = abcissa[1] - abcissa[0] ;
     2398        //cout << "abcissa range  : [" << abcissa[0] << "," << abcissa[nchan-1] << "]" << endl ;
     2399        if ( minfreq <= abcissa[0] )
     2400          nminchan = 0 ;
     2401        else {
     2402          //double cfreq = ( minfreq - abcissa[0] ) / resol ;
     2403          double cfreq = ( minfreq - abcissa[0] + 0.5 * resol ) / resol ;
     2404          nminchan = int(cfreq) + ( ( cfreq - int(cfreq) <= 0.5 ) ? 0 : 1 ) ;
     2405        }
     2406        if ( maxfreq >= abcissa[abcissa.size()-1] )
     2407          nmaxchan = abcissa.size() - 1 ;
     2408        else {
     2409          //double cfreq = ( abcissa[abcissa.size()-1] - maxfreq ) / resol ;
     2410          double cfreq = ( abcissa[abcissa.size()-1] - maxfreq + 0.5 * resol ) / resol ;
     2411          nmaxchan = abcissa.size() - 1 - int(cfreq) - ( ( cfreq - int(cfreq) >= 0.5 ) ? 1 : 0 ) ;
     2412        }
     2413        //cout << "channel range (" << irow << "): [" << nminchan << "," << nmaxchan << "]" << endl ;
     2414        if ( nmaxchan > nminchan ) {
     2415          newin[itable]->reshapeSpectrum( nminchan, nmaxchan, irow ) ;
     2416          int newchan = nmaxchan - nminchan + 1 ;
     2417          if ( count( freqIdUpdate.begin(), freqIdUpdate.end(), freqIdVec[itable][irow] ) == 0 && newchan < nchan )
     2418            freqIdUpdate.push_back( freqIdVec[itable][irow] ) ;
     2419        }
     2420        else {
     2421          throw(AipsError("Failed to pick up common part of frequency range.")) ;
     2422        }
     2423      }
     2424      for ( uInt i = 0 ; i < freqIdUpdate.size() ; i++ ) {
     2425        uInt freqId = freqIdUpdate[i] ;
     2426        Double refpix ;
     2427        Double refval ;
     2428        Double increment ;
     2429       
     2430        // update row
     2431        newin[itable]->frequencies().getEntry( refpix, refval, increment, freqId ) ;
     2432        refval = refval - ( refpix - nminchan ) * increment ;
     2433        refpix = 0 ;
     2434        newin[itable]->frequencies().setEntry( refpix, refval, increment, freqId ) ;
     2435      }   
     2436    }
     2437
     2438   
     2439    // reset spectra and flagtra: align spectral resolution
     2440    //cout << "Align spectral resolution" << endl ;
     2441    vector<double> gmaxdnu( freqgrp.size(), 0.0 ) ;
     2442    vector<uInt> gmemid( freqgrp.size(), 0 ) ;
     2443    for ( uInt igrp = 0 ; igrp < ifgrp.size() ; igrp++ ) {
     2444      double maxdnu = 0.0 ;       // maximum (coarsest) frequency resolution
     2445      int minchan = INT_MAX ;     // minimum channel number
     2446      Double refpixref = -1 ;     // reference of 'reference pixel'
     2447      Double refvalref = -1 ;     // reference of 'reference frequency'
     2448      Double refinc = -1 ;        // reference frequency resolution
     2449      uInt refreqid ;
     2450      uInt reftable = INT_MAX;
     2451      // process only if group member > 1
     2452      if ( ifgrp[igrp].size() > 2 ) {
     2453        // find minchan and maxdnu in each group
     2454        for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
     2455          uInt tableid = ifgrp[igrp][2*imem] ;
     2456          uInt rowid = ifgrp[igrp][2*imem+1] ;
     2457          vector<uInt>::iterator iter = find( freqIdVec[tableid].begin(), freqIdVec[tableid].end(), rowid ) ;
     2458          if ( iter != freqIdVec[tableid].end() ) {
     2459            uInt index = distance( freqIdVec[tableid].begin(), iter ) ;
     2460            vector<double> abcissa = newin[tableid]->getAbcissa( index ) ;
     2461            int nchan = abcissa.size() ;
     2462            double dnu = abcissa[1] - abcissa[0] ;
     2463            //cout << "GROUP " << igrp << " (" << tableid << "," << rowid << "): nchan = " << nchan << " (minchan = " << minchan << ")" << endl ;
     2464            if ( nchan < minchan ) {
     2465              minchan = nchan ;
     2466              maxdnu = dnu ;
     2467              newin[tableid]->frequencies().getEntry( refpixref, refvalref, refinc, rowid ) ;
     2468              refreqid = rowid ;
     2469              reftable = tableid ;
     2470            }
     2471          }
     2472        }
     2473        // regrid spectra in each group
     2474        cout << "GROUP " << igrp << endl ;
     2475        cout << "   Channel number is adjusted to " << minchan << endl ;
     2476        cout << "   Corresponding frequency resolution is " << maxdnu << "Hz" << endl ;
     2477        for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
     2478          uInt tableid = ifgrp[igrp][2*imem] ;
     2479          uInt rowid = ifgrp[igrp][2*imem+1] ;
     2480          freqIDCol.attach( newin[tableid]->table(), "FREQ_ID" ) ;
     2481          //cout << "tableid = " << tableid << " rowid = " << rowid << ": " << endl ;
     2482          //cout << "   regridChannel applied to " ;
     2483          if ( tableid != reftable )
     2484            refreqid = newin[tableid]->frequencies().addEntry( refpixref, refvalref, refinc ) ;
     2485          for ( uInt irow = 0 ; irow < newin[tableid]->table().nrow() ; irow++ ) {
     2486            uInt tfreqid = freqIdVec[tableid][irow] ;
     2487            if ( tfreqid == rowid ) {     
     2488              //cout << irow << " " ;
     2489              newin[tableid]->regridChannel( minchan, maxdnu, irow ) ;
     2490              freqIDCol.put( irow, refreqid ) ;
     2491              freqIdVec[tableid][irow] = refreqid ;
     2492            }
     2493          }
     2494          //cout << endl ;
     2495        }
     2496      }
     2497      else {
     2498        uInt tableid = ifgrp[igrp][0] ;
     2499        uInt rowid = ifgrp[igrp][1] ;
     2500        vector<uInt>::iterator iter = find( freqIdVec[tableid].begin(), freqIdVec[tableid].end(), rowid ) ;
     2501        if ( iter != freqIdVec[tableid].end() ) {
     2502          uInt index = distance( freqIdVec[tableid].begin(), iter ) ;
     2503          vector<double> abcissa = newin[tableid]->getAbcissa( index ) ;
     2504          minchan = abcissa.size() ;
     2505          maxdnu = abcissa[1] - abcissa[0] ;
     2506        }
     2507      }
     2508      for ( uInt i = 0 ; i < freqgrp.size() ; i++ ) {
     2509        if ( count( freqgrp[i].begin(), freqgrp[i].end(), igrp ) > 0 ) {
     2510          if ( maxdnu > gmaxdnu[i] ) {
     2511            gmaxdnu[i] = maxdnu ;
     2512            gmemid[i] = igrp ;
     2513          }
     2514          break ;
     2515        }
     2516      }
     2517    }
     2518    cout << endl ;
     2519
     2520    // set back coordinfo
     2521    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
     2522      vector<string> coordinfo = newin[itable]->getCoordInfo() ;
     2523      coordinfo[0] = oldinfo[itable] ;
     2524      newin[itable]->setCoordInfo( coordinfo ) ;
     2525    }     
     2526
     2527    // accumulate all rows into the first table
     2528    // NOTE: assumed in.size() = 1
     2529    vector< CountedPtr<Scantable> > tmp( 1 ) ;
     2530    if ( newin.size() == 1 )
     2531      tmp[0] = newin[0] ;
     2532    else
     2533      tmp[0] = merge( newin ) ;
     2534
     2535    //return tmp[0] ;
     2536
     2537    // average
     2538    CountedPtr<Scantable> tmpout = average( tmp, mask, weight, avmode ) ;
     2539
     2540    //return tmpout ;
     2541
     2542    // combine frequency group
     2543    cout << "Combine spectra based on frequency grouping" << endl ;
     2544    cout << "IFNO is renumbered as frequency group ID (see above)" << endl ;
     2545    vector<string> coordinfo = tmpout->getCoordInfo() ;
     2546    oldinfo[0] = coordinfo[0] ;
     2547    coordinfo[0] = "Hz" ;
     2548    tmpout->setCoordInfo( coordinfo ) ;
     2549    // create proformas of output table
     2550    stringstream taqlstream ;
     2551    taqlstream << "SELECT FROM $1 WHERE IFNO IN [" ;
     2552    for ( uInt i = 0 ; i < gmemid.size() ; i++ ) {
     2553      taqlstream << gmemid[i] ;
     2554      if ( i < gmemid.size() - 1 )
     2555        taqlstream << "," ;
     2556      else
     2557        taqlstream << "]" ;
     2558    }
     2559    string taql = taqlstream.str() ;
     2560    //cout << "taql = " << taql << endl ;
     2561    STSelector selector = STSelector() ;
     2562    selector.setTaQL( taql ) ;
     2563    oldInsitu = insitu_ ;
     2564    setInsitu( false ) ;
     2565    out = getScantable( tmpout, false ) ;
     2566    setInsitu( oldInsitu ) ;
     2567    out->setSelection( selector ) ;
     2568    // regrid rows
     2569    ifnoCol.attach( tmpout->table(), "IFNO" ) ;
     2570    for ( uInt irow = 0 ; irow < tmpout->table().nrow() ; irow++ ) {
     2571      uInt ifno = ifnoCol( irow ) ;
     2572      for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
     2573        if ( count( freqgrp[igrp].begin(), freqgrp[igrp].end(), ifno ) > 0 ) {
     2574          vector<double> abcissa = tmpout->getAbcissa( irow ) ;
     2575          double bw = ( abcissa[1] - abcissa[0] ) * abcissa.size() ;
     2576          int nchan = (int)( bw / gmaxdnu[igrp] ) ;
     2577          tmpout->regridChannel( nchan, gmaxdnu[igrp], irow ) ;
     2578          break ;
     2579        }
     2580      }
     2581    }
     2582    // combine spectra
     2583    ArrayColumn<Float> specColOut ;
     2584    specColOut.attach( out->table(), "SPECTRA" ) ;
     2585    ArrayColumn<uChar> flagColOut ;
     2586    flagColOut.attach( out->table(), "FLAGTRA" ) ;
     2587    ScalarColumn<uInt> ifnoColOut ;
     2588    ifnoColOut.attach( out->table(), "IFNO" ) ;
     2589    ScalarColumn<uInt> polnoColOut ;
     2590    polnoColOut.attach( out->table(), "POLNO" ) ;
     2591    ScalarColumn<uInt> freqidColOut ;
     2592    freqidColOut.attach( out->table(), "FREQ_ID" ) ;
     2593    Table &tab = tmpout->table() ;
     2594    TableIterator iter( tab, "POLNO" ) ;
     2595    vector< vector<uInt> > sizes( freqgrp.size() ) ;
     2596    while( !iter.pastEnd() ) {
     2597      vector< vector<Float> > specout( freqgrp.size() ) ;
     2598      vector< vector<uChar> > flagout( freqgrp.size() ) ;
     2599      ArrayColumn<Float> specCols ;
     2600      specCols.attach( iter.table(), "SPECTRA" ) ;
     2601      ArrayColumn<uChar> flagCols ;
     2602      flagCols.attach( iter.table(), "FLAGTRA" ) ;
     2603      ifnoCol.attach( iter.table(), "IFNO" ) ;
     2604      ScalarColumn<uInt> polnos ;
     2605      polnos.attach( iter.table(), "POLNO" ) ;
     2606      uInt polno = polnos( 0 ) ;
     2607      //cout << "POLNO iteration: " << polno << endl ;
     2608      for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
     2609        sizes[igrp].resize( freqgrp[igrp].size() ) ;
     2610        for ( uInt imem = 0 ; imem < freqgrp[igrp].size() ; imem++ ) {
     2611          for ( uInt irow = 0 ; irow < iter.table().nrow() ; irow++ ) {
     2612            uInt ifno = ifnoCol( irow ) ;
     2613            if ( ifno == freqgrp[igrp][imem] ) {
     2614              Vector<Float> spec = specCols( irow ) ;
     2615              Vector<uChar> flag = flagCols( irow ) ;
     2616              vector<Float> svec ;
     2617              spec.tovector( svec ) ;
     2618              vector<uChar> fvec ;
     2619              flag.tovector( fvec ) ;
     2620              //cout << "spec.size() = " << svec.size() << " fvec.size() = " << fvec.size() << endl ;
     2621              specout[igrp].insert( specout[igrp].end(), svec.begin(), svec.end() ) ;
     2622              flagout[igrp].insert( flagout[igrp].end(), fvec.begin(), fvec.end() ) ;
     2623              //cout << "specout[" << igrp << "].size() = " << specout[igrp].size() << endl ;
     2624              sizes[igrp][imem] = spec.nelements() ;
     2625            }
     2626          }
     2627        }
     2628        for ( uInt irow = 0 ; irow < out->table().nrow() ; irow++ ) {
     2629          uInt ifout = ifnoColOut( irow ) ;
     2630          uInt polout = polnoColOut( irow ) ;
     2631          if ( ifout == gmemid[igrp] && polout == polno ) {
     2632            // set SPECTRA and FRAGTRA
     2633            Vector<Float> newspec( specout[igrp] ) ;
     2634            Vector<uChar> newflag( flagout[igrp] ) ;
     2635            specColOut.put( irow, newspec ) ;
     2636            flagColOut.put( irow, newflag ) ;
     2637            // IFNO renumbering
     2638            ifnoColOut.put( irow, igrp ) ;
     2639          }
     2640        }
     2641      }
     2642      iter++ ;
     2643    }
     2644    // update FREQUENCIES subtable
     2645    vector<bool> updated( freqgrp.size(), false ) ;
     2646    for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
     2647      uInt index = 0 ;
     2648      uInt pixShift = 0 ;
     2649      while ( freqgrp[igrp][index] != gmemid[igrp] ) {
     2650        pixShift += sizes[igrp][index++] ;
     2651      }
     2652      for ( uInt irow = 0 ; irow < out->table().nrow() ; irow++ ) {
     2653        if ( ifnoColOut( irow ) == gmemid[igrp] && !updated[igrp] ) {
     2654          uInt freqidOut = freqidColOut( irow ) ;
     2655          //cout << "freqgrp " << igrp << " freqidOut = " << freqidOut << endl ;
     2656          double refpix ;
     2657          double refval ;
     2658          double increm ;
     2659          out->frequencies().getEntry( refpix, refval, increm, freqidOut ) ;
     2660          refpix += pixShift ;
     2661          out->frequencies().setEntry( refpix, refval, increm, freqidOut ) ;
     2662          updated[igrp] = true ;
     2663        }
     2664      }
     2665    }
     2666
     2667    //out = tmpout ;
     2668
     2669    coordinfo = tmpout->getCoordInfo() ;
     2670    coordinfo[0] = oldinfo[0] ;
     2671    tmpout->setCoordInfo( coordinfo ) ;
     2672  }
     2673  else {
     2674    // simple average
     2675    out =  average( in, mask, weight, avmode ) ;
     2676  }
     2677 
     2678  return out ;
     2679}
  • branches/alma/src/STMath.h

    r1388 r1446  
    6565    * average a vector of Scantables
    6666    * @param in the vector of Scantables to average
    67     * @param mask an optional mask to apply on specific weights
     67    * @param an optional mask to apply on specific weights
    6868    * @param weight weighting scheme
    6969    * @param avmode the mode ov averaging. Per "SCAN" or "ALL".
     
    134134  /**
    135135    * Calibrate total power scans (translated from GBTIDL)
    136     * @param calon uncalibrated Scantable with CAL noise signal
     136    * @param calon uncalibrated Scantable with CAL noise signal 
    137137    * @param caloff uncalibrated Scantable with no CAL signal
    138138    * @param tcal optional scalar Tcal, CAL temperature (K)
    139     * @return casa::CountedPtr<Scantable> which holds a calibrated Scantable
     139    * @return casa::CountedPtr<Scantable> which holds a calibrated Scantable 
    140140    * (spectrum - average of the two CAL on and off spectra;
    141141    * tsys - mean Tsys = <caloff>*Tcal/<calon-caloff> + Tcal/2)
    142     */
     142    */             
    143143  casa::CountedPtr<Scantable> dototalpower( const casa::CountedPtr<Scantable>& calon,
    144144                                            const casa::CountedPtr<Scantable>& caloff,
     
    151151    * @param smoothref optional Boxcar smooth width of the reference scans
    152152    * default: no smoothing (=1)
    153     * @param tsysv optional scalar Tsys value at the zenith, required to
    154     * set tau, as well
     153    * @param tsysv optional scalar Tsys value at the zenith, required to 
     154    * set tau, as well 
    155155    * @param tau optional scalar Tau value
    156156    * @return casa::CountedPtr<Scantable> which holds combined scans
     
    163163                                        casa::Float tau=0.0 );
    164164
    165  /**
     165  /**
    166166    * Calibrate GBT Nod scan pairs (translated from GBTIDL)
    167167    * @param s Scantable which contains Nod scans
     
    170170    * @param tsysv optional scalar Tsys value at the zenith, required to
    171171    * set tau, as well
    172     * @param tau optional scalar Tau value
     172    * @param tau optional scalar Tau value 
    173173    * @param tcal optional scalar Tcal, CAL temperature (K)
    174174    * @return casa::CountedPtr<Scantable> which holds calibrated scans
     
    266266              double width);
    267267
     268  // test for average spectra with different channel/resolution
     269  casa::CountedPtr<Scantable>
     270    new_average( const std::vector<casa::CountedPtr<Scantable> >& in,
     271                 const bool& compel,
     272                 const std::vector<bool>& mask = std::vector<bool>(),
     273                 const std::string& weight = "NONE",
     274                 const std::string& avmode = "SCAN" )
     275    throw (casa::AipsError) ;
     276
     277
    268278private:
    269279  casa::CountedPtr<Scantable>  applyToPol( const casa::CountedPtr<Scantable>& in,
  • branches/alma/src/STMathWrapper.h

    r1388 r1446  
    188188  { return ScantableWrapper(STMath::lagFlag(in.getCP(), frequency, width)); }
    189189
     190  // test for average spectra with different channel/resolution
     191  ScantableWrapper
     192    new_average( const std::vector<ScantableWrapper>& in,
     193                 const bool& compel,
     194                 const std::vector<bool>& mask,
     195                 const std::string& weight,
     196                 const std::string& avmode )
     197  {
     198    std::vector<casa::CountedPtr<Scantable> > sts;
     199    for (unsigned int i=0; i<in.size(); ++i) sts.push_back(in[i].getCP());
     200    return ScantableWrapper(STMath::new_average(sts, compel, mask, weight, avmode));
     201  }
     202
    190203};
    191204
  • branches/alma/src/STMolecules.cpp

    r870 r1446  
    1414#include <tables/Tables/SetupNewTab.h>
    1515#include <tables/Tables/ScaColDesc.h>
     16#include <tables/Tables/ArrColDesc.h>
    1617#include <tables/Tables/TableRecord.h>
    1718#include <tables/Tables/TableParse.h>
    1819#include <tables/Tables/TableRow.h>
    1920#include <casa/Containers/RecordField.h>
     21
     22#include <tables/Tables/TableProxy.h>
    2023
    2124#include "STMolecules.h"
     
    5962{
    6063  // add to base class table
    61   table_.addColumn(ScalarColumnDesc<Double>("RESTFREQUENCY"));
    62   table_.addColumn(ScalarColumnDesc<String>("NAME"));
    63   table_.addColumn(ScalarColumnDesc<String>("FORMATTEDNAME"));
     64  //table_.addColumn(ScalarColumnDesc<Double>("RESTFREQUENCY"));
     65  table_.addColumn(ArrayColumnDesc<Double>("RESTFREQUENCY"));
     66  //table_.addColumn(ScalarColumnDesc<String>("NAME"));
     67  table_.addColumn(ArrayColumnDesc<String>("NAME"));
     68  //table_.addColumn(ScalarColumnDesc<String>("FORMATTEDNAME"));
     69  table_.addColumn(ArrayColumnDesc<String>("FORMATTEDNAME"));
    6470  table_.rwKeywordSet().define("UNIT", String("Hz"));
    6571  // new cached columns
     
    6975}
    7076
     77/***
    7178uInt STMolecules::addEntry( Double restfreq, const String& name,
    7279                            const String& formattedname )
     
    94101  return resultid;
    95102}
    96 
     103***/
     104uInt STMolecules::addEntry( Vector<Double> restfreq, const Vector<String>& name,
     105                            const Vector<String>& formattedname )
     106{
     107// How to handle this...?
     108  Table result =
     109    table_( nelements(table_.col("RESTFREQUENCY")) == restfreq.size() &&
     110            all(table_.col("RESTFREQUENCY")== restfreq) );
     111  uInt resultid = 0;
     112  if ( result.nrow() > 0) {
     113    ROScalarColumn<uInt> c(result, "ID");
     114    c.get(0, resultid);
     115  } else {
     116    uInt rno = table_.nrow();
     117    table_.addRow();
     118    // get last assigned _id and increment
     119    if ( rno > 0 ) {
     120      idCol_.get(rno-1, resultid);
     121      resultid++;
     122    }
     123    restfreqCol_.put(rno, restfreq);
     124    nameCol_.put(rno, name);
     125    formattednameCol_.put(rno, formattedname);
     126    idCol_.put(rno, resultid);
     127  }
     128  return resultid;
     129}
     130
     131
     132
     133
     134/***
    97135void STMolecules::getEntry( Double& restfreq, String& name,
    98136                            String& formattedname, uInt id ) const
     
    104142  ROTableRow row(t);
    105143  // get first row - there should only be one matching id
     144
    106145  const TableRecord& rec = row.get(0);
    107146  restfreq = rec.asDouble("RESTFREQUENCY");
     
    109148  formattedname = rec.asString("FORMATTEDNAME");
    110149}
     150***/
     151void STMolecules::getEntry( Vector<Double>& restfreq, Vector<String>& name,
     152                            Vector<String>& formattedname, uInt id ) const
     153{
     154  Table t = table_(table_.col("ID") == Int(id) );
     155  if (t.nrow() == 0 ) {
     156    throw(AipsError("STMolecules::getEntry - id out of range"));
     157  }
     158  ROTableRow row(t);
     159  // get first row - there should only be one matching id
     160
     161  const TableRecord& rec = row.get(0);
     162  //restfreq = rec.asDouble("RESTFREQUENCY");
     163  restfreq = rec.asArrayDouble("RESTFREQUENCY");
     164  //name = rec.asString("NAME");
     165  name = rec.asArrayString("NAME");
     166  //formattedname = rec.asString("FORMATTEDNAME");
     167  formattedname = rec.asArrayString("FORMATTEDNAME");
     168}
    111169
    112170std::vector< double > asap::STMolecules::getRestFrequencies( ) const
    113171{
    114172  std::vector<double> out;
     173  //TableProxy itsTable(table_);
     174  //Record rec;
    115175  Vector<Double> rfs = restfreqCol_.getColumn();
    116176  rfs.tovector(out);
     177  //rec = itsTable.getVarColumn("RESTFREQUENCY", 0, 1, 1);
    117178  return out;
    118179}
    119180
    120 double asap::STMolecules::getRestFrequency( uInt id ) const
    121 {
     181std::vector< double > asap::STMolecules::getRestFrequency( uInt id ) const
     182{
     183  std::vector<double> out;
    122184  Table t = table_(table_.col("ID") == Int(id) );
    123185  if (t.nrow() == 0 ) {
     
    126188  ROTableRow row(t);
    127189  const TableRecord& rec = row.get(0);
    128   return double(rec.asDouble("RESTFREQUENCY"));
    129 }
    130 
     190  //return double(rec.asDouble("RESTFREQUENCY"));
     191  Vector<Double> rfs = rec.asArrayDouble("RESTFREQUENCY");
     192  rfs.tovector(out);
     193  return out;
     194}
     195
     196int asap::STMolecules::nrow() const
     197{
     198  return int(table_.nrow());
     199}
    131200
    132201}//namespace
  • branches/alma/src/STMolecules.h

    r1353 r1446  
    1717#include <tables/Tables/Table.h>
    1818#include <tables/Tables/ScalarColumn.h>
     19#include <tables/Tables/ArrayColumn.h>
     20#include <casa/Arrays/Array.h>
    1921
    2022#include "STSubTable.h"
     
    3739  STMolecules& operator=(const STMolecules& other);
    3840
     41/***
    3942  casa::uInt addEntry( casa::Double restfreq, const casa::String& name="",
    4043                       const casa::String& formattedname="");
     44***/
    4145
     46  casa::uInt addEntry( casa::Vector<casa::Double> restfreq, const casa::Vector<casa::String>& name=casa::Vector<casa::String>(0),
     47                       const casa::Vector<casa::String>& formattedname=casa::Vector<casa::String>(0));
     48
     49/***
    4250  void getEntry( casa::Double& restfreq, casa::String& name,
    4351                 casa::String& formattedname, casa::uInt id) const;
     52***/
     53  void getEntry( casa::Vector<casa::Double>& restfreq, casa::Vector<casa::String>& name,
     54                 casa::Vector<casa::String>& formattedname, casa::uInt id) const;
    4455
    4556  std::vector<double> getRestFrequencies() const;
    46   double getRestFrequency( casa::uInt id ) const;
     57  std::vector<double> getRestFrequency( casa::uInt id ) const;
    4758  const casa::String& name() const { return name_; }
     59  int nrow() const;
    4860
    4961private:
     
    5264  //casa::Table table_;
    5365  //casa::ScalarColumn<casa::uInt> freqidCol_;
    54   casa::ScalarColumn<casa::Double> restfreqCol_;
    55   casa::ScalarColumn<casa::String> nameCol_;
    56   casa::ScalarColumn<casa::String> formattednameCol_; // e.g. latex
     66  //casa::ScalarColumn<casa::Double> restfreqCol_;
     67  casa::ArrayColumn<casa::Double> restfreqCol_;
     68  //casa::ScalarColumn<casa::String> nameCol_;
     69  casa::ArrayColumn<casa::String> nameCol_;
     70  //casa::ScalarColumn<casa::String> formattednameCol_; // e.g. latex
     71  casa::ArrayColumn<casa::String> formattednameCol_; // e.g. latex
    5772
    5873};
  • branches/alma/src/STWriter.cpp

    r1387 r1446  
    165165      Table btable = beamit.table();
    166166      // position only varies by beam
     167      // No, we want to pointing data which varies by cycle!
    167168      MDirection::ScalarColumn dirCol(btable, "DIRECTION");
    168169      Vector<Double> direction = dirCol(0).getAngle("rad").getValue();
     
    181182      while (!cycit.pastEnd() ) {
    182183        Table ctable = cycit.table();
     184        //MDirection::ScalarColumn dirCol(ctable, "DIRECTION");
     185        //Vector<Double> direction = dirCol(0).getAngle("rad").getValue();
    183186        TableIterator ifit(ctable, "IFNO");
    184187        Int ifno = 1;
     
    191194          ifno = rec.asuInt("IFNO")+1;
    192195          uInt nchan = specCol(0).nelements();
    193           Double cdelt,crval,crpix, restfreq;
     196          //Double cdelt,crval,crpix, restfreq;
     197          Double cdelt,crval,crpix;
     198          Vector<Double> restfreq;
    194199          Float focusAxi, focusTan, focusRot,
    195200                temperature, pressure, humidity, windSpeed, windAz;
    196201          Float tmp0,tmp1,tmp2,tmp3,tmp4;
    197202          Vector<Float> tcalval;
    198           String stmp0,stmp1, tcalt;
     203          //String stmp0,stmp1, tcalt;
     204          String tcalt;
     205          Vector<String> stmp0, stmp1;
    199206          in->frequencies().getEntry(crpix,crval,cdelt, rec.asuInt("FREQ_ID"));
    200207          in->focus().getEntry(focusAxi, focusTan, focusRot,
     
    215222          Vector<Float> tsys = tsysFromTable(itable);
    216223          // dummy data
     224          //uInt npol;
     225          //if ( hdr.antennaname == "GBT" ) {
     226          //  npol = nPolUsed;
     227          //}
     228          //else {   
     229          //  npol = specs.ncolumn();
     230          //}
    217231          uInt npol = specs.ncolumn();
    218232
     
    231245                                      refFreqNew, nchan*abs(cdelt), cdelt,
    232246                                      restfreq,
    233                                       tcal,
     247                                      tcalval,
    234248                                      tcalt,
    235249                                      rec.asFloat("AZIMUTH"),
     
    253267          }
    254268          ++count;
    255           //++ifno;
     269          ++ifno;
    256270          ++ifit;
    257271        }
     
    259273        ++cycit;
    260274      }
    261       //++beamno;
     275      ++beamno;
    262276      ++beamit;
    263277    }
     
    272286  if ( format_ == "MS2" ) {
    273287    replacePtTab(table, filename);
    274   }
     288  } 
    275289  return 0;
    276290}
  • branches/alma/src/Scantable.cpp

    r1400 r1446  
    1111//
    1212#include <map>
     13#include <fstream>
    1314
    1415#include <casa/aips.h>
     
    2425#include <casa/Arrays/Vector.h>
    2526#include <casa/Arrays/VectorSTLIterator.h>
     27#include <casa/Arrays/Slice.h>
    2628#include <casa/BasicMath/Math.h>
    2729#include <casa/BasicSL/Constants.h>
     
    533535    return int(n);
    534536  } else {
    535     // take the first SCANNO,POLNO,BEAMNO,CYCLENO as nbeam shouldn't
    536     // vary with these
     537    // take the first SCANNO,POLNO,BEAMNO,CYCLENO as nbeam shouldn't vary with these
    537538    Table t = table_(table_.col("IFNO") == ifno);
    538539    if ( t.nrow() == 0 ) return 0;
     
    786787  table_.keywordSet().get("FluxUnit", tmp);
    787788  oss << setw(15) << "Flux Unit:" << tmp << endl;
    788   Vector<Double> vec(moleculeTable_.getRestFrequencies());
     789  //Vector<Double> vec(moleculeTable_.getRestFrequencies());
     790  int nid = moleculeTable_.nrow();
     791  Bool firstline = True;
    789792  oss << setw(15) << "Rest Freqs:";
    790   if (vec.nelements() > 0) {
    791       oss << setprecision(10) << vec << " [Hz]" << endl;
    792   } else {
    793       oss << "none" << endl;
     793  for (int i=0; i<nid; i++) {
     794      Table t = table_(table_.col("MOLECULE_ID") == i);
     795      if (t.nrow() >  0) {
     796          Vector<Double> vec(moleculeTable_.getRestFrequency(i));
     797          if (vec.nelements() > 0) {
     798               if (firstline) {
     799                   oss << setprecision(10) << vec << " [Hz]" << endl;
     800                   firstline=False;
     801               }
     802               else{
     803                   oss << setw(15)<<" " << setprecision(10) << vec << " [Hz]" << endl;
     804               }
     805          } else {
     806              oss << "none" << endl;
     807          }
     808      }
    794809  }
    795810
     
    843858        oss << setw(10) << "";
    844859        oss << setw(3) << std::right << irec.asuInt("IFNO") << std::left
    845             << setw(2) << "" << frequencies().print(irec.asuInt("FREQ_ID"))
    846             << endl;
     860            << setw(2) << "" << frequencies().print(irec.asuInt("FREQ_ID"));
    847861
    848862        ++iiter;
     
    891905  const MDirection& md = getDirection(whichrow);
    892906  const MEpoch& me = timeCol_(whichrow);
    893   Double rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
     907  //Double rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
     908  Vector<Double> rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
    894909  SpectralCoordinate spc =
    895910    freqTable_.getSpectralCoordinate(md, mp, me, rf, mfreqidCol_(whichrow));
     
    950965  const MDirection& md = getDirection(whichrow);
    951966  const MEpoch& me = timeCol_(whichrow);
    952   const Double& rf = mmolidCol_(whichrow);
     967  //const Double& rf = mmolidCol_(whichrow);
     968  const Vector<Double> rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
    953969  SpectralCoordinate spc =
    954970    freqTable_.getSpectralCoordinate(md, mp, me, rf, mfreqidCol_(whichrow));
     
    967983}
    968984
    969 void Scantable::setRestFrequencies( double rf, const std::string& name,
     985/**
     986void asap::Scantable::setRestFrequencies( double rf, const std::string& name,
    970987                                          const std::string& unit )
     988**/
     989void Scantable::setRestFrequencies( vector<double> rf, const vector<std::string>& name,
     990                                          const std::string& unit )
     991
    971992{
    972993  ///@todo lookup in line table to fill in name and formattedname
    973994  Unit u(unit);
    974   Quantum<Double> urf(rf, u);
    975   uInt id = moleculeTable_.addEntry(urf.getValue("Hz"), name, "");
     995  //Quantum<Double> urf(rf, u);
     996  Quantum<Vector<Double> >urf(rf, u);
     997  Vector<String> formattedname(0);
     998  //cerr<<"Scantable::setRestFrequnecies="<<urf<<endl;
     999 
     1000  //uInt id = moleculeTable_.addEntry(urf.getValue("Hz"), name, "");
     1001  uInt id = moleculeTable_.addEntry(urf.getValue("Hz"), mathutil::toVectorString(name), formattedname);
    9761002  TableVector<uInt> tabvec(table_, "MOLECULE_ID");
    9771003  tabvec = id;
    9781004}
    9791005
    980 void Scantable::setRestFrequencies( const std::string& name )
     1006/**
     1007void asap::Scantable::setRestFrequencies( const std::string& name )
    9811008{
    9821009  throw(AipsError("setRestFrequencies( const std::string& name ) NYI"));
     1010  ///@todo implement
     1011}
     1012**/
     1013void Scantable::setRestFrequencies( const vector<std::string>& name )
     1014{
     1015  throw(AipsError("setRestFrequencies( const vector<std::string>& name ) NYI"));
    9831016  ///@todo implement
    9841017}
     
    11171150}
    11181151
    1119 }
    1120  //namespace asap
     1152void asap::Scantable::reshapeSpectrum( int nmin, int nmax )
     1153  throw( casa::AipsError )
     1154{
     1155  // assumed that all rows have same nChan
     1156  Vector<Float> arr = specCol_( 0 ) ;
     1157  int nChan = arr.nelements() ;
     1158 
     1159  // if nmin < 0 or nmax < 0, nothing to do
     1160  if (  nmin < 0 ) {
     1161    throw( casa::indexError<int>( nmin, "asap::Scantable::reshapeSpectrum: Invalid range. Negative index is specified." ) ) ;
     1162    }
     1163  if (  nmax < 0  ) {
     1164    throw( casa::indexError<int>( nmax, "asap::Scantable::reshapeSpectrum: Invalid range. Negative index is specified." ) ) ;
     1165  }
     1166 
     1167  // if nmin > nmax, exchange values
     1168  if ( nmin > nmax ) {
     1169    int tmp = nmax ;
     1170    nmax = nmin ;
     1171    nmin = tmp ;
     1172    cout << "Swap values. Applied range is ["
     1173         << nmin << ", " << nmax << "]" << endl ;
     1174  }
     1175 
     1176  // if nmin exceeds nChan, nothing to do
     1177  if ( nmin >= nChan ) {
     1178    throw( casa::indexError<int>( nmin, "asap::Scantable::reshapeSpectrum: Invalid range. Specified minimum exceeds nChan." ) ) ;
     1179  }
     1180 
     1181  // if nmax exceeds nChan, reset nmax to nChan
     1182  if ( nmax >= nChan ) {
     1183    if ( nmin == 0 ) {
     1184      // nothing to do
     1185      cout << "Whole range is selected. Nothing to do." << endl ;
     1186      return ;
     1187    }
     1188    else {
     1189      cout << "Specified maximum exceeds nChan. Applied range is ["
     1190           << nmin << ", " << nChan-1 << "]." << endl ;
     1191      nmax = nChan - 1 ;
     1192    }
     1193  }
     1194 
     1195  // reshape specCol_ and flagCol_
     1196  for ( int irow = 0 ; irow < nrow() ; irow++ ) {
     1197    reshapeSpectrum( nmin, nmax, irow ) ;
     1198  }
     1199
     1200  // update FREQUENCIES subtable
     1201  Double refpix ;
     1202  Double refval ;
     1203  Double increment ;
     1204  int freqnrow = freqTable_.table().nrow() ;
     1205  Vector<uInt> oldId( freqnrow ) ;
     1206  Vector<uInt> newId( freqnrow ) ;
     1207  for ( int irow = 0 ; irow < freqnrow ; irow++ ) {
     1208    freqTable_.getEntry( refpix, refval, increment, irow ) ;
     1209    /***
     1210     * need to shift refpix to nmin
     1211     * note that channel nmin in old index will be channel 0 in new one
     1212     ***/
     1213    refval = refval - ( refpix - nmin ) * increment ;
     1214    refpix = 0 ; 
     1215    freqTable_.setEntry( refpix, refval, increment, irow ) ;
     1216  }
     1217 
     1218  // update nchan
     1219  int newsize = nmax - nmin + 1 ;
     1220  table_.rwKeywordSet().define( "nChan", newsize ) ;
     1221 
     1222  // update bandwidth
     1223  // assumed all spectra in the scantable have same bandwidth
     1224  table_.rwKeywordSet().define( "Bandwidth", increment * newsize ) ;
     1225 
     1226  return ;
     1227}
     1228
     1229void asap::Scantable::reshapeSpectrum( int nmin, int nmax, int irow )
     1230{
     1231  // reshape specCol_ and flagCol_
     1232  Vector<Float> oldspec = specCol_( irow ) ;
     1233  Vector<uChar> oldflag = flagsCol_( irow ) ;
     1234  uInt newsize = nmax - nmin + 1 ;
     1235  specCol_.put( irow, oldspec( Slice( nmin, newsize, 1 ) ) ) ;
     1236  flagsCol_.put( irow, oldflag( Slice( nmin, newsize, 1 ) ) ) ;
     1237
     1238  return ;
     1239}
     1240
     1241void asap::Scantable::regridChannel( int nChan, double dnu )
     1242{
     1243  cout << "Regrid abcissa with channel number " << nChan << " and spectral resoultion " << dnu << "Hz." << endl ;
     1244  // assumed that all rows have same nChan
     1245  Vector<Float> arr = specCol_( 0 ) ;
     1246  int oldsize = arr.nelements() ;
     1247
     1248  // if oldsize == nChan, nothing to do
     1249  if ( oldsize == nChan ) {
     1250    cout << "Specified channel number is same as current one. Nothing to do." << endl ;
     1251    return ;
     1252  }
     1253
     1254  // if oldChan < nChan, unphysical operation
     1255  if ( oldsize < nChan ) {
     1256    cout << "Unphysical operation. Nothing to do." << endl ;
     1257    return ;
     1258  }
     1259
     1260  // change channel number for specCol_ and flagCol_
     1261  Vector<Float> newspec( nChan, 0 ) ;
     1262  Vector<uChar> newflag( nChan, false ) ;
     1263  vector<string> coordinfo = getCoordInfo() ;
     1264  string oldinfo = coordinfo[0] ;
     1265  coordinfo[0] = "Hz" ;
     1266  setCoordInfo( coordinfo ) ;
     1267  for ( int irow = 0 ; irow < nrow() ; irow++ ) {
     1268    regridChannel( nChan, dnu, irow ) ;
     1269  }
     1270  coordinfo[0] = oldinfo ;
     1271  setCoordInfo( coordinfo ) ;
     1272
     1273
     1274  // NOTE: this method does not update metadata such as
     1275  //       FREQUENCIES subtable, nChan, Bandwidth, etc.
     1276
     1277  return ;
     1278}
     1279
     1280void asap::Scantable::regridChannel( int nChan, double dnu, int irow )
     1281{
     1282  // logging
     1283  //ofstream ofs( "average.log", std::ios::out | std::ios::app ) ;
     1284  //ofs << "IFNO = " << getIF( irow ) << " irow = " << irow << endl ;
     1285
     1286  Vector<Float> oldspec = specCol_( irow ) ;
     1287  Vector<uChar> oldflag = flagsCol_( irow ) ;
     1288  Vector<Float> newspec( nChan, 0 ) ;
     1289  Vector<uChar> newflag( nChan, false ) ;
     1290 
     1291  // regrid
     1292  vector<double> abcissa = getAbcissa( irow ) ;
     1293  int oldsize = abcissa.size() ;
     1294  double olddnu = abcissa[1] - abcissa[0] ;
     1295  int refChan = 0 ;
     1296  double frac = 0.0 ;
     1297  double wedge = 0.0 ;
     1298  double pile = 0.0 ;
     1299  double wsum = 0.0 ;
     1300  /***
     1301   * ichan = 0
     1302   ***/
     1303  //ofs << "olddnu = " << olddnu << ", dnu = " << dnu << endl ;
     1304  pile += dnu ;
     1305  wedge = olddnu * ( refChan + 1 ) ;
     1306  while ( wedge < pile ) {
     1307    newspec[0] += olddnu * oldspec[refChan] ;
     1308    newflag[0] = newflag[0] || oldflag[refChan] ;
     1309    //ofs << "channel " << refChan << " is included in new channel 0" << endl ;
     1310    refChan++ ;
     1311    wedge += olddnu ;
     1312    wsum += olddnu ;
     1313    //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ;
     1314  }
     1315  frac = ( wedge - pile ) / olddnu ;
     1316  wsum += ( 1.0 - frac ) * olddnu ;
     1317  newspec[0] += ( 1.0 - frac ) * olddnu * oldspec[refChan] ;
     1318  newflag[0] = newflag[0] || oldflag[refChan] ;
     1319  //ofs << "channel " << refChan << " is partly included in new channel 0" << " with fraction of " << ( 1.0 - frac ) << endl ;
     1320  //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ;
     1321  newspec[0] /= wsum ;
     1322  //ofs << "newspec[0] = " << newspec[0] << endl ;
     1323  //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
     1324
     1325  /***
     1326   * ichan = 1 - nChan-2
     1327   ***/
     1328  for ( int ichan = 1 ; ichan < nChan - 1 ; ichan++ ) {
     1329    pile += dnu ;
     1330    newspec[ichan] += frac * olddnu * oldspec[refChan] ;
     1331    newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
     1332    //ofs << "channel " << refChan << " is partly included in new channel " << ichan << " with fraction of " << frac << endl ;
     1333    refChan++ ;
     1334    wedge += olddnu ;
     1335    wsum = frac * olddnu ;
     1336    //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
     1337    while ( wedge < pile ) {
     1338      newspec[ichan] += olddnu * oldspec[refChan] ;
     1339      newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
     1340      //ofs << "channel " << refChan << " is included in new channel " << ichan << endl ;
     1341      refChan++ ;
     1342      wedge += olddnu ;
     1343      wsum += olddnu ;
     1344      //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
     1345    }
     1346    frac = ( wedge - pile ) / olddnu ;
     1347    wsum += ( 1.0 - frac ) * olddnu ;
     1348    newspec[ichan] += ( 1.0 - frac ) * olddnu * oldspec[refChan] ;
     1349    newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
     1350    //ofs << "channel " << refChan << " is partly included in new channel " << ichan << " with fraction of " << ( 1.0 - frac ) << endl ;
     1351    //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
     1352    //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
     1353    newspec[ichan] /= wsum ;
     1354    //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << endl ;
     1355  }
     1356
     1357  /***
     1358   * ichan = nChan-1
     1359   ***/
     1360  // NOTE: Assumed that all spectra have the same bandwidth
     1361  pile += dnu ;
     1362  newspec[nChan-1] += frac * olddnu * oldspec[refChan] ;
     1363  newflag[nChan-1] = newflag[nChan-1] || oldflag[refChan] ;
     1364  //ofs << "channel " << refChan << " is partly included in new channel " << nChan-1 << " with fraction of " << frac << endl ;
     1365  refChan++ ;
     1366  wedge += olddnu ;
     1367  wsum = frac * olddnu ;
     1368  //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
     1369  for ( int jchan = refChan ; jchan < oldsize ; jchan++ ) {
     1370    newspec[nChan-1] += olddnu * oldspec[jchan] ;
     1371    newflag[nChan-1] = newflag[nChan-1] || oldflag[jchan] ;
     1372    wsum += olddnu ;
     1373    //ofs << "channel " << jchan << " is included in new channel " << nChan-1 << " with fraction of " << frac << endl ;
     1374    //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
     1375  }
     1376  //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
     1377  //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
     1378  newspec[nChan-1] /= wsum ;
     1379  //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << endl ;
     1380 
     1381  specCol_.put( irow, newspec ) ;
     1382  flagsCol_.put( irow, newflag ) ;
     1383
     1384  // ofs.close() ;
     1385
     1386  return ;
     1387}
     1388
     1389}
     1390//namespace asap
  • branches/alma/src/Scantable.h

    r1400 r1446  
    2222#include <casa/Utilities/CountedPtr.h>
    2323
     24#include <casa/Exceptions/Error.h>
     25
    2426#include <coordinates/Coordinates/SpectralCoordinate.h>
    2527
     
    2931
    3032#include <measures/TableMeasures/ScalarMeasColumn.h>
     33
     34#include <casa/Arrays/Vector.h>
     35#include <casa/Quanta/Quantum.h>
    3136
    3237#include "Logger.h"
     
    167172
    168173        /**
    169          * get the direction type as a string, e.g. "J2000"
     174         * get the direction as a string
    170175         * @param[in] whichrow the row number
    171176         * return the direction string
     
    277282  std::vector<uint> getScanNos() { return getNumbers(scanCol_); }
    278283  int getScan(int whichrow) const { return scanCol_(whichrow); }
     284
     285  //TT addition
     286  std::vector<uint> getMolNos() {return getNumbers(mmolidCol_); }
    279287
    280288  /**
     
    352360  std::vector<double> getRestFrequencies() const
    353361    { return moleculeTable_.getRestFrequencies(); }
    354 
     362  std::vector<double> getRestFrequency(int id) const
     363    { return moleculeTable_.getRestFrequency(id); }
     364
     365  /**
    355366  void setRestFrequencies(double rf, const std::string& name = "",
    356367                          const std::string& = "Hz");
    357   void setRestFrequencies(const std::string& name);
     368  **/
     369  // Modified by Takeshi Nakazato 05/09/2008
     370  /***
     371  void setRestFrequencies(vector<double> rf, const vector<std::string>& name = "",
     372                          const std::string& = "Hz");
     373  ***/
     374  void setRestFrequencies(vector<double> rf,
     375                          const vector<std::string>& name = vector<std::string>(1,""),
     376                          const std::string& = "Hz");
     377
     378  //void setRestFrequencies(const std::string& name);
     379  void setRestFrequencies(const vector<std::string>& name);
    358380
    359381  void shift(int npix);
     
    393415   * against the information found in GBT_GO table for
    394416   * scan number orders to get correct pairs.
    395    * @param[in] scan list
    396    * @return status
     417   * @param[in] scan list 
     418   * @return status 
    397419   */
    398420  int checkScanInfo(const std::vector<int>& scanlist) const;
    399421
    400422  /**
    401    * Get the direction as a vector, for a specific row
     423   * Get the direction as a vector, for a specific row 
    402424   * @param[in] whichrow the row numbyyer
    403    * @return the direction in a vector
     425   * @return the direction in a vector 
    404426   */
    405427  std::vector<double> getDirectionVector(int whichrow) const;
    406428
     429  /**
     430   * Reshape spectrum
     431   * @param[in] nmin, nmax minimum and maximum channel
     432   * @param[in] irow       row number
     433   *
     434   * 30/07/2008 Takeshi Nakazato 
     435   **/
     436  void reshapeSpectrum( int nmin, int nmax ) throw( casa::AipsError );
     437  void reshapeSpectrum( int nmin, int nmax, int irow ) ;
     438
     439  /**
     440   * Change channel number under fixed bandwidth
     441   * @param[in] nchan, dnu new channel number and spectral resolution
     442   * @param[in] irow       row number
     443   *
     444   * 27/08/2008 Takeshi Nakazato
     445   **/
     446  void regridChannel( int nchan, double dnu ) ;
     447  void regridChannel( int nchan, double dnu, int irow ) ;
     448 
    407449private:
    408450
  • branches/alma/src/ScantableWrapper.h

    r1400 r1446  
    133133  std::vector<uint> getScanNos() { return table_->getScanNos(); }
    134134  int getScan(int whichrow) const {return table_->getScan(whichrow);}
     135  std::vector<uint> getMolNos() { return table_->getMolNos();}
    135136
    136137  STSelector getSelection() const { return table_->getSelection(); }
     
    157158  { table_->shift(npix); }
    158159
     160/**
     161  commented out by TT
    159162  void setRestFrequencies(double rf, const std::string& name,
    160163                          const std::string& unit)
    161164    { table_->setRestFrequencies(rf, name, unit); }
     165**/
     166  void setRestFrequencies(vector<double> rf, const vector<std::string>& name,
     167                          const std::string& unit)
     168    { table_->setRestFrequencies(rf, name, unit); }
     169
    162170/*
    163171  void setRestFrequencies(const std::string& name) {
     
    166174*/
    167175
     176/*
    168177  std::vector<double> getRestFrequencies() const
    169178    { return table_->getRestFrequencies(); }
     179*/
     180  std::vector<double> getRestFrequency(int id) const
     181    { return table_->getRestFrequency(id); }
    170182
    171183  void setCoordInfo(std::vector<string> theinfo) {
     
    208220  int checkScanInfo(const vector<int>& scanlist) const
    209221    { return table_->checkScanInfo(scanlist); }
    210 
     222 
    211223  std::vector<double>  getDirectionVector(int whichrow) const
    212224    { return table_->getDirectionVector(whichrow); }
     225
     226  void reshapeSpectrum( int nmin, int nmax )
     227  { table_->reshapeSpectrum( nmin, nmax ); }
    213228
    214229private:
  • branches/alma/src/Templates.cpp

    r1387 r1446  
    4141template class casa::PtrRep<asap::Scantable>;
    4242template class casa::SimpleCountedConstPtr<asap::Scantable>;
     43template class casa::SimpleCountedConstPtr<asap::STPol>;
    4344template class casa::SimpleCountedPtr<asap::Scantable>;
    4445template class casa::SimpleCountedPtr<asap::STPol>;
  • branches/alma/src/python_STMath.cpp

    r1387 r1446  
    7373        .def("_mx_extract", &STMathWrapper::mxExtract)
    7474        .def("_lag_flag", &STMathWrapper::lagFlag)
     75        // testing average spectra with different channel/resolution
     76        .def("_new_average", &STMathWrapper::new_average)
    7577          ;
    7678    };
  • branches/alma/src/python_Scantable.cpp

    r1387 r1446  
    5858    .def("getscannos", &ScantableWrapper::getScanNos)
    5959    .def("getcycle", &ScantableWrapper::getCycle)
     60    .def("getmolnos", &ScantableWrapper::getMolNos)
    6061    .def("nif", &ScantableWrapper::nif,
    6162         (boost::python::arg("scanno")=-1) )
     
    107108    .def("_summary",  &ScantableWrapper::summary,
    108109         (boost::python::arg("verbose")=true) )
    109     .def("_getrestfreqs",  &ScantableWrapper::getRestFrequencies)
     110    //.def("_getrestfreqs",  &ScantableWrapper::getRestFrequencies)
     111    .def("_getrestfreqs",  &ScantableWrapper::getRestFrequency)
    110112    .def("_setrestfreqs",  &ScantableWrapper::setRestFrequencies)
    111113    .def("shift_refpix", &ScantableWrapper::shift)
     
    123125    .def("_setsourcetype", &ScantableWrapper::setSourceType)
    124126    .def("_getdirectionvec", &ScantableWrapper::getDirectionVector)
     127    .def("_reshape", &ScantableWrapper::reshapeSpectrum,
     128         (boost::python::arg("nmin")=-1,
     129          boost::python::arg("nmax")=-1) )
    125130  ;
    126131};
Note: See TracChangeset for help on using the changeset viewer.