source: sasmodels/sasmodels/models/unified_power_Rg.py @ ef476e6

ticket-1257-vesicle-productticket_1156ticket_822_more_unit_tests
Last change on this file since ef476e6 was ef476e6, checked in by smk78, 5 years ago

Revert some of the changes done for #1059 committed yesterday

  • Property mode set to 100644
File size: 4.8 KB
Line 
1r"""
2Definition
3----------
4
5This model employs the empirical multiple level unified Exponential/Power-law
6fit method developed by Beaucage. Four functions are included so that 1, 2, 3,
7or 4 levels can be used. In addition a 0 level has been added which simply
8calculates
9
10.. math::
11
12    I(q) = \text{scale} / q + \text{background}
13
14The Beaucage method is able to reasonably approximate the scattering from
15many different types of particles, including fractal clusters, random coils
16(Debye equation), ellipsoidal particles, etc.
17
18The model works best for mass fractal systems characterized by Porod exponents
19between 5/3 and 3. It should not be used for surface fractal systems. Hammouda
20(2010) has pointed out a deficiency in the way this model handles the
21transitioning between the Guinier and Porod regimes and which can create
22artefacts that appear as kinks in the fitted model function.
23
24Also see the :ref:`guinier-porod` model.
25
26The empirical fit function is:
27
28.. math::
29
30    I(q) = \text{background}
31    + \sum_{i=1}^N \Bigl[
32        G_i \exp\Bigl(-\frac{q^2R_{gi}^2}{3}\Bigr)
33       + B_i \exp\Bigl(-\frac{q^2R_{g(i+1)}^2}{3}\Bigr)
34             \Bigl(\frac{1}{q_i^*}\Bigr)^{P_i} \Bigr]
35
36where
37
38.. math::
39
40    q_i^* = q \left[\operatorname{erf}
41            \left(\frac{q R_{gi}}{\sqrt{6}}\right)
42        \right]^{-3}
43
44
45For each level, the four parameters $G_i$, $R_{gi}$, $B_i$ and $P_i$ must
46be chosen.  Beaucage has an additional factor $k$ in the definition of
47$q_i^*$ which is ignored here.
48
49For example, to approximate the scattering from random coils (Debye equation),
50set $R_{gi}$ as the Guinier radius, $P_i = 2$, and $B_i = 2 G_i / R_{gi}$
51
52See the references for further information on choosing the parameters.
53
54For 2D data: The 2D scattering intensity is calculated in the same way as 1D,
55where the $q$ vector is defined as
56
57.. math::
58
59    q = \sqrt{q_x^2 + q_y^2}
60
61
62References
63----------
64
65.. [#] G Beaucage, *J. Appl. Cryst.*, 28 (1995) 717-728
66.. [#] G Beaucage, *J. Appl. Cryst.*, 29 (1996) 134-146
67.. [#] B Hammouda, *Analysis of the Beaucage model, J. Appl. Cryst.*, (2010), 43, 1474-1478
68
69Source
70------
71
72`unified_power_Rg.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/unified_power_Rg.py>`_
73
74Authorship and Verification
75----------------------------
76
77* **Author:**
78* **Last Modified by:**
79* **Last Reviewed by:**
80* **Source added by :** Steve King **Date:** March 25, 2019
81"""
82
83from __future__ import division
84
85import numpy as np
86from numpy import inf, exp, sqrt, errstate
87from scipy.special import erf, gamma
88
89category = "shape-independent"
90name = "unified_power_Rg"
91title = "Unified Power Rg"
92description = """
93        The Beaucage model employs the empirical multiple level unified
94        Exponential/Power-law fit method developed by G. Beaucage. Four functions
95        are included so that 1, 2, 3, or 4 levels can be used.
96        """
97
98# pylint: disable=bad-whitespace, line-too-long
99parameters = [
100    ["level",     "",     1,      [1, 6], "", "Level number"],
101    ["rg[level]", "Ang",  15.8,   [0, inf], "", "Radius of gyration"],
102    ["power[level]", "",  4,      [-inf, inf], "", "Power"],
103    ["B[level]",  "1/cm", 4.5e-6, [-inf, inf], "", ""],
104    ["G[level]",  "1/cm", 400,    [0, inf], "", ""],
105    ]
106# pylint: enable=bad-whitespace, line-too-long
107
108def Iq(q, level, rg, power, B, G):
109    """Return I(q) for unified power Rg model."""
110    level = int(level + 0.5)
111    if level == 0:
112        with errstate(divide='ignore'):
113            return 1./q
114
115    with errstate(divide='ignore', invalid='ignore'):
116        result = np.zeros(q.shape, 'd')
117        for i in range(level):
118            exp_now = exp(-(q*rg[i])**2/3.)
119            pow_now = (erf(q*rg[i]/sqrt(6.))**3/q)**power[i]
120            if i < level-1:
121                exp_next = exp(-(q*rg[i+1])**2/3.)
122            else:
123                exp_next = 1
124            result += G[i]*exp_now + B[i]*exp_next*pow_now
125
126    result[q == 0] = np.sum(G[:level])
127    return result
128
129Iq.vectorized = True
130
131def random():
132    """Return a random parameter set for the model."""
133    level = np.minimum(np.random.poisson(0.5) + 1, 6)
134    n = level
135    power = np.random.uniform(1.6, 3, n)
136    rg = 10**np.random.uniform(1, 5, n)
137    G = np.random.uniform(0.1, 10, n)**2 * 10**np.random.uniform(0.3, 3, n)
138    B = G * power / rg**power * gamma(power/2)
139    scale = 10**np.random.uniform(1, 4)
140    pars = dict(
141        #background=0,
142        scale=scale,
143        level=level,
144    )
145    pars.update(("power%d"%(k+1), v) for k, v in enumerate(power))
146    pars.update(("rg%d"%(k+1), v) for k, v in enumerate(rg))
147    pars.update(("B%d"%(k+1), v) for k, v in enumerate(B))
148    pars.update(("G%d"%(k+1), v) for k, v in enumerate(G))
149    return pars
150
151demo = dict(
152    level=2,
153    rg=[15.8, 21],
154    power=[4, 2],
155    B=[4.5e-6, 0.0006],
156    G=[400, 3],
157    scale=1.,
158    background=0.,
159)
Note: See TracBrowser for help on using the repository browser.