source: sasmodels/sasmodels/models/onion.py @ 63c6a08

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 63c6a08 was 63c6a08, checked in by Paul Kienzle <pkienzle@…>, 8 years ago

fix some latex problems

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