source: sasmodels/sasmodels/models/onion.py

Last change on this file was d57b06c, checked in by Paul Kienzle <pkienzle@…>, 15 months ago

Merge remote-tracking branch 'origin/master' into ticket-1263-source-link

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