source: sasmodels/sasmodels/models/onion.py @ 71b751d

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 71b751d was 71b751d, checked in by Paul Kienzle <pkienzle@…>, 23 months ago

update remaining form factors to use Fq interface

  • 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 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    &= 3BV(r_\text{shell}) e^A h(\alpha_\text{out},\beta_\text{out})
85        - 3BV(r_{\text{shell}-1}) h(\alpha_\text{in},\beta_\text{in})
86        + 3CV(r_{\text{shell}}) \frac{j_1(\beta_\text{out})}{\beta_\text{out}}
87        - 3CV(r_{\text{shell}-1}) \frac{j_1(\beta_\text{in})}{\beta_\text{in}}
88
89for
90
91.. math::
92    :nowrap:
93
94    \begin{align*}
95    B&=\frac{\rho_\text{out} - \rho_\text{in}}{e^A-1}
96         & C &= \frac{\rho_\text{in}e^A - \rho_\text{out}}{e^A-1} \\
97    \alpha_\text{in} &= A\frac{r_{\text{shell}-1}}{\Delta t_\text{shell}}
98         & \alpha_\text{out} &= A\frac{r_\text{shell}}{\Delta t_\text{shell}} \\
99    \beta_\text{in} &= qr_{\text{shell}-1}
100        & \beta_\text{out} &= qr_\text{shell} \\
101    \end{align*}
102
103where $h$ is
104
105 .. math::
106
107    h(x,y) = \frac{x \sin(y) - y\cos(y)}{(x^2+y^2)y}
108               - \frac{(x^2-y^2)\sin(y) - 2xy\cos(y)}{(x^2+y^2)^2y}
109
110
111For $A \sim 0$, e.g., $A = -0.0001$, this function converges to that of the
112linear SLD profile with
113$\rho_\text{shell}(r) \approx A(r-r_{\text{shell}-1})/\Delta t_\text{shell})+B$,
114so this case is equivalent to
115
116.. math::
117    :nowrap:
118
119    \begin{align*}
120    f_\text{shell}
121    &=
122      3 V(r_\text{shell}) \frac{\Delta\rho_\text{shell}}{\Delta t_\text{shell}}
123        \left[\frac{
124                2 \cos(qr_\text{out})
125                    + qr_\text{out} \sin(qr_\text{out})
126            }{
127                (qr_\text{out})^4
128            }\right] \\
129     &{}
130      -3 V(r_\text{shell}) \frac{\Delta\rho_\text{shell}}{\Delta t_\text{shell}}
131        \left[\frac{
132                    2\cos(qr_\text{in})
133                +qr_\text{in}\sin(qr_\text{in})
134            }{
135                (qr_\text{in})^4
136            }\right] \\
137    &{}
138      +3\rho_\text{out}V(r_\text{shell}) \frac{j_1(qr_\text{out})}{qr_\text{out}}
139      -3\rho_\text{in}V(r_{\text{shell}-1}) \frac{j_1(qr_\text{in})}{qr_\text{in}}
140    \end{align*}
141
142For $A = 0$, the exponential function has no dependence on the radius (so that
143$\rho_\text{out}$ is ignored in this case) and becomes flat. We set the constant
144to $\rho_\text{in}$ for convenience, and thus the form factor contributed by
145the shells is
146
147.. math::
148
149    f_\text{shell} =
150        3\rho_\text{in}V(r_\text{shell})
151           \frac{j_1(qr_\text{out})}{qr_\text{out}}
152        - 3\rho_\text{in}V(r_{\text{shell}-1})
153            \frac{j_1(qr_\text{in})}{qr_\text{in}}
154
155.. figure:: img/onion_geometry.png
156
157    Example of an onion model profile.
158
159The 2D scattering intensity is the same as $P(q)$ above, regardless of the
160orientation of the $q$ vector which is defined as
161
162.. math::
163
164    q = \sqrt{q_x^2 + q_y^2}
165
166NB: The outer most radius is used as the effective radius for $S(q)$
167when $P(q) S(q)$ is applied.
168
169References
170----------
171
172L A Feigin and D I Svergun,
173*Structure Analysis by Small-Angle X-Ray and Neutron Scattering*,
174Plenum Press, New York, 1987.
175"""
176
177#
178# Give a polynomial $\rho(r) = Ar^3 + Br^2 + Cr + D$ for density,
179#
180# .. math::
181#
182#    f = 4 \pi \int_a^b \rho(r) \sin(qr)/(qr) \mathrm{d}r  = h(b) - h(a)
183#
184# where
185#
186# .. math::
187#
188#    h(r) = \frac{4 \pi}{q^6}\left[
189#        (q^3(4Ar^3 + 3Br^2 + 2Cr + D) - q(24Ar + 6B)) \sin(qr)
190#      - (q^4(Ar^4 + Br^3 + Cr^2 + Dr) - q^2(12Ar^2 + 6Br + 2C) + 24A) \cos(qr)
191#    \right]
192#
193# Use the monotonic spline to get the polynomial coefficients for each shell.
194#
195# Order 0
196#
197# .. math::
198#
199#    h(r) = \frac{4 \pi}{q^3} \left[
200#       - \cos(qr) (Ar) q
201#       + \sin(qr) (A)
202#    \right]
203#
204# Order 1
205#
206# .. math::
207#
208#   h(r) = \frac{4 \pi}{q^4} \left[
209#       - \cos(qr) ( Ar^2 + Br) q^2
210#       + \sin(qr) ( Ar   + B ) q
211#       + \cos(qr) (2A        )
212#   \right]
213#
214# Order 2
215#
216# .. math::
217#  h(r) = \frac{4 \pi}{q^5} \left[
218#        - \cos(qr) ( Ar^3 +  Br^2 + Cr) q^3
219#        + \sin(qr) (3Ar^2 + 2Br   + C ) q^2
220#        + \cos(qr) (6Ar   + 2B        ) q
221#        - \sin(qr) (6A                )
222#
223# Order 3
224#
225#    h(r) = \frac{4 \pi}{q^6}\left[
226#      - \cos(qr) (  Ar^4 +  Br^3 +  Cr^2 + Dr) q^4
227#      + \sin(qr) ( 4Ar^3 + 3Br^2 + 2Cr   + D ) q^3
228#      + \cos(qr) (12Ar^2 + 6Br   + 2C        ) q^2
229#      - \sin(qr) (24Ar   + 6B                ) q
230#      - \cos(qr) (24A                        )
231#    \right]
232#
233# Order p
234#
235#    h(r) = \frac{4 \pi}{q^{2}}
236#      \sum_{k=0}^p -\frac{d^k\cos(qr)}{dr^k} \frac{d^k r\rho(r)}{dr^k} (qr)^{-k}
237#
238# Given the equation
239#
240#    f = sum_(k=0)^(n-1) h_k(r_(k+1)) - h_k(r_k)
241#
242# we can rearrange the terms so that
243#
244#    f = sum_0^(n-1) h_k(r_(k+1)) - sum_0^(n-1) h_k(r_k)
245#      = sum_1^n h_(k-1)(r_k) - sum_0^(n-1) h_k(r_k)
246#      = h_(n-1)(r_n) - h_0(r_0) + sum_1^(n-1) [h_(k-1)(r_k) - h_k(r_k)]
247#      = h_(n-1)(r_n) - h_0(r_0) - sum_1^(n-1) h_(Delta k)(r_k)
248#
249# where
250#
251#    h_(Delta k)(r) = h(Delta rho_k, r)
252#
253# for
254#
255#    Delta rho_k = (A_k-A_(k-1)) r^p + (B_k-B_(k-1)) r^(p-1) + ...
256#
257# Using l'H\^opital's Rule 6 times on the order 3 polynomial,
258#
259#   lim_(q->0) h(r) = (140D r^3 + 180C r^4 + 144B r^5 + 120A r^6)/720
260#
261
262from __future__ import division
263
264from math import fabs, exp, expm1
265
266import numpy as np
267from numpy import inf, nan
268
269name = "onion"
270title = "Onion shell model with constant, linear or exponential density"
271
272description = """\
273Form factor of mutishells normalized by the volume. Here each shell is
274described by an exponential function;
275
276        I) For A_shell != 0,
277                f(r) = B*exp(A_shell*(r-r_in)/thick_shell)+C
278        where
279                B=(sld_out-sld_in)/(exp(A_shell)-1)
280                C=sld_in-B.
281        Note that in the above case, the function becomes a linear function
282        as A_shell --> 0+ or 0-.
283
284        II) For the exact point of A_shell == 0,
285                f(r) = sld_in ,i.e., it crosses over flat function
286        Note that the 'sld_out' becaomes NULL in this case.
287
288        background:background,
289        rad_core0: radius of sphere(core)
290        thick_shell#:the thickness of the shell#
291        sld_core0: the SLD of the sphere
292        sld_solv: the SLD of the solvent
293        sld_shell: the SLD of the shell#
294        A_shell#: the coefficient in the exponential function
295"""
296
297category = "shape:sphere"
298
299# TODO: n is a volume parameter that is not polydisperse
300
301# pylint: disable=bad-whitespace, line-too-long
302#   ["name", "units", default, [lower, upper], "type","description"],
303parameters = [
304    ["sld_core", "1e-6/Ang^2", 1.0, [-inf, inf], "sld", "Core scattering length density"],
305    ["radius_core", "Ang", 200., [0, inf], "volume", "Radius of the core"],
306    ["sld_solvent", "1e-6/Ang^2", 6.4, [-inf, inf], "sld", "Solvent scattering length density"],
307    ["n_shells", "", 1, [0, 10], "volume", "number of shells"],
308    ["sld_in[n_shells]", "1e-6/Ang^2", 1.7, [-inf, inf], "sld", "scattering length density at the inner radius of shell k"],
309    ["sld_out[n_shells]", "1e-6/Ang^2", 2.0, [-inf, inf], "sld", "scattering length density at the outer radius of shell k"],
310    ["thickness[n_shells]", "Ang", 40., [0, inf], "volume", "Thickness of shell k"],
311    ["A[n_shells]", "", 1.0, [-inf, inf], "", "Decay rate of shell k"],
312    ]
313# pylint: enable=bad-whitespace, line-too-long
314
315source = ["lib/sas_3j1x_x.c", "onion.c"]
316single = False
317have_Fq = True
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    n_shells = int(n_shells+0.5)
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_in[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(radius_core, n_shells, thickness):
369    """Effective radius"""
370    n = int(n_shells[0]+0.5)
371    return np.sum(thickness[:n], axis=0) + radius_core
372
373demo = {
374    "sld_solvent": 2.2,
375    "sld_core": 1.0,
376    "radius_core": 100,
377    "n_shells": 4,
378    "sld_in": [0.5, 1.5, 0.9, 2.0],
379    "sld_out": [nan, 0.9, 1.2, 1.6],
380    "thickness": [50, 75, 150, 75],
381    "A": [0, -1, 1e-4, 1],
382    # Could also specify them individually as
383    # "A1": 0, "A2": -1, "A3": 1e-4, "A4": 1,
384    #"radius_core_pd_n": 10,
385    #"radius_core_pd": 0.4,
386    #"thickness4_pd_n": 10,
387    #"thickness4_pd": 0.4,
388    }
Note: See TracBrowser for help on using the repository browser.