source: sasmodels/sasmodels/models/hollow_rectangular_prism.py @ 0e55afe

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 0e55afe was 0e55afe, checked in by Paul Kienzle <pkienzle@…>, 6 years ago

Merge branch 'master' into ticket-786

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