source: sasmodels/sasmodels/models/core_shell_parallelepiped.py @ 96153e4

ticket-608-user-defined-weightsticket_1156
Last change on this file since 96153e4 was 96153e4, checked in by butler, 10 months ago

more corrections and normalizations

Addreses #896. Brings parallelepiped and core-shell parallelepiped more
into line with each other and corrects angle definitions refering to
orietnational distribution docs for details. Still need to sort out
with Paul Kienzle the proper order the angles are applied.

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