# Changeset 8115d82 in sasmodels

Ignore:
Timestamp:
Feb 10, 2016 5:53:39 AM (8 years ago)
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
34d6cab
Parents:
a36c6d3
Message:

Two Power Law code converted from old model and ready for testing

Location:
sasmodels/models
Files:
 ra36c6d3 r""" This model calculates an empirical functional form for SAS data characterized by two Lorentzian-type functions. by two power laws. Definition .. math:: I(q) = \frac{A}{1 +(Q\xi_1)^n} + \frac{C}{1 +(Q\xi_2)^m} + \text{B} I(q) = \begin{cases} A \pow(qval,-m1) & q <= qc \\ C \pow(qval,-m2) & q > qc \end{cases} where $A$ = Lorentzian scale factor #1, $C$ = Lorentzian scale #2, $\xi_1$ and $\xi_2$ are the corresponding correlation lengths, and $n$ and $m$ are the respective power law exponents (set $n = m = 2$ for Ornstein-Zernicke behaviour). where $qc$ = the location of the crossover from one slope to the other, $A$ = the scaling coefficent that sets the overall intensity of the lower Q power law region, $m1$ = power law exponent at low Q $m2$ = power law exponent at high Q The scaling of the second power law region (coefficent C) is then automatically scaled to match the first by following formula .. math:: C = frac{A}{\pow(qc/-m1)\pow(qc/-m2)} ...note Be sure to enter the power law exponents as positive values! For 2D data the scattering intensity is calculated in the same way as 1D, .. figure:: img/two_lorentzian.jpg .. figure:: img/two_power_law_1d.jpg 1D plot using the default values (w/500 data point). """ from numpy import power from numpy import sqrt from sas.models.BaseComponent import BaseComponent from numpy import power import math name = "two_power_law" title = "Two Power Law intensity calculation" description = """ I(q) = coef_A*pow(qval,-1.0*power1) + background for q<=qc =C*pow(qval,-1.0*power2) + background for q>qc where C=coef_A*pow(qc,-1.0*power1)/pow(qc,-1.0*power2). class TwoPowerLawModel(BaseComponent): """ Class that evaluates a TwoPowerLawModel. coef_A = scaling coefficent qc = crossover location [1/A] power_1 (=m1) = power law exponent at low Q power_2 (=m2) = power law exponent at high Q background = Incoherent background [1/cm] """ category = "shape-independent" Calculate:: # pylint: disable=bad-whitespace, line-too-long #            ["name", "units", default, [lower, upper], "type", "description"], parameters = [["coef_A",  "",     1.0, [-inf, inf], "", "scaling coefficent in low Q region"], ["qc", "1/Ang", 0.04, [0, inf], "", "crossover location"], ["power_1",    "",    1.0, [0, inf], "", "power law exponent at low Q"], ["power_2",  "",      4.0, [0, inf], "", "power law exponent at high Q"], ] # pylint: enable=bad-whitespace, line-too-long I(q) = coef_A pow(qval,-power1) for q<=qc I(q) = C pow(qval,-power2) for q>qc def Iq(q, coef_A=1.0, qc=0.04, power_1=1.0, power_2=4.0, ): where C=coef_A pow(qc,-power1)/pow(qc,-power2). List of default parameters: * coef_A = coefficient * power1 = (-) Power @ low Q * power2 = (-) Power @ high Q * qc = crossover Q-value * background = incoherent background """ def __init__(self): """ Initialization """ # Initialize BaseComponent first, then sphere BaseComponent.__init__(self) ## Name of the model self.name = "TwoPowerLaw" self.description="""I(q) = coef_A*pow(qval,-1.0*power1) for q<=qc =C*pow(qval,-1.0*power2) for q>qc where C=coef_A*pow(qc,-1.0*power1)/pow(qc,-1.0*power2). List of default parameters: coef_A = coefficient power1 = (-) Power @ low Q power2 = (-) Power @ high Q qc = crossover Q-value background = incoherent background """ ## Define parameters self.params = {} self.params['coef_A']  = 1.0 self.params['power1']     = 1.0 self.params['power2']  = 4.0 self.params['qc']     = 0.04 self.params['background']     = 0.0 ## Parameter details [units, min, max] self.details = {} self.details['coef_A'] = ['', None, None] self.details['power1'] =  ['', None, None] self.details['power2']  =  ['', None, None] self.details['qc']  =   ['1/A', None, None] self.details['background']   =  ['[1/cm]', None, None] :param q:                   Input q-value (float or [float, float]) :param coef_A:              Scaling coefficent at low Q :param qc:                  Crossover location :param power_1:             Exponent of power law function at low Q :param power_2:             Exponent of power law function at high Q :return:                    Calculated intensity """ #list of parameter that cannot be fitted self.fixed= [] def _twopowerlaw(self, x): """ Model definition """ qc= self.params['qc'] if(x<=qc): inten = self.params['coef_A']*power(x,-1.0*self.params['power1']) else: scale = self.params['coef_A']*power(qc,-1.0*self.params['power1']) \ / power(qc,-1.0*self.params['power2']) inten = scale*power(x,-1.0*self.params['power2']) inten += self.params['background'] # pylint: disable=bad-whitespace if(g<=qc): intensity = coef_A*power(q,-1.0*power_1) else: coef_C = coef_A*power(qc,-1.0*power_1)/power(qc,-1.0*power_2) intensity = coef_C*power(q,-1.0*power_2) return inten def run(self, x = 0.0): """ Evaluate the model @param x: input q-value (float or [float, float] as [r, theta]) @return: (guinier value) """ if x.__class__.__name__ == 'list': return self._twopowerlaw(x[0]) elif x.__class__.__name__ == 'tuple': raise ValueError, "Tuples are not allowed as input to BaseComponent models" else: return self._twopowerlaw(x) def runXY(self, x = 0.0): """ Evaluate the model @param x: input q-value (float or [float, float] as [qx, qy]) @return: guinier value """ if x.__class__.__name__ == 'list': q = math.sqrt(x[0]**2 + x[1]**2) return self._twopowerlaw(q) elif x.__class__.__name__ == 'tuple': raise ValueError, "Tuples are not allowed as input to BaseComponent models" else: return self._twopowerlaw(x) return intensity Iq.vectorized = True  # Iq accepts an array of q values def Iqxy(qx, qy, *args): """ :param qx:   Input q_x-value :param qy:   Input q_y-value :param args: Remaining arguments :return:     2D-Intensity """ return Iq(sqrt(qx**2 + qy**2), *args) Iqxy.vectorized = True  # Iqxy accepts an array of qx, qy values demo = dict(scale=1, background=0.1, coef_A=1, qc=0.04, power_1=1.0, power_2=4.0) oldname = "TwoPowerLawModel" oldpars = dict(background='background', coef_A='coef_A', qc='qc', power_1='power1', power_2='power2') tests = [ # Accuracy tests based on content in test/utest_extra_models.py [{'coeff_A':    1.0, 'qc':         0.04, 'power_1':    1.0, 'power_2':    4.0, 'background': 0.1, }, 0.001, 1000], [{'coeff_A':    1.0, 'qc':         0.04, 'power_1':    1.0, 'power_2':    4.0, 'background': 0.1, }, 0.150141, 0.125945], [{'coeff_A':    1.0, 'qc':         0.04, 'power_1':    1.0, 'power_2':    4.0, 'background': 0.1, }, [0.442528,0.0], 0.00166884], ]