# source:sasmodels/sasmodels/models/onion.py@db1d9d5

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

merge with master

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