source: sasmodels/sasmodels/models/hollow_rectangular_prism.py @ a34b811

ticket-1257-vesicle-productticket_1156ticket_822_more_unit_tests
Last change on this file since a34b811 was a34b811, checked in by Paul Kienzle <pkienzle@…>, 5 years ago

use radius_effective/radius_effective_mode/radius_effective_modes consistently throughout the code

  • Property mode set to 100644
File size: 7.8 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.. [#] L. Onsager, *Ann. New York Acad. Sci.*, 51 (1949) 627-659
101
102Source
103------
104
105`hollow_rectangular_prism.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/hollow_rectangular_prism.py>`_
106
107`hollow_rectangular_prism.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/hollow_rectangular_prism.c>`_
108
109Authorship and Verification
110----------------------------
111
112* **Author:** Miguel Gonzales **Date:** February 26, 2016
113* **Last Modified by:** Paul Kienzle **Date:** December 14, 2017
114* **Last Reviewed by:** Paul Butler **Date:** September 06, 2018
115* **Source added by :** Steve King **Date:** March 25, 2019
116"""
117
118import numpy as np
119from numpy import inf
120
121name = "hollow_rectangular_prism"
122title = "Hollow rectangular parallelepiped with uniform scattering length density."
123description = """
124    I(q)= scale*V*(sld - sld_solvent)^2*P(q,theta,phi)+background
125        P(q,theta,phi) = (2/pi/V^2) * double integral from 0 to pi/2 of ...
126           (AP1-AP2)^2(q)*sin(theta)*dtheta*dphi
127        AP1 = S(q*C*cos(theta)/2) * S(q*A*sin(theta)*sin(phi)/2) * S(q*B*sin(theta)*cos(phi)/2)
128        AP2 = S(q*C'*cos(theta)) * S(q*A'*sin(theta)*sin(phi)) * S(q*B'*sin(theta)*cos(phi))
129        C' = (C/2-thickness)
130        B' = (B/2-thickness)
131        A' = (A/2-thickness)
132        S(x) = sin(x)/x
133"""
134category = "shape:parallelepiped"
135
136#             ["name", "units", default, [lower, upper], "type","description"],
137parameters = [["sld", "1e-6/Ang^2", 6.3, [-inf, inf], "sld",
138               "Parallelepiped scattering length density"],
139              ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "sld",
140               "Solvent scattering length density"],
141              ["length_a", "Ang", 35, [0, inf], "volume",
142               "Shortest, external, size of the parallelepiped"],
143              ["b2a_ratio", "Ang", 1, [0, inf], "volume",
144               "Ratio sides b/a"],
145              ["c2a_ratio", "Ang", 1, [0, inf], "volume",
146               "Ratio sides c/a"],
147              ["thickness", "Ang", 1, [0, inf], "volume",
148               "Thickness of parallelepiped"],
149              ["theta", "degrees", 0, [-360, 360], "orientation",
150               "c axis to beam angle"],
151              ["phi", "degrees", 0, [-360, 360], "orientation",
152               "rotation about beam"],
153              ["psi", "degrees", 0, [-360, 360], "orientation",
154               "rotation about c axis"],
155             ]
156
157source = ["lib/gauss76.c", "hollow_rectangular_prism.c"]
158have_Fq = True
159radius_effective_modes = [
160    "equivalent cylinder excluded volume", "equivalent outer volume sphere",
161    "half length_a", "half length_b", "half length_c",
162    "equivalent outer circular cross-section",
163    "half ab diagonal", "half diagonal",
164    ]
165
166def random():
167    """Return a random parameter set for the model."""
168    a, b, c = 10**np.random.uniform(1, 4.7, size=3)
169    # Thickness is limited to 1/2 the smallest dimension
170    # Use a distribution with a preference for thin shell or thin core
171    # Avoid core,shell radii < 1
172    min_dim = 0.5*min(a, b, c)
173    thickness = np.random.beta(0.5, 0.5)*(min_dim-2) + 1
174    #print(a, b, c, thickness, thickness/min_dim)
175    pars = dict(
176        length_a=a,
177        b2a_ratio=b/a,
178        c2a_ratio=c/a,
179        thickness=thickness,
180    )
181    return pars
182
183
184# parameters for demo
185demo = dict(scale=1, background=0,
186            sld=6.3, sld_solvent=1.0,
187            length_a=35, b2a_ratio=1, c2a_ratio=1, thickness=1,
188            length_a_pd=0.1, length_a_pd_n=10,
189            b2a_ratio_pd=0.1, b2a_ratio_pd_n=1,
190            c2a_ratio_pd=0.1, c2a_ratio_pd_n=1)
191
192tests = [[{}, 0.2, 0.76687283098],
193         [{}, [0.2], [0.76687283098]],
194        ]
Note: See TracBrowser for help on using the repository browser.