source: sasmodels/sasmodels/models/polymer_excl_volume.py @ 168052c

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 168052c was 168052c, checked in by piotr, 8 years ago

pylint prettification

  • Property mode set to 100644
File size: 5.2 KB
Line 
1r"""
2This model describes the scattering from polymer chains subject to excluded
3volume effects and has been used as a template for describing mass fractals.
4
5Definition
6----------
7
8The form factor was originally presented in the following integral form
9(Benoit, 1957)
10
11.. math::
12
13    P(Q)=2\int_0^{1}dx(1-x)exp\left[-\frac{Q^2a^2}{6}n^{2v}x^{2v}\right]
14
15where $\nu$ is the excluded volume parameter
16(which is related to the Porod exponent $m$ as $\nu=1/m$ ),
17$a$ is the statistical segment length of the polymer chain,
18and $n$ is the degree of polymerization.
19This integral was later put into an almost analytical form as follows
20(Hammouda, 1993)
21
22.. math::
23
24    P(Q)=\frac{1}{\nu U^{1/2\nu}}\gamma\left(\frac{1}{2\nu},U\right) -
25    \frac{1}{\nu U^{1/\nu}}\gamma\left(\frac{1}{\nu},U\right)
26
27where $\gamma(x,U)$ is the incomplete gamma function
28
29.. math::
30
31    \gamma(x,U)=\int_0^{U}dt\ exp(-t)t^{x-1}
32
33and the variable $U$ is given in terms of the scattering vector $Q$ as
34
35.. math::
36
37    U=\frac{Q^2a^2n^{2\nu}}{6} = \frac{Q^2R_{g}^2(2\nu+1)(2\nu+2)}{6}
38
39The square of the radius-of-gyration is defined as
40
41.. math::
42
43    R_{g}^2 = \frac{a^2n^{2\nu}}{(2\nu+1)(2\nu+2)}
44
45Note that this model applies only in the mass fractal range (ie, $5/3<=m<=3$ )
46and **does not apply** to surface fractals ( $3<m<=4$ ).
47It also does not reproduce the rigid rod limit (m=1) because it assumes chain
48flexibility from the outset. It may cover a portion of the semi-flexible chain
49range ( $1<m<5/3$ ).
50
51A low-Q expansion yields the Guinier form and a high-Q expansion yields the
52Porod form which is given by
53
54.. math::
55
56    P(Q\rightarrow \infty) = \frac{1}{\nu U^{1/2\nu}}\Gamma\left(
57    \frac{1}{2\nu}\right) - \frac{1}{\nu U^{1/\nu}}\Gamma\left(
58    \frac{1}{\nu}\right)
59
60Here $\Gamma(x) = \gamma(x,\infty)$ is the gamma function.
61
62The asymptotic limit is dominated by the first term
63
64.. math::
65
66    P(Q\rightarrow \infty) \sim \frac{1}{\nu U^{1/2\nu}}\Gamma\left(\frac{1}{2\nu}\right) =
67    \frac{m}{\left(QR_{g}\right)^m}\left[\frac{6}{(2\nu +1)(2\nu +2)} \right]^{m/2}
68    \Gamma (m/2)
69
70The special case when $\nu=0.5$ (or $m=2/\nu=2$ ) corresponds to Gaussian chains for
71which the form factor is given by the familiar Debye function.
72
73.. math::
74
75    P(Q) = \frac{2}{Q^4R_{g}^4} \left[exp(-Q^2R_{g}^2) - 1 + Q^2R_{g}^2 \right]
76
77For 2D data: The 2D scattering intensity is calculated in the same way as 1D,
78where the $q$ vector is defined as
79
80.. math::
81
82    q = \sqrt{q_x^2 + q_y^2}
83
84This example dataset is produced using 200 data points, $qmin=0.001Ang^{-1}$,
85$qmax=0.2Ang^{-1}$ and the default values
86
87.. figure:: img/polymer_excl_volume_1d.jpg
88
89    1D plot using the default values (w/500 data point).
90
91
92References
93----------
94
95H Benoit, *Comptes Rendus*, 245 (1957) 2244-2247
96
97B Hammouda, *SANS from Homogeneous Polymer Mixtures - A Unified Overview,
98Advances in Polym. Sci.* 106(1993) 87-133
99
100"""
101
102from math import sqrt
103from numpy import inf, power
104from scipy.special import gammainc, gamma
105
106name = "polymer_excl_volume"
107title = "Polymer Excluded Volume model"
108description = """Compute the scattering intensity from polymers with excluded
109                volume effects.
110                rg:         radius of gyration
111                porod_exp:  Porod exponent
112              """
113category = "shape-independent"
114
115# pylint: disable=bad-whitespace, line-too-long
116#             ["name", "units", default, [lower, upper], "type", "description"],
117parameters = [["rg",        "Ang", 60.0, [0, inf],    "", "Radius of Gyration"],
118              ["porod_exp", "",     3.0, [-inf, inf], "", "Porod exponent"],
119             ]
120# pylint: enable=bad-whitespace, line-too-long
121
122
123def Iq(q,
124       rg=60.0,
125       porod_exp=3.0):
126    """
127    :param q:         Input q-value (float or [float, float])
128    :param rg:        Radius of gyration
129    :param porod_exp: Porod exponent
130    :return:          Calculated intensity
131    """
132    nu = 1.0/porod_exp
133    u = q*q*rg*rg*(2.0*nu+1.0) * (2.0*nu+2.0)/6.0
134    o2nu = 1.0/(2.0*nu)
135
136    intensity = ((1.0/(nu*power(u, o2nu))) *
137                 (gamma(o2nu)*gammainc(o2nu, u) -
138                  1.0/power(u, o2nu) * gamma(porod_exp) *
139                  gammainc(porod_exp, u))) * (q > 0) + 1.0*(q <= 0)
140
141    return intensity
142
143Iq.vectorized = True  # Iq accepts an array of q values
144
145
146def Iqxy(qx, qy, *args):
147    """
148    :param qx:   Input q_x-value
149    :param qy:   Input q_y-value
150    :param args: Remaining arguments
151    :return:     2D-Intensity
152    """
153
154    return Iq(sqrt(qx**2 + qy**2), *args)
155
156Iqxy.vectorized = True  # Iqxy accepts an array of qx, qy values
157
158
159demo = dict(scale=1, background=0.0,
160            rg=60.0,
161            porod_exp=3.0)
162
163oldname = "PolymerExclVolume"
164oldpars = dict(background='background', scale='scale',
165               rg='rg',
166               porod_exp='m')
167
168tests = [
169    # Accuracy tests based on content in test/polyexclvol_default_igor.txt
170    [{'rg': 60, 'porod_exp': 3.0}, 0.001, 0.998801],
171    [{'rg': 60, 'porod_exp': 3.0}, 0.105363, 0.0162751],
172    [{'rg': 60, 'porod_exp': 3.0}, 0.665075, 6.56261e-05],
173
174    # Additional tests with larger range of parameters
175    [{'rg': 10, 'porod_exp': 4.0}, 0.1, 0.723436675809],
176    [{'rg': 2.2, 'porod_exp': 22.0, 'background': 100.0}, 5.0, 100.0],
177    [{'rg': 1.1, 'porod_exp': 1, 'background': 10.0, 'scale': 1.25},
178     20000., 10.0000712097]
179    ]
Note: See TracBrowser for help on using the repository browser.