source: sasmodels/sasmodels/models/onion.py @ 99658f6

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 99658f6 was ee60aa7, checked in by Paul Kienzle <pkienzle@…>, 6 years ago

clean up effective radius functions; improve mono_gauss_coil accuracy; start moving VR into C

  • Property mode set to 100644
File size: 11.3 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
318effective_radius_type = ["outer radius"]
319
320profile_axes = ['Radius (A)', 'SLD (1e-6/A^2)']
321def profile(sld_core, radius_core, sld_solvent, n_shells,
322            sld_in, sld_out, thickness, A):
323    """
324    Returns shape profile with x=radius, y=SLD.
325    """
326    n_shells = int(n_shells+0.5)
327    total_radius = 1.25*(sum(thickness[:n_shells]) + radius_core + 1)
328    dz = total_radius/400  # 400 points for a smooth plot
329
330    z = []
331    rho = []
332
333    # add in the core
334    z.append(0)
335    rho.append(sld_core)
336    z.append(radius_core)
337    rho.append(sld_core)
338
339    # add in the shells
340    for k in range(int(n_shells)):
341        # Left side of each shells
342        z_current = z[-1]
343        z.append(z_current)
344        rho.append(sld_in[k])
345
346        if fabs(A[k]) < 1.0e-16:
347            # flat shell
348            z.append(z_current + thickness[k])
349            rho.append(sld_in[k])
350        else:
351            # exponential shell
352            # num_steps must be at least 1, so use floor()+1 rather than ceil
353            # to protect against a thickness0.
354            num_steps = np.floor(thickness[k]/dz) + 1
355            slope = (sld_out[k] - sld_in[k]) / expm1(A[k])
356            const = (sld_in[k] - slope)
357            for z_shell in np.linspace(0, thickness[k], num_steps+1):
358                z.append(z_current+z_shell)
359                rho.append(slope*exp(A[k]*z_shell/thickness[k]) + const)
360
361    # add in the solvent
362    z.append(z[-1])
363    rho.append(sld_solvent)
364    z.append(total_radius)
365    rho.append(sld_solvent)
366
367    return np.asarray(z), np.asarray(rho)
368
369demo = {
370    "sld_solvent": 2.2,
371    "sld_core": 1.0,
372    "radius_core": 100,
373    "n_shells": 4,
374    "sld_in": [0.5, 1.5, 0.9, 2.0],
375    "sld_out": [nan, 0.9, 1.2, 1.6],
376    "thickness": [50, 75, 150, 75],
377    "A": [0, -1, 1e-4, 1],
378    # Could also specify them individually as
379    # "A1": 0, "A2": -1, "A3": 1e-4, "A4": 1,
380    #"radius_core_pd_n": 10,
381    #"radius_core_pd": 0.4,
382    #"thickness4_pd_n": 10,
383    #"thickness4_pd": 0.4,
384    }
Note: See TracBrowser for help on using the repository browser.