source: sasmodels/sasmodels/models/flexible_cylinder.py @ e220acc

ticket-1257-vesicle-productticket_1156ticket_822_more_unit_tests
Last change on this file since e220acc was db1d9d5, checked in by Paul Kienzle <pkienzle@…>, 5 years ago

merge with master

  • Property mode set to 100644
File size: 7.0 KB
Line 
1r"""
2This model provides the form factor, $P(q)$, for a flexible cylinder
3where the form factor is normalized by the volume of the cylinder.
4**Inter-cylinder interactions are NOT provided for.**
5
6.. math::
7
8    P(q) = \text{scale} \left<F^2\right>/V + \text{background}
9
10where the averaging $\left<\ldots\right>$ is applied only for the 1D
11calculation
12
13The 2D scattering intensity is the same as 1D, regardless of the orientation of
14the q vector which is defined as
15
16.. math::
17
18    q = \sqrt{q_x^2 + q_y^2}
19
20Definitions
21-----------
22
23.. figure:: img/flexible_cylinder_geometry.jpg
24
25
26The chain of contour length, $L$, (the total length) can be described as a
27chain of some number of locally stiff segments of length $l_p$, the persistence
28length (the length along the cylinder over which the flexible cylinder can be
29considered a rigid rod).
30The Kuhn length $(b = 2*l_p)$ is also used to describe the stiffness of a chain.
31
32In the parameters, the sld and sld\_solvent represent the SLD of the cylinder
33and solvent respectively.
34
35Our model uses the form factor calculations in reference [1] as implemented in a
36c-library provided by the NIST Center for Neutron Research (Kline, 2006). This states:
37
38    'Method 3 With Excluded Volume' is used.
39    The model is a parametrization of simulations of a discrete representation
40    of the worm-like chain model of Kratky and Porod applied in the
41    pseudocontinuous limit.
42    See equations (13,26-27) in the original reference for the details.
43
44.. note::
45
46    There are several typos in the original reference that have been corrected
47    by WRC [2]. Details of the corrections are in the reference below. Most notably
48
49    - Equation (13): the term $(1 - w(QR))$ should swap position with $w(QR)$
50
51    - Equations (23) and (24) are incorrect; WRC has entered these into
52      Mathematica and solved analytically. The results were then converted to
53      code.
54
55    - Equation (27) should be $q0 = max(a3/(Rg^2)^{1/2},3)$ instead of
56      $max(a3*b(Rg^2)^{1/2},3)$
57
58    - The scattering function is negative for a range of parameter values and
59      q-values that are experimentally accessible. A correction function has been
60      added to give the proper behavior.
61
62
63**This is a model with complex behaviour depending on the ratio of** $L/b$ **and the
64reader is strongly encouraged to read reference [1] before use.**
65
66References
67----------
68
69.. [#] J S Pedersen and P Schurtenberger. *Scattering functions of semiflexible polymers with and without excluded volume effects.* Macromolecules, 29 (1996) 7602-7612
70
71Correction of the formula can be found in
72
73.. [#] W R Chen, P D Butler and L J Magid, *Incorporating Intermicellar Interactions in the Fitting of SANS Data from Cationic Wormlike Micelles.* Langmuir, 22(15) 2006 6539-6548
74
75Source
76------
77
78`flexible_cylinder.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/flexible_cylinder.py>`_
79
80`flexible_cylinder.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/flexible_cylinder.c>`_
81
82`wrc_cyl.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/lib/wrc_cyl.c>`_
83
84Authorship and Verification
85----------------------------
86
87* **Author:**
88* **Last Modified by:**
89* **Last Reviewed by:** Steve King **Date:** March 26, 2019
90* **Source added by :** Steve King **Date:** March 25, 2019
91"""
92
93import numpy as np
94from numpy import inf
95
96name = "flexible_cylinder"
97title = "Flexible cylinder where the form factor is normalized by the volume " \
98        "of the cylinder."
99description = """Note : scale and contrast = (sld - sld_solvent) are both
100                multiplicative factors in the model and are perfectly
101                correlated. One or both of these parameters must be held fixed
102                during model fitting.
103              """
104
105category = "shape:cylinder"
106single = False  # double precision only!
107
108# pylint: disable=bad-whitespace, line-too-long
109#             ["name", "units", default, [lower, upper], "type", "description"],
110parameters = [
111    ["length",      "Ang",       1000.0, [0, inf],    "volume", "Length of the flexible cylinder"],
112    ["kuhn_length", "Ang",        100.0, [0, inf],    "volume", "Kuhn length of the flexible cylinder"],
113    ["radius",      "Ang",         20.0, [0, inf],    "volume", "Radius of the flexible cylinder"],
114    ["sld",         "1e-6/Ang^2",   1.0, [-inf, inf], "sld",    "Cylinder scattering length density"],
115    ["sld_solvent", "1e-6/Ang^2",   6.3, [-inf, inf], "sld",    "Solvent scattering length density"],
116    ]
117# pylint: enable=bad-whitespace, line-too-long
118source = ["lib/polevl.c", "lib/sas_J1.c", "lib/wrc_cyl.c", "flexible_cylinder.c"]
119
120def random():
121    """Return a random parameter set for the model."""
122    length = 10**np.random.uniform(2, 6)
123    radius = 10**np.random.uniform(1, 3)
124    kuhn_length = 10**np.random.uniform(-2, 0)*length
125    pars = dict(
126        length=length,
127        radius=radius,
128        kuhn_length=kuhn_length,
129    )
130    return pars
131
132tests = [
133    # Accuracy tests based on content in test/utest_other_models.py
134    [{'length':     1000.0,  # test T1
135      'kuhn_length': 100.0,
136      'radius':       20.0,
137      'sld':           1.0,
138      'sld_solvent':   6.3,
139      'background':    0.0001,
140     }, 0.001, 3509.2187],
141
142    # Additional tests with larger range of parameters
143    [{'length':    1000.0,  # test T2
144      'kuhn_length': 100.0,
145      'radius':       20.0,
146      'sld':           1.0,
147      'sld_solvent':   6.3,
148      'background':    0.0001,
149     }, 1.0, 0.000595345],
150    [{'length':        10.0,  # test T3
151      'kuhn_length': 800.0,
152      'radius':        2.0,
153      'sld':           6.0,
154      'sld_solvent':  12.3,
155      'background':    0.001,
156     }, 0.1, 1.55228],
157    [{'length':        100.0,  # test T4
158      'kuhn_length': 800.0,
159      'radius':       50.0,
160      'sld':           0.1,
161      'sld_solvent':   5.1,
162      'background':    0.0,
163     }, 1.0, 0.000938456]
164    ]
165
166# There are a few branches in the code that ought to have test values:
167#
168# For length > 4 * kuhn_length
169#        if length > 10 * kuhn_length then C is scaled by 3.06 (L/b)^(-0.44)
170#        q*kuhn_length <= 3.1  => Sexv_new
171#           dS/dQ < 0 has different behaviour from dS/dQ >= 0
172#  T2    q*kuhn_length > 3.1   => a_long
173#
174# For length <= 4 * kuhn_length
175#        q*kuhn_length <= max(1.9/Rg_short, 3.0)  => Sdebye((q*Rg)^2)
176#           q*Rg < 0.5 uses Pade approx, q*Rg > 1.0 uses math lib
177#  T3,T4 q*kuhn_length > max(1.9/Rg_short, 3.0)   => a_short
178#
179# Note that the transitions between branches may be abrupt.  You can see a
180# several percent change around length=10*kuhn_length and length=4*kuhn_length
181# using the following:
182#
183#    sascomp flexible_cylinder -calc=double -sets=10 length=10*kuhn_length,10.000001*kuhn_length
184#    sascomp flexible_cylinder -calc=double -sets=10 length=4*kuhn_length,4.000001*kuhn_length
185#
186# The transition between low q and high q around q*kuhn_length = 3 seems
187# to be good to 4 digits or better.  This was tested by computing the value
188# on each branches near the transition point and reporting the relative error
189# for kuhn lengths of 10, 100 and 1000 and a variety of length:kuhn_length
190# ratios.
Note: See TracBrowser for help on using the repository browser.