Changeset 1701 for branches/alma


Ignore:
Timestamp:
02/16/10 16:29:35 (15 years ago)
Author:
WataruKawasaki
Message:

New Development: Yes

JIRA Issue: Yes (CAS-1800 + CAS-1807)

Ready to Release: Yes

Interface Changes: Yes

What Interface Changed: added new methods to scantable and fitter.

Test Programs:

Put in Release Notes: No

Module(s): sdfit, sdflag

Description: added new methods 'scantable.clip' and 'fitter.set_lorentz_parameters'.


Location:
branches/alma/python
Files:
2 edited

Legend:

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

    r1676 r1701  
    8282            lpoly:   use polynomial of the order given with linear least squares fit
    8383            gauss:   fit the number of gaussian specified
     84            lorentz: fit the number of lorentzian specified
    8485        Example:
    85             fitter.set_function(gauss=2) # will fit two gaussians
    8686            fitter.set_function(poly=3)  # will fit a 3rd order polynomial via nonlinear method
    8787            fitter.set_function(lpoly=3)  # will fit a 3rd order polynomial via linear method
     88            fitter.set_function(gauss=2) # will fit two gaussians
     89            fitter.set_function(lorentz=2) # will fit two lorentzians
    8890        """
    8991        #default poly order 0
     
    105107            self.components = [ 3 for i in range(n) ]
    106108            self.uselinear = False
     109        elif kwargs.has_key('lorentz'):
     110            n = kwargs.get('lorentz')
     111            self.fitfunc = 'lorentz'
     112            self.fitfuncs = [ 'lorentz' for i in range(n) ]
     113            self.components = [ 3 for i in range(n) ]
     114            self.uselinear = False
    107115        else:
    108116            msg = "Invalid function type."
     
    160168                asaplog.push(out,False)
    161169        self.fitter.setdata(self.x, self.y, self.mask)
    162         if self.fitfunc == 'gauss':
     170        if self.fitfunc == 'gauss' or self.fitfunc == 'lorentz':
    163171            ps = self.fitter.getparameters()
    164172            if len(ps) == 0 or estimate:
     
    243251            else:
    244252                raise RuntimeError(msg)
    245         if self.fitfunc == "gauss" and component is not None:
     253        if (self.fitfunc == "gauss" or self.fitfunc == 'lorentz') and component is not None:
    246254            if not self.fitted and sum(self.fitter.getparameters()) == 0:
    247255                pars = _n_bools(len(self.components)*3, False)
     
    301309                raise ValueError(msg)
    302310
     311    def set_lorentz_parameters(self, peak, centre, fwhm,
     312                             peakfixed=0, centrefixed=0,
     313                             fwhmfixed=0,
     314                             component=0):
     315        """
     316        Set the Parameters of a 'Lorentzian' component, set with set_function.
     317        Parameters:
     318            peak, centre, fwhm:  The gaussian parameters
     319            peakfixed,
     320            centrefixed,
     321            fwhmfixed:           Optional parameters to indicate if
     322                                 the paramters should be held fixed during
     323                                 the fitting process. The default is to keep
     324                                 all parameters flexible.
     325            component:           The number of the component (Default is the
     326                                 component 0)
     327        """
     328        if self.fitfunc != "lorentz":
     329            msg = "Function only operates on Lorentzian components."
     330            if rcParams['verbose']:
     331                #print msg
     332                asaplog.push(msg)
     333                print_log('ERROR')
     334                return
     335            else:
     336                raise ValueError(msg)
     337        if 0 <= component < len(self.components):
     338            d = {'params':[peak, centre, fwhm],
     339                 'fixed':[peakfixed, centrefixed, fwhmfixed]}
     340            self.set_parameters(d, component)
     341        else:
     342            msg = "Please select a valid  component."
     343            if rcParams['verbose']:
     344                #print msg
     345                asaplog.push(msg)
     346                print_log('ERROR')
     347                return
     348            else:
     349                raise ValueError(msg)
     350
    303351    def get_area(self, component=None):
    304352        """
    305         Return the area under the fitted gaussian component.
    306         Parameters:
    307               component:   the gaussian component selection,
     353        Return the area under the fitted gaussian/lorentzian component.
     354        Parameters:
     355              component:   the gaussian/lorentzian component selection,
    308356                           default (None) is the sum of all components
    309357        Note:
    310               This will only work for gaussian fits.
     358              This will only work for gaussian/lorentzian fits.
    311359        """
    312360        if not self.fitted: return
    313         if self.fitfunc == "gauss":
     361        if self.fitfunc == "gauss" or self.fitfunc == "lorentz":
    314362            pars = list(self.fitter.getparameters())
    315363            from math import log,pi,sqrt
    316             fac = sqrt(pi/log(16.0))
     364            if self.fitfunc == "gauss":
     365                fac = sqrt(pi/log(16.0))
     366            elif self.fitfunc == "lorentz":
     367                fac = pi/2.0
    317368            areas = []
    318369            for i in range(len(self.components)):
     
    346397        cerrs = errs
    347398        if component is not None:
    348             if self.fitfunc == "gauss":
     399            if self.fitfunc == "gauss" or self.fitfunc == "lorentz":
    349400                i = 3*component
    350401                if i < len(errs):
     
    373424        area = []
    374425        if component is not None:
    375             if self.fitfunc == "gauss":
     426            if self.fitfunc == "gauss" or self.fitfunc == "lorentz":
    376427                i = 3*component
    377428                cpars = pars[i:i+3]
     
    388439            cfixed = fixed
    389440            cerrs = errs
    390             if self.fitfunc == "gauss":
     441            if self.fitfunc == "gauss" or self.fitfunc == "lorentz":
    391442                for c in range(len(self.components)):
    392443                  a = self.get_area(c)
     
    413464                c+=1
    414465            out = out[:-1]  # remove trailing ','
    415         elif self.fitfunc == 'gauss':
     466        elif self.fitfunc == 'gauss' or self.fitfunc == 'lorentz':
    416467            i = 0
    417468            c = 0
  • branches/alma/python/scantable.py

    r1685 r1701  
    862862        self._add_history("flag_row", varlist)
    863863       
     864    def clip(self, uthres=None, dthres=None, clipoutside=True, unflag=False):
     865        """
     866        Flag the selected data outside a specified range (in channel-base)
     867        Parameters:
     868            uthres:      upper threshold.
     869            dthres:      lower threshold
     870            clipoutside: True for flagging data outside the range [dthres:uthres].
     871                         False for glagging data inside the range.
     872            unflag     : if True, unflag the data.
     873        """
     874        varlist = vars()
     875        try:
     876            self._clip(uthres, dthres, clipoutside, unflag)
     877        except RuntimeError, msg:
     878            if rcParams['verbose']:
     879                print_log()
     880                asaplog.push(str(msg))
     881                print_log('ERROR')
     882                return
     883            else: raise
     884        self._add_history("clip", varlist)
    864885       
    865886    def lag_flag(self, frequency, width=0.0, unit="GHz", insitu=None):
Note: See TracChangeset for help on using the changeset viewer.