source: sasmodels/sasmodels/models/two_power_law.py

Last change on this file was c1e44e5, checked in by Paul Kienzle <pkienzle@…>, 5 years ago

Add local link to source files. Refs #1263.

  • Property mode set to 100644
File size: 4.4 KB
Line 
1r"""
2Definition
3----------
4
5The scattering intensity $I(q)$ is calculated as
6
7.. math::
8
9    I(q) = \begin{cases}
10    A q^{-m1} + \text{background} & q <= q_c \\
11    C q^{-m2} + \text{background} & q > q_c
12    \end{cases}
13
14where $q_c$ = the location of the crossover from one slope to the other,
15$A$ = the scaling coefficent that sets the overall intensity of the lower Q
16power law region, $m1$ = power law exponent at low Q, and $m2$ = power law
17exponent at high Q.  The scaling of the second power law region (coefficent C)
18is then automatically scaled to match the first by following formula:
19
20.. math::
21    C = \frac{A q_c^{m2}}{q_c^{m1}}
22
23.. note::
24    Be sure to enter the power law exponents as positive values!
25
26For 2D data the scattering intensity is calculated in the same way as 1D,
27where the $q$ vector is defined as
28
29.. math::
30
31    q = \sqrt{q_x^2 + q_y^2}
32
33
34References
35----------
36
37None.
38
39Authorship and Verification
40----------------------------
41
42* **Author:** NIST IGOR/DANSE **Date:** pre 2010
43* **Last Modified by:** Wojciech Wpotrzebowski **Date:** February 18, 2016
44* **Last Reviewed by:** Paul Butler **Date:** March 21, 2016
45"""
46
47import numpy as np
48from numpy import inf, power, empty, errstate
49
50name = "two_power_law"
51title = "This model calculates an empirical functional form for SAS data \
52characterized by two power laws."
53description = """
54            I(q) = coef_A*pow(qval,-1.0*power1) + background for q<=q_c
55            =C*pow(qval,-1.0*power2) + background for q>q_c
56            where C=coef_A*pow(q_c,-1.0*power1)/pow(q_c,-1.0*power2).
57
58            coef_A = scaling coefficent
59            q_c = crossover location [1/A]
60            power_1 (=m1) = power law exponent at low Q
61            power_2 (=m2) = power law exponent at high Q
62            background = Incoherent background [1/cm]
63        """
64category = "shape-independent"
65
66# pylint: disable=bad-whitespace, line-too-long
67#   ["name", "units", default, [lower, upper], "type", "description"],
68parameters = [
69    ["coefficent_1", "",       1.0, [-inf, inf], "", "coefficent A in low Q region"],
70    ["crossover",    "1/Ang",  0.04,[0, inf],    "", "crossover location"],
71    ["power_1",      "",       1.0, [0, inf],    "", "power law exponent at low Q"],
72    ["power_2",      "",       4.0, [0, inf],    "", "power law exponent at high Q"],
73    ]
74# pylint: enable=bad-whitespace, line-too-long
75
76
77def Iq(q,
78       coefficent_1=1.0,
79       crossover=0.04,
80       power_1=1.0,
81       power_2=4.0,
82      ):
83    """
84    :param q:                   Input q-value (float or [float, float])
85    :param coefficent_1:        Scaling coefficent at low Q
86    :param crossover:           Crossover location
87    :param power_1:             Exponent of power law function at low Q
88    :param power_2:             Exponent of power law function at high Q
89    :return:                    Calculated intensity
90    """
91    result = empty(q.shape, 'd')
92    index = (q <= crossover)
93    with errstate(divide='ignore'):
94        coefficent_2 = coefficent_1 * power(crossover, power_2 - power_1)
95        result[index] = coefficent_1 * power(q[index], -power_1)
96        result[~index] = coefficent_2 * power(q[~index], -power_2)
97    return result
98
99Iq.vectorized = True  # Iq accepts an array of q values
100
101def random():
102    """Return a random parameter set for the model."""
103    coefficient_1 = 1
104    crossover = 10**np.random.uniform(-3, -1)
105    power_1 = np.random.uniform(1, 6)
106    power_2 = np.random.uniform(1, 6)
107    pars = dict(
108        scale=1, #background=0,
109        coefficient_1=coefficient_1,
110        crossover=crossover,
111        power_1=power_1,
112        power_2=power_2,
113    )
114    return pars
115
116demo = dict(scale=1, background=0.0,
117            coefficent_1=1.0,
118            crossover=0.04,
119            power_1=1.0,
120            power_2=4.0)
121
122tests = [
123    # Accuracy tests based on content in test/utest_extra_models.py
124    [{'coefficent_1':     1.0,
125      'crossover':  0.04,
126      'power_1':    1.0,
127      'power_2':    4.0,
128      'background': 0.0,
129     }, 0.001, 1000],
130
131    [{'coefficent_1':     1.0,
132      'crossover':  0.04,
133      'power_1':    1.0,
134      'power_2':    4.0,
135      'background': 0.0,
136     }, 0.150141, 0.125945],
137
138    [{'coefficent_1':    1.0,
139      'crossover':  0.04,
140      'power_1':    1.0,
141      'power_2':    4.0,
142      'background': 0.0,
143     }, 0.442528, 0.00166884],
144
145    [{'coefficent_1':    1.0,
146      'crossover':  0.04,
147      'power_1':    1.0,
148      'power_2':    4.0,
149      'background': 0.0,
150     }, (0.442528, 0.00166884), 0.00166884],
151
152]
Note: See TracBrowser for help on using the repository browser.