source: sasmodels/sasmodels/models/onion.py @ 46ed760

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

enable multiplicity models, though they are not yet correct

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