source: sasmodels/sasmodels/models/core_shell_parallelepiped.py @ 331870d

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 331870d was 331870d, checked in by gonzalezm, 5 years ago

A few typos corrected

  • Property mode set to 100644
File size: 9.9 KB
Line 
1r"""
2Definition
3----------
4
5Calculates the form factor for a rectangular solid with a core-shell structure.
6The thickness and the scattering length density of the shell or
7"rim" can be different on each (pair) of faces.
8
9The form factor is normalized by the particle volume $V$ such that
10
11.. math::
12
13    I(q) = \frac{\text{scale}}{V} \langle P(q,\alpha,\beta) \rangle
14    + \text{background}
15
16where $\langle \ldots \rangle$ is an average over all possible orientations
17of the rectangular solid, and the usual $\Delta \rho^2 \ V^2$ term cannot be
18pulled out of the form factor term due to the multiple slds in the model.
19
20The core of the solid is defined by the dimensions $A$, $B$, $C$ such that
21$A < B < C$.
22
23.. figure:: img/parallelepiped_geometry.jpg
24
25   Core of the core shell parallelepiped with the corresponding definition
26   of sides.
27
28
29There are rectangular "slabs" of thickness $t_A$ that add to the $A$ dimension
30(on the $BC$ faces). There are similar slabs on the $AC$ $(=t_B)$ and $AB$
31$(=t_C)$ faces. The projection in the $AB$ plane is
32
33.. figure:: img/core_shell_parallelepiped_projection.jpg
34
35   AB cut through the core-shell parallelipiped showing the cross secion of
36   four of the six shell slabs. As can be seen, this model leaves **"gaps"**
37   at the corners of the solid.
38
39
40The total volume of the solid is thus given as
41
42.. math::
43
44    V = ABC + 2t_ABC + 2t_BAC + 2t_CAB
45
46The intensity calculated follows the :ref:`parallelepiped` model, with the
47core-shell intensity being calculated as the square of the sum of the
48amplitudes of the core and the slabs on the edges. The scattering amplitude is
49computed for a particular orientation of the core-shell parallelepiped with
50respect to the scattering vector and then averaged over all possible
51orientations, where $\alpha$ is the angle between the $z$ axis and the $C$ axis
52of the parallelepiped, and $\beta$ is the angle between the projection of the
53particle in the $xy$ detector plane and the $y$ axis.
54
55.. math::
56
57    P(q)=\frac {\int_{0}^{\pi/2}\int_{0}^{\pi/2}F^2(q,\alpha,\beta) \ sin\alpha
58    \ d\alpha \ d\beta} {\int_{0}^{\pi/2} \ sin\alpha \ d\alpha \ d\beta}
59
60and
61
62.. math::
63
64    F(q,\alpha,\beta)
65    &= (\rho_\text{core}-\rho_\text{solvent})
66       S(Q_A, A) S(Q_B, B) S(Q_C, C) \\
67    &+ (\rho_\text{A}-\rho_\text{solvent})
68        \left[S(Q_A, A+2t_A) - S(Q_A, A)\right] S(Q_B, B) S(Q_C, C) \\
69    &+ (\rho_\text{B}-\rho_\text{solvent})
70        S(Q_A, A) \left[S(Q_B, B+2t_B) - S(Q_B, B)\right] S(Q_C, C) \\
71    &+ (\rho_\text{C}-\rho_\text{solvent})
72        S(Q_A, A) S(Q_B, B) \left[S(Q_C, C+2t_C) - S(Q_C, C)\right]
73
74with
75
76.. math::
77
78    S(Q_X, L) = L \frac{\sin (\tfrac{1}{2} Q_X L)}{\tfrac{1}{2} Q_X L}
79
80and
81
82.. math::
83
84    Q_A &= q \sin\alpha \sin\beta \\
85    Q_B &= q \sin\alpha \cos\beta \\
86    Q_C &= q \cos\alpha
87
88
89where $\rho_\text{core}$, $\rho_\text{A}$, $\rho_\text{B}$ and $\rho_\text{C}$
90are the scattering lengths of the parallelepiped core, and the rectangular
91slabs of thickness $t_A$, $t_B$ and $t_C$, respectively. $\rho_\text{solvent}$
92is the scattering length of the solvent.
93
94.. note::
95
96   the code actually implements two substitutions: $d(cos\alpha)$ is
97   substituted for -$sin\alpha \ d\alpha$ (note that in the
98   :ref:`parallelepiped` code this is explicitly implemented with
99   $\sigma = cos\alpha$), and $\beta$ is set to $\beta = u \pi/2$ so that
100   $du = \pi/2 \ d\beta$.  Thus both integrals go from 0 to 1 rather than 0
101   to $\pi/2$.
102
103FITTING NOTES
104~~~~~~~~~~~~~
105
106If the scale is set equal to the particle volume fraction, $\phi$, the returned
107value is the scattered intensity per unit volume, $I(q) = \phi P(q)$. However,
108**no interparticle interference effects are included in this calculation.**
109
110There are many parameters in this model. Hold as many fixed as possible with
111known values, or you will certainly end up at a solution that is unphysical.
112
113The returned value is in units of |cm^-1|, on absolute scale.
114
115NB: The 2nd virial coefficient of the core_shell_parallelepiped is calculated
116based on the the averaged effective radius $(=\sqrt{(A+2t_A)(B+2t_B)/\pi})$
117and length $(C+2t_C)$ values, after appropriately sorting the three dimensions
118to give an oblate or prolate particle, to give an effective radius
119for $S(q)$ when $P(q) * S(q)$ is applied.
120
121For 2d data the orientation of the particle is required, described using
122angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below. For further
123details of the calculation and angular dispersions see :ref:`orientation`.
124The angle $\Psi$ is the rotational angle around the *long_c* axis. For example,
125$\Psi = 0$ when the *short_b* axis is parallel to the *x*-axis of the detector.
126
127.. note:: For 2d, constraints must be applied during fitting to ensure that the
128   inequality $A < B < C$ is not violated, and hence the correct definition
129   of angles is preserved. The calculation will not report an error,
130   but the results may be not correct.
131
132.. figure:: img/parallelepiped_angle_definition.png
133
134    Definition of the angles for oriented core-shell parallelepipeds.
135    Note that rotation $\theta$, initially in the $xz$ plane, is carried
136    out first, then rotation $\phi$ about the $z$ axis, finally rotation
137    $\Psi$ is now around the axis of the particle. The neutron or X-ray
138    beam is along the $z$ axis.
139
140.. figure:: img/parallelepiped_angle_projection.png
141
142    Examples of the angles for oriented core-shell parallelepipeds against the
143    detector plane.
144
145References
146----------
147
148.. [#] P Mittelbach and G Porod, *Acta Physica Austriaca*, 14 (1961) 185-211
149    Equations (1), (13-14). (in German)
150.. [#] D Singh (2009). *Small angle scattering studies of self assembly in
151   lipid mixtures*, Johns Hopkins University Thesis (2009) 223-225. `Available
152   from Proquest <http://search.proquest.com/docview/304915826?accountid
153   =26379>`_
154
155Authorship and Verification
156----------------------------
157
158* **Author:** NIST IGOR/DANSE **Date:** pre 2010
159* **Converted to sasmodels by:** Miguel Gonzalez **Date:** February 26, 2016
160* **Last Modified by:** Paul Kienzle **Date:** October 17, 2017
161* Cross-checked against hollow rectangular prism and rectangular prism for
162  equal thickness overlapping sides, and by Monte Carlo sampling of points
163  within the shape for non-uniform, non-overlapping sides.
164"""
165
166import numpy as np
167from numpy import pi, inf, sqrt, cos, sin
168
169name = "core_shell_parallelepiped"
170title = "Rectangular solid with a core-shell structure."
171description = """
172     P(q)=
173"""
174category = "shape:parallelepiped"
175
176#             ["name", "units", default, [lower, upper], "type","description"],
177parameters = [["sld_core", "1e-6/Ang^2", 1, [-inf, inf], "sld",
178               "Parallelepiped core scattering length density"],
179              ["sld_a", "1e-6/Ang^2", 2, [-inf, inf], "sld",
180               "Parallelepiped A rim scattering length density"],
181              ["sld_b", "1e-6/Ang^2", 4, [-inf, inf], "sld",
182               "Parallelepiped B rim scattering length density"],
183              ["sld_c", "1e-6/Ang^2", 2, [-inf, inf], "sld",
184               "Parallelepiped C rim scattering length density"],
185              ["sld_solvent", "1e-6/Ang^2", 6, [-inf, inf], "sld",
186               "Solvent scattering length density"],
187              ["length_a", "Ang", 35, [0, inf], "volume",
188               "Shorter side of the parallelepiped"],
189              ["length_b", "Ang", 75, [0, inf], "volume",
190               "Second side of the parallelepiped"],
191              ["length_c", "Ang", 400, [0, inf], "volume",
192               "Larger side of the parallelepiped"],
193              ["thick_rim_a", "Ang", 10, [0, inf], "volume",
194               "Thickness of A rim"],
195              ["thick_rim_b", "Ang", 10, [0, inf], "volume",
196               "Thickness of B rim"],
197              ["thick_rim_c", "Ang", 10, [0, inf], "volume",
198               "Thickness of C rim"],
199              ["theta", "degrees", 0, [-360, 360], "orientation",
200               "c axis to beam angle"],
201              ["phi", "degrees", 0, [-360, 360], "orientation",
202               "rotation about beam"],
203              ["psi", "degrees", 0, [-360, 360], "orientation",
204               "rotation about c axis"],
205             ]
206
207source = ["lib/gauss76.c", "core_shell_parallelepiped.c"]
208
209
210def ER(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c):
211    """
212        Return equivalent radius (ER)
213    """
214    from .parallelepiped import ER as ER_p
215
216    a = length_a + 2*thick_rim_a
217    b = length_b + 2*thick_rim_b
218    c = length_c + 2*thick_rim_c
219    return ER_p(a, b, c)
220
221# VR defaults to 1.0
222
223def random():
224    outer = 10**np.random.uniform(1, 4.7, size=3)
225    thick = np.random.beta(0.5, 0.5, size=3)*(outer-2) + 1
226    length = outer - thick
227    pars = dict(
228        length_a=length[0],
229        length_b=length[1],
230        length_c=length[2],
231        thick_rim_a=thick[0],
232        thick_rim_b=thick[1],
233        thick_rim_c=thick[2],
234    )
235    return pars
236
237# parameters for demo
238demo = dict(scale=1, background=0.0,
239            sld_core=1, sld_a=2, sld_b=4, sld_c=2, sld_solvent=6,
240            length_a=35, length_b=75, length_c=400,
241            thick_rim_a=10, thick_rim_b=10, thick_rim_c=10,
242            theta=0, phi=0, psi=0,
243            length_a_pd=0.1, length_a_pd_n=1,
244            length_b_pd=0.1, length_b_pd_n=1,
245            length_c_pd=0.1, length_c_pd_n=1,
246            thick_rim_a_pd=0.1, thick_rim_a_pd_n=1,
247            thick_rim_b_pd=0.1, thick_rim_b_pd_n=1,
248            thick_rim_c_pd=0.1, thick_rim_c_pd_n=1,
249            theta_pd=10, theta_pd_n=1,
250            phi_pd=10, phi_pd_n=1,
251            psi_pd=10, psi_pd_n=1)
252
253# rkh 7/4/17 add random unit test for 2d, note make all params different,
254# 2d values not tested against other codes or models
255if 0:  # pak: model rewrite; need to update tests
256    qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.)
257    tests = [[{}, 0.2, 0.533149288477],
258             [{}, [0.2], [0.533149288477]],
259             [{'theta':10.0, 'phi':20.0}, (qx, qy), 0.0853299803222],
260             [{'theta':10.0, 'phi':20.0}, [(qx, qy)], [0.0853299803222]],
261            ]
262    del qx, qy  # not necessary to delete, but cleaner
Note: See TracBrowser for help on using the repository browser.