source: sasmodels/sasmodels/models/core_shell_parallelepiped.py @ 5bc6d21

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 5bc6d21 was 5bc6d21, checked in by butler, 2 years ago

more standardization of text

moved validation to its own section (is that the standard?) and moved
parallelipiped frong stuff around so start with Definition as per
standard.

  • 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
145
146Validation
147----------
148
149Cross-checked against hollow rectangular prism and rectangular prism for equal
150thickness overlapping sides, and by Monte Carlo sampling of points within the
151shape for non-uniform, non-overlapping sides.
152
153
154References
155----------
156
157.. [#] P Mittelbach and G Porod, *Acta Physica Austriaca*, 14 (1961) 185-211
158    Equations (1), (13-14). (in German)
159.. [#] D Singh (2009). *Small angle scattering studies of self assembly in
160   lipid mixtures*, Johns Hopkins University Thesis (2009) 223-225. `Available
161   from Proquest <http://search.proquest.com/docview/304915826?accountid
162   =26379>`_
163
164Authorship and Verification
165----------------------------
166
167* **Author:** NIST IGOR/DANSE **Date:** pre 2010
168* **Converted to sasmodels by:** Miguel Gonzalez **Date:** February 26, 2016
169* **Last Modified by:** Paul Kienzle **Date:** October 17, 2017
170"""
171
172import numpy as np
173from numpy import pi, inf, sqrt, cos, sin
174
175name = "core_shell_parallelepiped"
176title = "Rectangular solid with a core-shell structure."
177description = """
178     P(q)=
179"""
180category = "shape:parallelepiped"
181
182#             ["name", "units", default, [lower, upper], "type","description"],
183parameters = [["sld_core", "1e-6/Ang^2", 1, [-inf, inf], "sld",
184               "Parallelepiped core scattering length density"],
185              ["sld_a", "1e-6/Ang^2", 2, [-inf, inf], "sld",
186               "Parallelepiped A rim scattering length density"],
187              ["sld_b", "1e-6/Ang^2", 4, [-inf, inf], "sld",
188               "Parallelepiped B rim scattering length density"],
189              ["sld_c", "1e-6/Ang^2", 2, [-inf, inf], "sld",
190               "Parallelepiped C rim scattering length density"],
191              ["sld_solvent", "1e-6/Ang^2", 6, [-inf, inf], "sld",
192               "Solvent scattering length density"],
193              ["length_a", "Ang", 35, [0, inf], "volume",
194               "Shorter side of the parallelepiped"],
195              ["length_b", "Ang", 75, [0, inf], "volume",
196               "Second side of the parallelepiped"],
197              ["length_c", "Ang", 400, [0, inf], "volume",
198               "Larger side of the parallelepiped"],
199              ["thick_rim_a", "Ang", 10, [0, inf], "volume",
200               "Thickness of A rim"],
201              ["thick_rim_b", "Ang", 10, [0, inf], "volume",
202               "Thickness of B rim"],
203              ["thick_rim_c", "Ang", 10, [0, inf], "volume",
204               "Thickness of C rim"],
205              ["theta", "degrees", 0, [-360, 360], "orientation",
206               "c axis to beam angle"],
207              ["phi", "degrees", 0, [-360, 360], "orientation",
208               "rotation about beam"],
209              ["psi", "degrees", 0, [-360, 360], "orientation",
210               "rotation about c axis"],
211             ]
212
213source = ["lib/gauss76.c", "core_shell_parallelepiped.c"]
214
215
216def ER(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c):
217    """
218        Return equivalent radius (ER)
219    """
220    from .parallelepiped import ER as ER_p
221
222    a = length_a + 2*thick_rim_a
223    b = length_b + 2*thick_rim_b
224    c = length_c + 2*thick_rim_c
225    return ER_p(a, b, c)
226
227# VR defaults to 1.0
228
229def random():
230    outer = 10**np.random.uniform(1, 4.7, size=3)
231    thick = np.random.beta(0.5, 0.5, size=3)*(outer-2) + 1
232    length = outer - thick
233    pars = dict(
234        length_a=length[0],
235        length_b=length[1],
236        length_c=length[2],
237        thick_rim_a=thick[0],
238        thick_rim_b=thick[1],
239        thick_rim_c=thick[2],
240    )
241    return pars
242
243# parameters for demo
244demo = dict(scale=1, background=0.0,
245            sld_core=1, sld_a=2, sld_b=4, sld_c=2, sld_solvent=6,
246            length_a=35, length_b=75, length_c=400,
247            thick_rim_a=10, thick_rim_b=10, thick_rim_c=10,
248            theta=0, phi=0, psi=0,
249            length_a_pd=0.1, length_a_pd_n=1,
250            length_b_pd=0.1, length_b_pd_n=1,
251            length_c_pd=0.1, length_c_pd_n=1,
252            thick_rim_a_pd=0.1, thick_rim_a_pd_n=1,
253            thick_rim_b_pd=0.1, thick_rim_b_pd_n=1,
254            thick_rim_c_pd=0.1, thick_rim_c_pd_n=1,
255            theta_pd=10, theta_pd_n=1,
256            phi_pd=10, phi_pd_n=1,
257            psi_pd=10, psi_pd_n=1)
258
259# rkh 7/4/17 add random unit test for 2d, note make all params different,
260# 2d values not tested against other codes or models
261if 0:  # pak: model rewrite; need to update tests
262    qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.)
263    tests = [[{}, 0.2, 0.533149288477],
264             [{}, [0.2], [0.533149288477]],
265             [{'theta':10.0, 'phi':20.0}, (qx, qy), 0.0853299803222],
266             [{'theta':10.0, 'phi':20.0}, [(qx, qy)], [0.0853299803222]],
267            ]
268    del qx, qy  # not necessary to delete, but cleaner
Note: See TracBrowser for help on using the repository browser.