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@…>, 5 years ago

merge with master

  • Property mode set to 100644
File size: 12.0 KB
RevLine 
[fdb1487]1r"""
2This model provides the form factor, $P(q)$, for a multi-shell sphere where
[bccb40f]3the scattering length density (SLD) of each shell is described by an
[fdb1487]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
[db1d9d5]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$.
[fdb1487]12
13Definition
14----------
15
16The 1D scattering intensity is calculated in the following way
17
18.. math::
19
[63c6a08]20    P(q) = [f]^2 / V_\text{particle}
[fdb1487]21
22where
23
24.. math::
[63c6a08]25    :nowrap:
[fdb1487]26
[63c6a08]27    \begin{align*}
28    f &= f_\text{core}
[fdb1487]29            + \left(\sum_{\text{shell}=1}^N f_\text{shell}\right)
30            + f_\text{solvent}
[63c6a08]31    \end{align*}
[fdb1487]32
33The shells are spherically symmetric with particle density $\rho(r)$ and
34constant SLD within the core and solvent, so
35
36.. math::
[63c6a08]37    :nowrap:
38
39    \begin{align*}
[fdb1487]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}
[63c6a08]52    \end{align*}
[fdb1487]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
[db1d9d5]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)$.
[fdb1487]64
[db1d9d5]65Now consider the SLD of a shell defined by
[fdb1487]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
[db1d9d5]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$):
[fdb1487]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
[ca04add]92        \right] \frac{\sin(qr)}{qr}\,r^2\,\mathrm{d}r \\
[fdb1487]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
[db1d9d5]98where
[fdb1487]99
100.. math::
101    :nowrap:
102
103    \begin{align*}
104    B&=\frac{\rho_\text{out} - \rho_\text{in}}{e^A-1}
[ca04add]105         & C &= \frac{\rho_\text{in}e^A - \rho_\text{out}}{e^A-1} \\
[db1d9d5]106
[fdb1487]107    \alpha_\text{in} &= A\frac{r_{\text{shell}-1}}{\Delta t_\text{shell}}
[ca04add]108         & \alpha_\text{out} &= A\frac{r_\text{shell}}{\Delta t_\text{shell}} \\
[db1d9d5]109
[fdb1487]110    \beta_\text{in} &= qr_{\text{shell}-1}
[ca04add]111        & \beta_\text{out} &= qr_\text{shell} \\
[fdb1487]112    \end{align*}
113
[db1d9d5]114and
[fdb1487]115
116 .. math::
117
[db1d9d5]118     h(x,y) = \frac{x \sin(y) - y\cos(y)}{(x^2+y^2)y}
[fdb1487]119               - \frac{(x^2-y^2)\sin(y) - 2xy\cos(y)}{(x^2+y^2)^2y}
120
121
[db1d9d5]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
[fdb1487]131
132.. math::
[63c6a08]133    :nowrap:
[fdb1487]134
[63c6a08]135    \begin{align*}
[fdb1487]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}}
[63c6a08]156    \end{align*}
[fdb1487]157
[db1d9d5]158
159**Constant SLD** ($A = 0$):
160
161When $A = 0$ the exponential function has no dependence on the radius (meaning
[bccb40f]162$\rho_\text{out}$ is ignored in this case) and becomes flat. We set the constant
[fdb1487]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
[0507e09]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
192`onion.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/onion.py>`_
193
194`onion.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/onion.c>`_
195
196Authorship and Verification
197----------------------------
198
[a34b811]199* **Author:**
200* **Last Modified by:**
[db1d9d5]201* **Last Reviewed by:** Steve King **Date:** March 28, 2019
[0507e09]202* **Source added by :** Steve King **Date:** March 25, 2019
[fdb1487]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
[2d81cfe]292from math import fabs, exp, expm1
293
[fdb1487]294import numpy as np
295from numpy import inf, nan
296
297name = "onion"
298title = "Onion shell model with constant, linear or exponential density"
299
300description = """\
[db1d9d5]301Form factor of multishells normalized by the volume. Here each shell is
[fdb1487]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
[db1d9d5]314        Note that the 'sld_out' becomes NULL in this case.
[fdb1487]315
316        background:background,
317        rad_core0: radius of sphere(core)
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
[b0696e1]327# TODO: n is a volume parameter that is not polydisperse
[fdb1487]328
[db1d9d5]329# NOTE: Joachim Wuttke has suggested an alternative parameterisation
330#       in Ticket #1107
331
[40a87fa]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"],
[db1d9d5]338    ["n_shells", "", 1, [0, 10], "volume", "number of shells (must be integer)"],
[40a87fa]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
[fdb1487]345
[925ad6e]346source = ["lib/sas_3j1x_x.c", "onion.c"]
[a0494e9]347single = False
[71b751d]348have_Fq = True
[a34b811]349radius_effective_modes = ["outer radius"]
[6a8fdfe]350
[ce896fd]351profile_axes = ['Radius (A)', 'SLD (1e-6/A^2)']
[e187b25]352def profile(sld_core, radius_core, sld_solvent, n_shells,
353            sld_in, sld_out, thickness, A):
[fdb1487]354    """
[fa5fd8d]355    Returns shape profile with x=radius, y=SLD.
[fdb1487]356    """
[768c0c4]357    n_shells = int(n_shells+0.5)
[e187b25]358    total_radius = 1.25*(sum(thickness[:n_shells]) + radius_core + 1)
[40a87fa]359    dz = total_radius/400  # 400 points for a smooth plot
[fdb1487]360
[40a87fa]361    z = []
[e187b25]362    rho = []
[fdb1487]363
364    # add in the core
[40a87fa]365    z.append(0)
[e187b25]366    rho.append(sld_core)
[40a87fa]367    z.append(radius_core)
[e187b25]368    rho.append(sld_core)
[fdb1487]369
370    # add in the shells
[3cd1001]371    for k in range(int(n_shells)):
[fdb1487]372        # Left side of each shells
[40a87fa]373        z_current = z[-1]
374        z.append(z_current)
[e187b25]375        rho.append(sld_in[k])
[fdb1487]376
377        if fabs(A[k]) < 1.0e-16:
378            # flat shell
[40a87fa]379            z.append(z_current + thickness[k])
[bccb40f]380            rho.append(sld_in[k])
[fdb1487]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.
[40a87fa]385            num_steps = np.floor(thickness[k]/dz) + 1
[e187b25]386            slope = (sld_out[k] - sld_in[k]) / expm1(A[k])
387            const = (sld_in[k] - slope)
[40a87fa]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)
[2d81cfe]391
[fdb1487]392    # add in the solvent
[40a87fa]393    z.append(z[-1])
[e187b25]394    rho.append(sld_solvent)
[40a87fa]395    z.append(total_radius)
[e187b25]396    rho.append(sld_solvent)
[fdb1487]397
[40a87fa]398    return np.asarray(z), np.asarray(rho)
[fdb1487]399
[b297ba9]400# TODO: no random parameter function for onion model
401
[fdb1487]402demo = {
[ce896fd]403    "sld_solvent": 2.2,
404    "sld_core": 1.0,
[9762341]405    "radius_core": 100,
[a0494e9]406    "n_shells": 4,
[ce896fd]407    "sld_in": [0.5, 1.5, 0.9, 2.0],
408    "sld_out": [nan, 0.9, 1.2, 1.6],
[fdb1487]409    "thickness": [50, 75, 150, 75],
410    "A": [0, -1, 1e-4, 1],
411    # Could also specify them individually as
[a0494e9]412    # "A1": 0, "A2": -1, "A3": 1e-4, "A4": 1,
[9762341]413    #"radius_core_pd_n": 10,
414    #"radius_core_pd": 0.4,
[d119f34]415    #"thickness4_pd_n": 10,
416    #"thickness4_pd": 0.4,
[fdb1487]417    }
Note: See TracBrowser for help on using the repository browser.