-                      ra36c6d3
+                      r8115d82
 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
 D 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],
+]

Note: See TracChangeset for help on using the changeset viewer.

SasView

Changeset 8115d82 in sasmodels

Legend:

sasmodels/models/two_power_law.py

Download in other formats: