source: sasmodels/sasmodels/models/onion.py @ a34b811

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

use radius_effective/radius_effective_mode/radius_effective_modes consistently throughout the code

  • Property mode set to 100644
File size: 11.8 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
172.. [#] L A Feigin and D I Svergun, *Structure Analysis by Small-Angle X-Ray and Neutron Scattering*, Plenum Press, New York, 1987.
173
174Source
175------
176
177`onion.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/onion.py>`_
178
179`onion.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/onion.c>`_
180
181Authorship and Verification
182----------------------------
183
184* **Author:**
185* **Last Modified by:**
186* **Last Reviewed by:**
187* **Source added by :** Steve King **Date:** March 25, 2019
188"""
189
190#
191# Give a polynomial $\rho(r) = Ar^3 + Br^2 + Cr + D$ for density,
192#
193# .. math::
194#
195#    f = 4 \pi \int_a^b \rho(r) \sin(qr)/(qr) \mathrm{d}r  = h(b) - h(a)
196#
197# where
198#
199# .. math::
200#
201#    h(r) = \frac{4 \pi}{q^6}\left[
202#        (q^3(4Ar^3 + 3Br^2 + 2Cr + D) - q(24Ar + 6B)) \sin(qr)
203#      - (q^4(Ar^4 + Br^3 + Cr^2 + Dr) - q^2(12Ar^2 + 6Br + 2C) + 24A) \cos(qr)
204#    \right]
205#
206# Use the monotonic spline to get the polynomial coefficients for each shell.
207#
208# Order 0
209#
210# .. math::
211#
212#    h(r) = \frac{4 \pi}{q^3} \left[
213#       - \cos(qr) (Ar) q
214#       + \sin(qr) (A)
215#    \right]
216#
217# Order 1
218#
219# .. math::
220#
221#   h(r) = \frac{4 \pi}{q^4} \left[
222#       - \cos(qr) ( Ar^2 + Br) q^2
223#       + \sin(qr) ( Ar   + B ) q
224#       + \cos(qr) (2A        )
225#   \right]
226#
227# Order 2
228#
229# .. math::
230#  h(r) = \frac{4 \pi}{q^5} \left[
231#        - \cos(qr) ( Ar^3 +  Br^2 + Cr) q^3
232#        + \sin(qr) (3Ar^2 + 2Br   + C ) q^2
233#        + \cos(qr) (6Ar   + 2B        ) q
234#        - \sin(qr) (6A                )
235#
236# Order 3
237#
238#    h(r) = \frac{4 \pi}{q^6}\left[
239#      - \cos(qr) (  Ar^4 +  Br^3 +  Cr^2 + Dr) q^4
240#      + \sin(qr) ( 4Ar^3 + 3Br^2 + 2Cr   + D ) q^3
241#      + \cos(qr) (12Ar^2 + 6Br   + 2C        ) q^2
242#      - \sin(qr) (24Ar   + 6B                ) q
243#      - \cos(qr) (24A                        )
244#    \right]
245#
246# Order p
247#
248#    h(r) = \frac{4 \pi}{q^{2}}
249#      \sum_{k=0}^p -\frac{d^k\cos(qr)}{dr^k} \frac{d^k r\rho(r)}{dr^k} (qr)^{-k}
250#
251# Given the equation
252#
253#    f = sum_(k=0)^(n-1) h_k(r_(k+1)) - h_k(r_k)
254#
255# we can rearrange the terms so that
256#
257#    f = sum_0^(n-1) h_k(r_(k+1)) - sum_0^(n-1) h_k(r_k)
258#      = sum_1^n h_(k-1)(r_k) - sum_0^(n-1) h_k(r_k)
259#      = h_(n-1)(r_n) - h_0(r_0) + sum_1^(n-1) [h_(k-1)(r_k) - h_k(r_k)]
260#      = h_(n-1)(r_n) - h_0(r_0) - sum_1^(n-1) h_(Delta k)(r_k)
261#
262# where
263#
264#    h_(Delta k)(r) = h(Delta rho_k, r)
265#
266# for
267#
268#    Delta rho_k = (A_k-A_(k-1)) r^p + (B_k-B_(k-1)) r^(p-1) + ...
269#
270# Using l'H\^opital's Rule 6 times on the order 3 polynomial,
271#
272#   lim_(q->0) h(r) = (140D r^3 + 180C r^4 + 144B r^5 + 120A r^6)/720
273#
274
275from __future__ import division
276
277from math import fabs, exp, expm1
278
279import numpy as np
280from numpy import inf, nan
281
282name = "onion"
283title = "Onion shell model with constant, linear or exponential density"
284
285description = """\
286Form factor of mutishells normalized by the volume. Here each shell is
287described by an exponential function;
288
289        I) For A_shell != 0,
290                f(r) = B*exp(A_shell*(r-r_in)/thick_shell)+C
291        where
292                B=(sld_out-sld_in)/(exp(A_shell)-1)
293                C=sld_in-B.
294        Note that in the above case, the function becomes a linear function
295        as A_shell --> 0+ or 0-.
296
297        II) For the exact point of A_shell == 0,
298                f(r) = sld_in ,i.e., it crosses over flat function
299        Note that the 'sld_out' becaomes NULL in this case.
300
301        background:background,
302        rad_core0: radius of sphere(core)
303        thick_shell#:the thickness of the shell#
304        sld_core0: the SLD of the sphere
305        sld_solv: the SLD of the solvent
306        sld_shell: the SLD of the shell#
307        A_shell#: the coefficient in the exponential function
308"""
309
310category = "shape:sphere"
311
312# TODO: n is a volume parameter that is not polydisperse
313
314# pylint: disable=bad-whitespace, line-too-long
315#   ["name", "units", default, [lower, upper], "type","description"],
316parameters = [
317    ["sld_core", "1e-6/Ang^2", 1.0, [-inf, inf], "sld", "Core scattering length density"],
318    ["radius_core", "Ang", 200., [0, inf], "volume", "Radius of the core"],
319    ["sld_solvent", "1e-6/Ang^2", 6.4, [-inf, inf], "sld", "Solvent scattering length density"],
320    ["n_shells", "", 1, [0, 10], "volume", "number of shells"],
321    ["sld_in[n_shells]", "1e-6/Ang^2", 1.7, [-inf, inf], "sld", "scattering length density at the inner radius of shell k"],
322    ["sld_out[n_shells]", "1e-6/Ang^2", 2.0, [-inf, inf], "sld", "scattering length density at the outer radius of shell k"],
323    ["thickness[n_shells]", "Ang", 40., [0, inf], "volume", "Thickness of shell k"],
324    ["A[n_shells]", "", 1.0, [-inf, inf], "", "Decay rate of shell k"],
325    ]
326# pylint: enable=bad-whitespace, line-too-long
327
328source = ["lib/sas_3j1x_x.c", "onion.c"]
329single = False
330have_Fq = True
331radius_effective_modes = ["outer radius"]
332
333profile_axes = ['Radius (A)', 'SLD (1e-6/A^2)']
334def profile(sld_core, radius_core, sld_solvent, n_shells,
335            sld_in, sld_out, thickness, A):
336    """
337    Returns shape profile with x=radius, y=SLD.
338    """
339    n_shells = int(n_shells+0.5)
340    total_radius = 1.25*(sum(thickness[:n_shells]) + radius_core + 1)
341    dz = total_radius/400  # 400 points for a smooth plot
342
343    z = []
344    rho = []
345
346    # add in the core
347    z.append(0)
348    rho.append(sld_core)
349    z.append(radius_core)
350    rho.append(sld_core)
351
352    # add in the shells
353    for k in range(int(n_shells)):
354        # Left side of each shells
355        z_current = z[-1]
356        z.append(z_current)
357        rho.append(sld_in[k])
358
359        if fabs(A[k]) < 1.0e-16:
360            # flat shell
361            z.append(z_current + thickness[k])
362            rho.append(sld_in[k])
363        else:
364            # exponential shell
365            # num_steps must be at least 1, so use floor()+1 rather than ceil
366            # to protect against a thickness0.
367            num_steps = np.floor(thickness[k]/dz) + 1
368            slope = (sld_out[k] - sld_in[k]) / expm1(A[k])
369            const = (sld_in[k] - slope)
370            for z_shell in np.linspace(0, thickness[k], num_steps+1):
371                z.append(z_current+z_shell)
372                rho.append(slope*exp(A[k]*z_shell/thickness[k]) + const)
373
374    # add in the solvent
375    z.append(z[-1])
376    rho.append(sld_solvent)
377    z.append(total_radius)
378    rho.append(sld_solvent)
379
380    return np.asarray(z), np.asarray(rho)
381
382# TODO: no random parameter function for onion model
383
384demo = {
385    "sld_solvent": 2.2,
386    "sld_core": 1.0,
387    "radius_core": 100,
388    "n_shells": 4,
389    "sld_in": [0.5, 1.5, 0.9, 2.0],
390    "sld_out": [nan, 0.9, 1.2, 1.6],
391    "thickness": [50, 75, 150, 75],
392    "A": [0, -1, 1e-4, 1],
393    # Could also specify them individually as
394    # "A1": 0, "A2": -1, "A3": 1e-4, "A4": 1,
395    #"radius_core_pd_n": 10,
396    #"radius_core_pd": 0.4,
397    #"thickness4_pd_n": 10,
398    #"thickness4_pd": 0.4,
399    }
Note: See TracBrowser for help on using the repository browser.