source: sasmodels/sasmodels/models/stacked_disks.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: 11.0 KB
RevLine 
[66d119f]1r"""
2Definition
3----------
4
[98ce141]5This model provides the form factor, $P(q)$, for stacked discs (tactoids)
6with a core/layer structure which is constructed itself as $P(q) S(Q)$
7multiplying a $P(q)$ for individual core/layer disks by a structure factor
8$S(q)$ proposed by Kratky and Porod in 1949\ [#CIT1949]_ assuming the next
9neighbor distance (d-spacing) in the stack of parallel discs obeys a Gaussian
10distribution. As such the normalization of this "composite" form factor is
11relative to the individual disk volume, not the volume of the stack of disks.
[2d81cfe]12This model is appropriate for example for non non exfoliated clay particles
13such as Laponite.
[98ce141]14
[5111921]15.. figure:: img/stacked_disks_geometry.png
[66d119f]16
[98ce141]17   Geometry of a single core/layer disk
18
[66d119f]19The scattered intensity $I(q)$ is calculated as
20
21.. math::
22
23    I(q) = N\int_{0}^{\pi /2}\left[ \Delta \rho_t
[98ce141]24    \left( V_t f_t(q,\alpha) - V_c f_c(q,\alpha)\right) + \Delta
25    \rho_c V_c f_c(q,\alpha)\right]^2 S(q,\alpha)\sin{\alpha}\ d\alpha
26    + \text{background}
[66d119f]27
28where the contrast
29
30.. math::
31
[40a87fa]32    \Delta \rho_i = \rho_i - \rho_\text{solvent}
[66d119f]33
[2d81cfe]34and $N$ is the number of individual (single) discs per unit volume, $\alpha$
35is the angle between the axis of the disc and $q$, and $V_t$ and $V_c$ are the
[98ce141]36total volume and the core volume of a single disc, respectively, and
[66d119f]37
38.. math::
39
[98ce141]40    f_t(q,\alpha) =
[40a87fa]41    \left(\frac{\sin(q(d+h)\cos{\alpha})}{q(d+h)\cos{\alpha}}\right)
[66d119f]42    \left(\frac{2J_1(qR\sin{\alpha})}{qR\sin{\alpha}} \right)
43
[98ce141]44    f_c(q,\alpha) =
[40a87fa]45    \left(\frac{\sin(qh)\cos{\alpha})}{qh\cos{\alpha}}\right)
46    \left(\frac{2J_1(qR\sin{\alpha})}{qR\sin{\alpha}}\right)
[66d119f]47
[a807206]48where $d$ = thickness of the layer (*thick_layer*),
49$2h$ = core thickness (*thick_core*), and $R$ = radius of the disc (*radius*).
[66d119f]50
51.. math::
52
[98ce141]53    S(q,\alpha) = 1 + \frac{1}{2}\sum_{k=1}^n(n-k)\cos{(kDq\cos{\alpha})}
54    \exp\left[ -k(q)^2(D\cos{\alpha}~\sigma_d)^2/2\right]
[66d119f]55
[40a87fa]56where $n$ is the total number of the disc stacked (*n_stacking*),
57$D = 2(d+h)$ is the next neighbor center-to-center distance (d-spacing),
[7c57861]58and $\sigma_d$ = the Gaussian standard deviation of the d-spacing (*sigma_d*).
[98ce141]59Note that $D\cos(\alpha)$ is the component of $D$ parallel to $q$ and the last
[f3073b0]60term in the equation above is effectively a Debye-Waller factor term.
[66d119f]61
62.. note::
[98ce141]63
64    1. Each assembly in the stack is layer/core/layer, so the spacing of the
[40a87fa]65    cores is core plus two layers. The 2nd virial coefficient of the cylinder
66    is calculated based on the *radius* and *length*
[a807206]67    = *n_stacking* * (*thick_core* + 2 * *thick_layer*)
[66d119f]68    values, and used as the effective radius for $S(Q)$ when $P(Q) * S(Q)$
69    is applied.
70
[98ce141]71    2. the resolution smearing calculation uses 76 Gaussian quadrature points
72    to properly smear the model since the function is HIGHLY oscillatory,
73    especially around the q-values that correspond to the repeat distance of
74    the layers.
75
[110f69c]762d scattering from oriented stacks is calculated in the same way as for
77cylinders, for further details of the calculation and angular dispersions
78see :ref:`orientation`.
[66d119f]79
[9802ab3]80.. figure:: img/cylinder_angle_definition.png
[66d119f]81
[eda8b30]82    Angles $\theta$ and $\phi$ orient the stack of discs relative
[110f69c]83    to the beam line coordinates, where the beam is along the $z$ axis.
84    Rotation $\theta$, initially in the $xz$ plane, is carried out first,
85    then rotation $\phi$ about the $z$ axis. Orientation distributions are
86    described as rotations about two perpendicular axes $\delta_1$ and
87    $\delta_2$ in the frame of the cylinder itself, which when
88    $\theta = \phi = 0$ are parallel to the $Y$ and $X$ axes.
[66d119f]89
90
[98ce141]91Our model is derived from the form factor calculations implemented in a
[07300ea]92c-library provided by the NIST Center for Neutron Research\ [#CIT_Kline]_
[66d119f]93
[e664a11]94References
95----------
[66d119f]96
[98ce141]97.. [#CIT1949] O Kratky and G Porod, *J. Colloid Science*, 4, (1949) 35
98.. [#CIT_Kline] S R Kline, *J Appl. Cryst.*, 39 (2006) 895
99.. [#] J S Higgins and H C Benoit, *Polymers and Neutron Scattering*,
100   Clarendon, Oxford, 1994
101.. [#] A Guinier and G Fournet, *Small-Angle Scattering of X-Rays*,
102   John Wiley and Sons, New York, 1955
[53215cf]103
[0507e09]104Source
105------
106
107`stacked_disks.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/stacked_disks.py>`_
108
109`stacked_disks.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/stacked_disks.c>`_
110
[98ce141]111Authorship and Verification
112----------------------------
[53215cf]113
[98ce141]114* **Author:** NIST IGOR/DANSE **Date:** pre 2010
[ef07e95]115* **Last Modified by:** Paul Butler and Paul Kienzle **Date:** November 26, 2016
116* **Last Reviewed by:** Paul Butler and Paul Kienzle **Date:** November 26, 2016
[0507e09]117* **Source added by :** Steve King **Date:** March 25, 2019
[66d119f]118"""
119
[2d81cfe]120import numpy as np
[0b56f38]121from numpy import inf, sin, cos, pi
[66d119f]122
123name = "stacked_disks"
[98ce141]124title = "Form factor for a stacked set of non exfoliated core/shell disks"
[66d119f]125description = """\
126    One layer of disk consists of a core, a top layer, and a bottom layer.
127    radius =  the radius of the disk
[a807206]128    thick_core = thickness of the core
129    thick_layer = thickness of a layer
[e664a11]130    sld_core = the SLD of the core
131    sld_layer = the SLD of the layers
[66d119f]132    n_stacking = the number of the disks
[7c57861]133    sigma_d =  Gaussian STD of d-spacing
[e664a11]134    sld_solvent = the SLD of the solvent
[66d119f]135    """
136category = "shape:cylinder"
137
138# pylint: disable=bad-whitespace, line-too-long
139#   ["name", "units", default, [lower, upper], "type","description"],
140parameters = [
[a807206]141    ["thick_core",  "Ang",        10.0, [0, inf],    "volume",      "Thickness of the core disk"],
142    ["thick_layer", "Ang",        10.0, [0, inf],    "volume",      "Thickness of layer each side of core"],
[66d119f]143    ["radius",      "Ang",        15.0, [0, inf],    "volume",      "Radius of the stacked disk"],
[edf1e8b]144    ["n_stacking",  "",            1.0, [1, inf],    "volume",      "Number of stacked layer/core/layer disks"],
[7c57861]145    ["sigma_d",     "Ang",         0,   [0, inf],    "",            "Sigma of nearest neighbor spacing"],
[42356c8]146    ["sld_core",    "1e-6/Ang^2",  4,   [-inf, inf], "sld",         "Core scattering length density"],
147    ["sld_layer",   "1e-6/Ang^2",  0.0, [-inf, inf], "sld",         "Layer scattering length density"],
148    ["sld_solvent", "1e-6/Ang^2",  5.0, [-inf, inf], "sld",         "Solvent scattering length density"],
[9b79f29]149    ["theta",       "degrees",     0,   [-360, 360], "orientation", "Orientation of the stacked disk axis w/respect incoming beam"],
150    ["phi",         "degrees",     0,   [-360, 360], "orientation", "Rotation about beam"],
[66d119f]151    ]
152# pylint: enable=bad-whitespace, line-too-long
153
[43b7eea]154source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "stacked_disks.c"]
[66d119f]155
[8f04da4]156def random():
[b297ba9]157    """Return a random parameter set for the model."""
[8f04da4]158    radius = 10**np.random.uniform(1, 4.7)
159    total_stack = 10**np.random.uniform(1, 4.7)
160    n_stacking = int(10**np.random.uniform(0, np.log10(total_stack)-1) + 0.5)
161    d = total_stack/n_stacking
162    thick_core = np.random.uniform(0, d-2)  # at least 1 A for each layer
163    thick_layer = (d - thick_core)/2
164    # Let polydispersity peak around 15%; 95% < 0.4; max=100%
165    sigma_d = d * np.random.beta(1.5, 7)
166    pars = dict(
167        thick_core=thick_core,
168        thick_layer=thick_layer,
169        radius=radius,
170        n_stacking=n_stacking,
171        sigma_d=sigma_d,
172    )
173    return pars
174
[66d119f]175demo = dict(background=0.001,
176            scale=0.01,
[a807206]177            thick_core=10.0,
178            thick_layer=10.0,
[66d119f]179            radius=15.0,
180            n_stacking=1,
[7c57861]181            sigma_d=0,
[e664a11]182            sld_core=4,
183            sld_layer=0.0,
184            sld_solvent=5.0,
[b7e8b94]185            theta=90,
[66d119f]186            phi=0)
[b7e8b94]187# After redefinition of spherical coordinates -
188# tests had in old coords theta=0, phi=0; new coords theta=90, phi=0
[0b56f38]189q = 0.1
190# april 6 2017, rkh added a 2d unit test, assume correct!
191qx = q*cos(pi/6.0)
192qy = q*sin(pi/6.0)
[b7e8b94]193# Accuracy tests based on content in test/utest_extra_models.py.
[f3073b0]194# Added 2 tests with n_stacked = 5 using SasView 3.1.2 - PDB;
195# for which alas q=0.001 values seem closer to n_stacked=1 not 5,
196# changed assuming my 4.1 code OK, RKH
[bb4b509]197tests = [
[b7e8b94]198    [{'thick_core': 10.0,
199      'thick_layer': 15.0,
200      'radius': 3000.0,
201      'n_stacking': 1.0,
202      'sigma_d': 0.0,
203      'sld_core': 4.0,
204      'sld_layer': -0.4,
[bb4b509]205      'sld_solvent': 5.0,
[b7e8b94]206      'theta': 90.0,
207      'phi': 0.0,
208      'scale': 0.01,
209      'background': 0.001,
210     }, 0.001, 5075.12],
211    [{'thick_core': 10.0,
212      'thick_layer': 15.0,
213      'radius': 3000.0,
214      'n_stacking': 5,
215      'sigma_d': 0.0,
216      'sld_core': 4.0,
217      'sld_layer': -0.4,
[bb4b509]218      'sld_solvent': 5.0,
[b7e8b94]219      'theta': 90.0,
220      'phi': 0.0,
221      'scale': 0.01,
222      'background': 0.001,
[2d81cfe]223      # n_stacking=1 not 5 ? slight change in value here 11jan2017,
224      # check other cpu types
225      #}, 0.001, 5065.12793824],
226      #}, 0.001, 5075.11570493],
[b7e8b94]227     }, 0.001, 25325.635693],
228    [{'thick_core': 10.0,
229      'thick_layer': 15.0,
[0b56f38]230      'radius': 100.0,
231      'n_stacking': 5,
232      'sigma_d': 0.0,
233      'sld_core': 4.0,
234      'sld_layer': -0.4,
[bb4b509]235      'sld_solvent': 5.0,
[0b56f38]236      'theta': 90.0,
237      'phi': 20.0,
238      'scale': 0.01,
[bb4b509]239      'background': 0.001,
240     }, (qx, qy), 0.0491167089952],
[0b56f38]241    [{'thick_core': 10.0,
242      'thick_layer': 15.0,
[b7e8b94]243      'radius': 3000.0,
244      'n_stacking': 5,
245      'sigma_d': 0.0,
246      'sld_core': 4.0,
247      'sld_layer': -0.4,
[bb4b509]248      'sld_solvent': 5.0,
[b7e8b94]249      'theta': 90.0,
250      'phi': 0.0,
251      'scale': 0.01,
252      'background': 0.001,
[2d81cfe]253      # n_stacking=1 not 5 ?  slight change in value here 11jan2017,
254      # check other cpu types
255      #}, 0.164, 0.0127673597265],
256      #}, 0.164, 0.01480812968],
[b7e8b94]257     }, 0.164, 0.0598367986509],
258
259    [{'thick_core': 10.0,
260      'thick_layer': 15.0,
261      'radius': 3000.0,
262      'n_stacking': 1.0,
263      'sigma_d': 0.0,
264      'sld_core': 4.0,
265      'sld_layer': -0.4,
[bb4b509]266      'sld_solvent': 5.0,
[b7e8b94]267      'theta': 90.0,
268      'phi': 0.0,
269      'scale': 0.01,
270      'background': 0.001,
[2d81cfe]271      # second test here was at q=90, changed it to q=5,
272      # note I(q) is then just value of flat background
[b7e8b94]273     }, [0.001, 5.0], [5075.12, 0.001]],
274
275    [{'thick_core': 10.0,
276      'thick_layer': 15.0,
277      'radius': 3000.0,
278      'n_stacking': 1.0,
279      'sigma_d': 0.0,
280      'sld_core': 4.0,
281      'sld_layer': -0.4,
[bb4b509]282      'sld_solvent': 5.0,
[b7e8b94]283      'theta': 90.0,
284      'phi': 0.0,
285      'scale': 0.01,
286      'background': 0.001,
287     }, ([0.4, 0.5]), [0.00105074, 0.00121761]],
[2d81cfe]288    #[{'thick_core': 10.0,
289    #  'thick_layer': 15.0,
290    #  'radius': 3000.0,
291    #  'n_stacking': 1.0,
292    #  'sigma_d': 0.0,
293    #  'sld_core': 4.0,
294    #  'sld_layer': -0.4,
295    #  'sld_solvent': 5.0,
296    #  'theta': 90.0,
297    #  'phi': 20.0,
298    #  'scale': 0.01,
299    #  'background': 0.001,
300    # 2017-05-18 PAK temporarily suppress output of qx,qy test; j1 is
301    #     not accurate for large qr
302    # }, (qx, qy), 0.0341738733124],
303    # }, (qx, qy), None],
[b7e8b94]304
305    [{'thick_core': 10.0,
306      'thick_layer': 15.0,
307      'radius': 3000.0,
308      'n_stacking': 1.0,
309      'sigma_d': 0.0,
310      'sld_core': 4.0,
311      'sld_layer': -0.4,
[bb4b509]312      'sld_solvent': 5.0,
[b7e8b94]313      'theta': 90.0,
314      'phi': 0.0,
315      'scale': 0.01,
316      'background': 0.001,
317     }, ([1.3, 1.57]), [0.0010039, 0.0010038]],
318    ]
[1b85b55]319# 11Jan2017   RKH checking unit test again, note they are all 1D, no 2D
Note: See TracBrowser for help on using the repository browser.