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

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since cf3d0ce was 99658f6, checked in by grethevj, 5 years ago

updated ER functions including cylinder excluded volume, to match 4.x

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