source: sasmodels/sasmodels/models/onion.py @ 7b68dc5

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 7b68dc5 was 7b68dc5, checked in by Paul Kienzle <pkienzle@…>, 6 years ago

core_multi_shell, onion: convert multiplicity to integer to avoid warnings

  • Property mode set to 100644
File size: 11.4 KB
Line 
1r"""
2This model provides the form factor, $P(q)$, for a multi-shell sphere where
3the scattering length density (SLD) of the each shell is described by an
4exponential, linear, or constant function. The form factor is normalized by
5the volume of the sphere where the SLD is not identical to the SLD of the
6solvent. We currently provide up to 9 shells with this model.
7
8NB: *radius* represents the core radius $r_0$ and
9*thickness[k]* represents the thickness of the shell, $r_{k+1} - r_k$.
10
11Definition
12----------
13
14The 1D scattering intensity is calculated in the following way
15
16.. math::
17
18    P(q) = [f]^2 / V_\text{particle}
19
20where
21
22.. math::
23    :nowrap:
24
25    \begin{align*}
26    f &= f_\text{core}
27            + \left(\sum_{\text{shell}=1}^N f_\text{shell}\right)
28            + f_\text{solvent}
29    \end{align*}
30
31The shells are spherically symmetric with particle density $\rho(r)$ and
32constant SLD within the core and solvent, so
33
34.. math::
35    :nowrap:
36
37    \begin{align*}
38    f_\text{core}
39        &= 4\pi\int_0^{r_\text{core}} \rho_\text{core}
40            \frac{\sin(qr)}{qr}\, r^2\,\mathrm{d}r
41        &= 3\rho_\text{core} V(r_\text{core})
42            \frac{j_1(qr_\text{core})}{qr_\text{core}} \\
43    f_\text{shell}
44        &= 4\pi\int_{r_{\text{shell}-1}}^{r_\text{shell}}
45            \rho_\text{shell}(r)\frac{\sin(qr)}{qr}\,r^2\,\mathrm{d}r \\
46    f_\text{solvent}
47        &= 4\pi\int_{r_N}^\infty
48            \rho_\text{solvent}\frac{\sin(qr)}{qr}\,r^2\,\mathrm{d}r
49        &= -3\rho_\text{solvent}V(r_N)\frac{j_1(q r_N)}{q r_N}
50    \end{align*}
51
52where the spherical bessel function $j_1$ is
53
54.. math::
55
56    j_1(x) = \frac{\sin(x)}{x^2} - \frac{\cos(x)}{x}
57
58and the volume is $V(r) = \frac{4\pi}{3}r^3$. The volume of the particle
59is determined by the radius of the outer shell, so $V_\text{particle} = V(r_N)$.
60
61Now lets consider the SLD of a shell defined by
62
63.. math::
64
65    \rho_\text{shell}(r) = \begin{cases}
66        B\exp\left(A(r-r_{\text{shell}-1})/\Delta t_\text{shell}\right)
67            + C & \mbox{for } A \neq 0 \\
68        \rho_\text{in} = \text{constant} & \mbox{for } A = 0
69    \end{cases}
70
71An example of a possible SLD profile is shown below where
72$\rho_\text{in}$ and $\Delta t_\text{shell}$ stand for the
73SLD of the inner side of the $k^\text{th}$ shell and the
74thickness of the $k^\text{th}$ shell in the equation above, respectively.
75
76For $A > 0$,
77
78.. math::
79
80    f_\text{shell} &= 4 \pi \int_{r_{\text{shell}-1}}^{r_\text{shell}}
81        \left[ B\exp
82            \left(A (r - r_{\text{shell}-1}) / \Delta t_\text{shell} \right) + C
83        \right] \frac{\sin(qr)}{qr}\,r^2\,\mathrm{d}r
84
85    &= 3BV(r_\text{shell}) e^A h(\alpha_\text{out},\beta_\text{out})
86        - 3BV(r_{\text{shell}-1}) h(\alpha_\text{in},\beta_\text{in})
87        + 3CV(r_{\text{shell}}) \frac{j_1(\beta_\text{out})}{\beta_\text{out}}
88        - 3CV(r_{\text{shell}-1}) \frac{j_1(\beta_\text{in})}{\beta_\text{in}}
89
90for
91
92.. math::
93    :nowrap:
94
95    \begin{align*}
96    B&=\frac{\rho_\text{out} - \rho_\text{in}}{e^A-1}
97         &C &= \frac{\rho_\text{in}e^A - \rho_\text{out}}{e^A-1} \\
98    \alpha_\text{in} &= A\frac{r_{\text{shell}-1}}{\Delta t_\text{shell}}
99         &\alpha_\text{out} &= A\frac{r_\text{shell}}{\Delta t_\text{shell}} \\
100    \beta_\text{in} &= qr_{\text{shell}-1}
101        &\beta_\text{out} &= qr_\text{shell} \\
102    \end{align*}
103
104where $h$ is
105
106 .. math::
107
108    h(x,y) = \frac{x \sin(y) - y\cos(y)}{(x^2+y^2)y}
109               - \frac{(x^2-y^2)\sin(y) - 2xy\cos(y)}{(x^2+y^2)^2y}
110
111
112For $A \sim 0$, e.g., $A = -0.0001$, this function converges to that of the
113linear SLD profile with
114$\rho_\text{shell}(r) \approx A(r-r_{\text{shell}-1})/\Delta t_\text{shell})+B$,
115so this case is equivalent to
116
117.. math::
118    :nowrap:
119
120    \begin{align*}
121    f_\text{shell}
122    &=
123      3 V(r_\text{shell}) \frac{\Delta\rho_\text{shell}}{\Delta t_\text{shell}}
124        \left[\frac{
125                2 \cos(qr_\text{out})
126                    + qr_\text{out} \sin(qr_\text{out})
127            }{
128                (qr_\text{out})^4
129            }\right] \\
130     &{}
131      -3 V(r_\text{shell}) \frac{\Delta\rho_\text{shell}}{\Delta t_\text{shell}}
132        \left[\frac{
133                    2\cos(qr_\text{in})
134                +qr_\text{in}\sin(qr_\text{in})
135            }{
136                (qr_\text{in})^4
137            }\right] \\
138    &{}
139      +3\rho_\text{out}V(r_\text{shell}) \frac{j_1(qr_\text{out})}{qr_\text{out}}
140      -3\rho_\text{in}V(r_{\text{shell}-1}) \frac{j_1(qr_\text{in})}{qr_\text{in}}
141    \end{align*}
142
143For $A = 0$, the exponential function has no dependence on the radius (so that
144$\rho_\text{out}$ is ignored this case) and becomes flat. We set the constant
145to $\rho_\text{in}$ for convenience, and thus the form factor contributed by
146the shells is
147
148.. math::
149
150    f_\text{shell} =
151        3\rho_\text{in}V(r_\text{shell})
152           \frac{j_1(qr_\text{out})}{qr_\text{out}}
153        - 3\rho_\text{in}V(r_{\text{shell}-1})
154            \frac{j_1(qr_\text{in})}{qr_\text{in}}
155
156.. figure:: img/onion_geometry.png
157
158    Example of an onion model profile.
159
160The 2D scattering intensity is the same as $P(q)$ above, regardless of the
161orientation of the $q$ vector which is defined as
162
163.. math::
164
165    q = \sqrt{q_x^2 + q_y^2}
166
167NB: The outer most radius is used as the effective radius for $S(q)$
168when $P(q) S(q)$ is applied.
169
170References
171----------
172
173L A Feigin and D I Svergun,
174*Structure Analysis by Small-Angle X-Ray and Neutron Scattering*,
175Plenum Press, New York, 1987.
176"""
177
178#
179# Give a polynomial $\rho(r) = Ar^3 + Br^2 + Cr + D$ for density,
180#
181# .. math::
182#
183#    f = 4 \pi \int_a^b \rho(r) \sin(qr)/(qr) \mathrm{d}r  = h(b) - h(a)
184#
185# where
186#
187# .. math::
188#
189#    h(r) = \frac{4 \pi}{q^6}\left[
190#        (q^3(4Ar^3 + 3Br^2 + 2Cr + D) - q(24Ar + 6B)) \sin(qr)
191#      - (q^4(Ar^4 + Br^3 + Cr^2 + Dr) - q^2(12Ar^2 + 6Br + 2C) + 24A) \cos(qr)
192#    \right]
193#
194# Use the monotonic spline to get the polynomial coefficients for each shell.
195#
196# Order 0
197#
198# .. math::
199#
200#    h(r) = \frac{4 \pi}{q^3} \left[
201#       - \cos(qr) (Ar) q
202#       + \sin(qr) (A)
203#    \right]
204#
205# Order 1
206#
207# .. math::
208#
209#   h(r) = \frac{4 \pi}{q^4} \left[
210#       - \cos(qr) ( Ar^2 + Br) q^2
211#       + \sin(qr) ( Ar   + B ) q
212#       + \cos(qr) (2A        )
213#   \right]
214#
215# Order 2
216#
217# .. math::
218#  h(r) = \frac{4 \pi}{q^5} \left[
219#        - \cos(qr) ( Ar^3 +  Br^2 + Cr) q^3
220#        + \sin(qr) (3Ar^2 + 2Br   + C ) q^2
221#        + \cos(qr) (6Ar   + 2B        ) q
222#        - \sin(qr) (6A                )
223#
224# Order 3
225#
226#    h(r) = \frac{4 \pi}{q^6}\left[
227#      - \cos(qr) (  Ar^4 +  Br^3 +  Cr^2 + Dr) q^4
228#      + \sin(qr) ( 4Ar^3 + 3Br^2 + 2Cr   + D ) q^3
229#      + \cos(qr) (12Ar^2 + 6Br   + 2C        ) q^2
230#      - \sin(qr) (24Ar   + 6B                ) q
231#      - \cos(qr) (24A                        )
232#    \right]
233#
234# Order p
235#
236#    h(r) = \frac{4 \pi}{q^{2}}
237#      \sum_{k=0}^p -\frac{d^k\cos(qr)}{dr^k} \frac{d^k r\rho(r)}{dr^k} (qr)^{-k}
238#
239# Given the equation
240#
241#    f = sum_(k=0)^(n-1) h_k(r_(k+1)) - h_k(r_k)
242#
243# we can rearrange the terms so that
244#
245#    f = sum_0^(n-1) h_k(r_(k+1)) - sum_0^(n-1) h_k(r_k)
246#      = sum_1^n h_(k-1)(r_k) - sum_0^(n-1) h_k(r_k)
247#      = h_(n-1)(r_n) - h_0(r_0) + sum_1^(n-1) [h_(k-1)(r_k) - h_k(r_k)]
248#      = h_(n-1)(r_n) - h_0(r_0) - sum_1^(n-1) h_(Delta k)(r_k)
249#
250# where
251#
252#    h_(Delta k)(r) = h(Delta rho_k, r)
253#
254# for
255#
256#    Delta rho_k = (A_k-A_(k-1)) r^p + (B_k-B_(k-1)) r^(p-1) + ...
257#
258# Using l'H\^opital's Rule 6 times on the order 3 polynomial,
259#
260#   lim_(q->0) h(r) = (140D r^3 + 180C r^4 + 144B r^5 + 120A r^6)/720
261#
262
263
264from __future__ import division
265
266import numpy as np
267from numpy import inf, nan
268from math import fabs, exp, expm1
269
270name = "onion"
271title = "Onion shell model with constant, linear or exponential density"
272
273description = """\
274Form factor of mutishells normalized by the volume. Here each shell is
275described by an exponential function;
276
277        I) For A_shell != 0,
278                f(r) = B*exp(A_shell*(r-r_in)/thick_shell)+C
279        where
280                B=(sld_out-sld_in)/(exp(A_shell)-1)
281                C=sld_in-B.
282        Note that in the above case, the function becomes a linear function
283        as A_shell --> 0+ or 0-.
284
285        II) For the exact point of A_shell == 0,
286                f(r) = sld_in ,i.e., it crosses over flat function
287        Note that the 'sld_out' becaomes NULL in this case.
288
289        background:background,
290        rad_core0: radius of sphere(core)
291        thick_shell#:the thickness of the shell#
292        sld_core0: the SLD of the sphere
293        sld_solv: the SLD of the solvent
294        sld_shell: the SLD of the shell#
295        A_shell#: the coefficient in the exponential function
296"""
297
298category = "shape:sphere"
299
300# TODO: n is a volume parameter that is not polydisperse
301
302# pylint: disable=bad-whitespace, line-too-long
303#   ["name", "units", default, [lower, upper], "type","description"],
304parameters = [
305    ["sld_core", "1e-6/Ang^2", 1.0, [-inf, inf], "sld", "Core scattering length density"],
306    ["radius_core", "Ang", 200., [0, inf], "volume", "Radius of the core"],
307    ["sld_solvent", "1e-6/Ang^2", 6.4, [-inf, inf], "sld", "Solvent scattering length density"],
308    ["n_shells", "", 1, [0, 10], "volume", "number of shells"],
309    ["sld_in[n_shells]", "1e-6/Ang^2", 1.7, [-inf, inf], "sld", "scattering length density at the inner radius of shell k"],
310    ["sld_out[n_shells]", "1e-6/Ang^2", 2.0, [-inf, inf], "sld", "scattering length density at the outer radius of shell k"],
311    ["thickness[n_shells]", "Ang", 40., [0, inf], "volume", "Thickness of shell k"],
312    ["A[n_shells]", "", 1.0, [-inf, inf], "", "Decay rate of shell k"],
313    ]
314# pylint: enable=bad-whitespace, line-too-long
315
316source = ["lib/sph_j1c.c", "onion.c"]
317single = False
318
319profile_axes = ['Radius (A)', 'SLD (1e-6/A^2)']
320def profile(sld_core, radius_core, sld_solvent, n_shells,
321            sld_in, sld_out, thickness, A):
322    """
323    Returns shape profile with x=radius, y=SLD.
324    """
325
326    total_radius = 1.25*(sum(thickness[:n_shells]) + radius_core + 1)
327    dz = total_radius/400  # 400 points for a smooth plot
328
329    z = []
330    rho = []
331
332    # add in the core
333    z.append(0)
334    rho.append(sld_core)
335    z.append(radius_core)
336    rho.append(sld_core)
337
338    # add in the shells
339    for k in range(int(n_shells)):
340        # Left side of each shells
341        z_current = z[-1]
342        z.append(z_current)
343        rho.append(sld_in[k])
344
345        if fabs(A[k]) < 1.0e-16:
346            # flat shell
347            z.append(z_current + thickness[k])
348            rho.append(sld_out[k])
349        else:
350            # exponential shell
351            # num_steps must be at least 1, so use floor()+1 rather than ceil
352            # to protect against a thickness0.
353            num_steps = np.floor(thickness[k]/dz) + 1
354            slope = (sld_out[k] - sld_in[k]) / expm1(A[k])
355            const = (sld_in[k] - slope)
356            for z_shell in np.linspace(0, thickness[k], num_steps+1):
357                z.append(z_current+z_shell)
358                rho.append(slope*exp(A[k]*z_shell/thickness[k]) + const)
359
360    # add in the solvent
361    z.append(z[-1])
362    rho.append(sld_solvent)
363    z.append(total_radius)
364    rho.append(sld_solvent)
365
366    return np.asarray(z), np.asarray(rho)
367
368def ER(core_radius, n, thickness):
369    """Effective radius"""
370    return np.sum(thickness[:int(n[0])], axis=0) + core_radius
371
372demo = {
373    "sld_solvent": 2.2,
374    "sld_core": 1.0,
375    "core_radius": 100,
376    "n_shells": 4,
377    "sld_in": [0.5, 1.5, 0.9, 2.0],
378    "sld_out": [nan, 0.9, 1.2, 1.6],
379    "thickness": [50, 75, 150, 75],
380    "A": [0, -1, 1e-4, 1],
381    # Could also specify them individually as
382    # "A1": 0, "A2": -1, "A3": 1e-4, "A4": 1,
383    #"core_radius_pd_n": 10,
384    #"core_radius_pd": 0.4,
385    #"thickness4_pd_n": 10,
386    #"thickness4_pd": 0.4,
387    }
388
Note: See TracBrowser for help on using the repository browser.