source: sasmodels/sasmodels/models/two_power_law.py @ 0507e09

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 0507e09 was 0507e09, checked in by smk78, 5 years ago

Added link to source code to each model. Closes #883

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