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

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

First edits to core shell parallelipiped docs

Part f ticket 896 effort to ensure the documentation matches the code
now that the code has been fixed and validated.

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