source: sasmodels/sasmodels/models/hollow_rectangular_prism.py @ 455aaa1

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 455aaa1 was 455aaa1, checked in by butler, 6 years ago

update hollow_rectangular_prism docs

Add statement about vol fraction being for shell only. Also fix typos,
clean up/normalize ref and authorship sections and fix/normalize
hollow_cylinder refs and authorship

  • Property mode set to 100644
File size: 8.0 KB
Line 
1# rectangular_prism model
2# Note: model title and parameter table are inserted automatically
3r"""
4Definition
5----------
6
7This model provides the form factor, $P(q)$, for a hollow rectangular
8parallelepiped with a wall of thickness $\Delta$. The 1D scattering intensity
9for this model is calculated by forming the difference of the amplitudes of two
10massive parallelepipeds differing in their outermost dimensions in each
11direction by the same length increment $2\Delta$ (\ [#Nayuk2012]_ Nayuk, 2012).
12
13As in the case of the massive parallelepiped model (:ref:`rectangular-prism`),
14the scattering amplitude is computed for a particular orientation of the
15parallelepiped with respect to the scattering vector and then averaged over all
16possible orientations, giving
17
18.. math::
19  P(q) =  \frac{1}{V^2} \frac{2}{\pi} \times \, \int_0^{\frac{\pi}{2}} \,
20  \int_0^{\frac{\pi}{2}} A_{P\Delta}^2(q) \, \sin\theta \, d\theta \, d\phi
21
22where $\theta$ is the angle between the $z$ axis and the longest axis
23of the parallelepiped, $\phi$ is the angle between the scattering vector
24(lying in the $xy$ plane) and the $y$ axis, and
25
26.. math::
27  :nowrap:
28
29  \begin{align*}
30  A_{P\Delta}(q) & =  A B C
31    \left[\frac{\sin \bigl( q \frac{C}{2} \cos\theta \bigr)}
32    {\left( q \frac{C}{2} \cos\theta \right)} \right]
33    \left[\frac{\sin \bigl( q \frac{A}{2} \sin\theta \sin\phi \bigr)}
34    {\left( q \frac{A}{2} \sin\theta \sin\phi \right)}\right]
35    \left[\frac{\sin \bigl( q \frac{B}{2} \sin\theta \cos\phi \bigr)}
36    {\left( q \frac{B}{2} \sin\theta \cos\phi \right)}\right] \\
37    & - 8
38    \left(\frac{A}{2}-\Delta\right) \left(\frac{B}{2}-\Delta\right) \left(\frac{C}{2}-\Delta\right)
39    \left[ \frac{\sin \bigl[ q \bigl(\frac{C}{2}-\Delta\bigr) \cos\theta \bigr]}
40    {q \bigl(\frac{C}{2}-\Delta\bigr) \cos\theta} \right]
41    \left[ \frac{\sin \bigl[ q \bigl(\frac{A}{2}-\Delta\bigr) \sin\theta \sin\phi \bigr]}
42    {q \bigl(\frac{A}{2}-\Delta\bigr) \sin\theta \sin\phi} \right]
43    \left[ \frac{\sin \bigl[ q \bigl(\frac{B}{2}-\Delta\bigr) \sin\theta \cos\phi \bigr]}
44    {q \bigl(\frac{B}{2}-\Delta\bigr) \sin\theta \cos\phi} \right]
45  \end{align*}
46
47where $A$, $B$ and $C$ are the external sides of the parallelepiped fulfilling
48$A \le B \le C$, and the volume $V$ of the parallelepiped is
49
50.. math::
51  V = A B C \, - \, (A - 2\Delta) (B - 2\Delta) (C - 2\Delta)
52
53The 1D scattering intensity is then calculated as
54
55.. math::
56  I(q) = \text{scale} \times V \times (\rho_\text{p} -
57  \rho_\text{solvent})^2 \times P(q) + \text{background}
58
59where $\rho_\text{p}$ is the scattering length density of the parallelepiped,
60$\rho_\text{solvent}$ is the scattering length density of the solvent,
61and (if the data are in absolute units) *scale* represents the volume fraction
62(which is unitless) of the rectangular shell of material (i.e. not including
63the volume of the solvent filled core).
64
65For 2d data the orientation of the particle is required, described using
66angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details
67of the calculation and angular dispersions see :ref:`orientation` .
68The angle $\Psi$ is the rotational angle around the long *C* axis. For example,
69$\Psi = 0$ when the *B* axis is parallel to the *x*-axis of the detector.
70
71For 2d, constraints must be applied during fitting to ensure that the inequality
72$A < B < C$ is not violated, and hence the correct definition of angles is
73preserved. The calculation will not report an error if the inequality is *not*
74preserved, but the results may be not correct.
75
76.. figure:: img/parallelepiped_angle_definition.png
77
78    Definition of the angles for oriented hollow rectangular prism.
79    Note that rotation $\theta$, initially in the $xz$ plane, is carried out first, then
80    rotation $\phi$ about the $z$ axis, finally rotation $\Psi$ is now around the axis of the prism.
81    The neutron or X-ray beam is along the $z$ axis.
82
83.. figure:: img/parallelepiped_angle_projection.png
84
85    Examples of the angles for oriented hollow rectangular prisms against the
86    detector plane.
87
88
89Validation
90----------
91
92Validation of the code was conducted by qualitatively comparing the output
93of the 1D model to the curves shown in (Nayuk, 2012).
94
95
96References
97----------
98
99.. [#Nayuk2012] R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854
100
101
102Authorship and Verification
103----------------------------
104
105* **Author:** Miguel Gonzales **Date:** February 26, 2016
106* **Last Modified by:** Paul Kienzle **Date:** December 14, 2017
107* **Last Reviewed by:** Paul Butler **Date:** September 06, 2018
108"""
109
110import numpy as np
111from numpy import pi, inf, sqrt
112
113name = "hollow_rectangular_prism"
114title = "Hollow rectangular parallelepiped with uniform scattering length density."
115description = """
116    I(q)= scale*V*(sld - sld_solvent)^2*P(q,theta,phi)+background
117        P(q,theta,phi) = (2/pi/V^2) * double integral from 0 to pi/2 of ...
118           (AP1-AP2)^2(q)*sin(theta)*dtheta*dphi
119        AP1 = S(q*C*cos(theta)/2) * S(q*A*sin(theta)*sin(phi)/2) * S(q*B*sin(theta)*cos(phi)/2)
120        AP2 = S(q*C'*cos(theta)) * S(q*A'*sin(theta)*sin(phi)) * S(q*B'*sin(theta)*cos(phi))
121        C' = (C/2-thickness)
122        B' = (B/2-thickness)
123        A' = (A/2-thickness)
124        S(x) = sin(x)/x
125"""
126category = "shape:parallelepiped"
127
128#             ["name", "units", default, [lower, upper], "type","description"],
129parameters = [["sld", "1e-6/Ang^2", 6.3, [-inf, inf], "sld",
130               "Parallelepiped scattering length density"],
131              ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "sld",
132               "Solvent scattering length density"],
133              ["length_a", "Ang", 35, [0, inf], "volume",
134               "Shorter side of the parallelepiped"],
135              ["b2a_ratio", "Ang", 1, [0, inf], "volume",
136               "Ratio sides b/a"],
137              ["c2a_ratio", "Ang", 1, [0, inf], "volume",
138               "Ratio sides c/a"],
139              ["thickness", "Ang", 1, [0, inf], "volume",
140               "Thickness of parallelepiped"],
141              ["theta", "degrees", 0, [-360, 360], "orientation",
142               "c axis to beam angle"],
143              ["phi", "degrees", 0, [-360, 360], "orientation",
144               "rotation about beam"],
145              ["psi", "degrees", 0, [-360, 360], "orientation",
146               "rotation about c axis"],
147             ]
148
149source = ["lib/gauss76.c", "hollow_rectangular_prism.c"]
150
151def ER(length_a, b2a_ratio, c2a_ratio, thickness):
152    """
153    Return equivalent radius (ER)
154    thickness parameter not used
155    """
156    b_side = length_a * b2a_ratio
157    c_side = length_a * c2a_ratio
158
159    # surface average radius (rough approximation)
160    surf_rad = sqrt(length_a * b_side / pi)
161
162    ddd = 0.75 * surf_rad * (2 * surf_rad * c_side + (c_side + surf_rad) * (c_side + pi * surf_rad))
163    return 0.5 * (ddd) ** (1. / 3.)
164
165def VR(length_a, b2a_ratio, c2a_ratio, thickness):
166    """
167    Return shell volume and total volume
168    """
169    b_side = length_a * b2a_ratio
170    c_side = length_a * c2a_ratio
171    a_core = length_a - 2.0*thickness
172    b_core = b_side - 2.0*thickness
173    c_core = c_side - 2.0*thickness
174    vol_core = a_core * b_core * c_core
175    vol_total = length_a * b_side * c_side
176    vol_shell = vol_total - vol_core
177    return vol_total, vol_shell
178
179
180def random():
181    a, b, c = 10**np.random.uniform(1, 4.7, size=3)
182    # Thickness is limited to 1/2 the smallest dimension
183    # Use a distribution with a preference for thin shell or thin core
184    # Avoid core,shell radii < 1
185    min_dim = 0.5*min(a, b, c)
186    thickness = np.random.beta(0.5, 0.5)*(min_dim-2) + 1
187    #print(a, b, c, thickness, thickness/min_dim)
188    pars = dict(
189        length_a=a,
190        b2a_ratio=b/a,
191        c2a_ratio=c/a,
192        thickness=thickness,
193    )
194    return pars
195
196
197# parameters for demo
198demo = dict(scale=1, background=0,
199            sld=6.3, sld_solvent=1.0,
200            length_a=35, b2a_ratio=1, c2a_ratio=1, thickness=1,
201            length_a_pd=0.1, length_a_pd_n=10,
202            b2a_ratio_pd=0.1, b2a_ratio_pd_n=1,
203            c2a_ratio_pd=0.1, c2a_ratio_pd_n=1)
204
205tests = [[{}, 0.2, 0.76687283098],
206         [{}, [0.2], [0.76687283098]],
207        ]
Note: See TracBrowser for help on using the repository browser.