# source:sasmodels/sasmodels/models/parallelepiped.py@dbf1a60

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

parallelepiped noting the change of integration varialbes in the
computation. Cleaned up and final corrections to the core shell
documentation and did some cleaning of parallelipiped. In particular
tried to bring a bit more consistency between the docs.

• Property mode set to 100644
File size: 10.1 KB
Line
1# parallelepiped model
2# Note: model title and parameter table are inserted automatically
3r"""
4The form factor is normalized by the particle volume.
5For information about polarised and magnetic scattering, see
6the :ref:magnetism documentation.
7
8Definition
9----------
10
11 This model calculates the scattering from a rectangular parallelepiped
12 (:numref:parallelepiped-image).
13 If you need to apply polydispersity, see also :ref:rectangular-prism.
14
15.. _parallelepiped-image:
16
17
18.. figure:: img/parallelepiped_geometry.jpg
19
20   Parallelepiped with the corresponding definition of sides.
21
22The three dimensions of the parallelepiped (strictly here a cuboid) may be
23given in *any* size order. To avoid multiple fit solutions, especially
24with Monte-Carlo fit methods, it may be advisable to restrict their ranges.
25There may be a number of closely similar "best fits", so some trial and
26error, or fixing of some dimensions at expected values, may help.
27
28The 1D scattering intensity $I(q)$ is calculated as:
29
30.. Comment by Miguel Gonzalez:
31   I am modifying the original text because I find the notation a little bit
32   confusing. I think that in most textbooks/papers, the notation P(Q) is
33   used for the form factor (adim, P(Q=0)=1), although F(q) seems also to
34   be used. But here (as for many other models), P(q) is used to represent
35   the scattering intensity (in cm-1 normally). It would be good to agree on
36   a common notation.
37
38.. math::
39
40    I(q) = \frac{\text{scale}}{V} (\Delta\rho \cdot V)^2
41           \left< P(q, \alpha, \beta) \right> + \text{background}
42
43where the volume $V = A B C$, the contrast is defined as
44$\Delta\rho = \rho_\text{p} - \rho_\text{solvent}$, $P(q, \alpha, \beta)$
45is the form factor corresponding to a parallelepiped oriented
46at an angle $\alpha$ (angle between the long axis C and $\vec q$), and $\beta$
47( the angle between the projection of the particle in the $xy$ detector plane
48and the $y$ axis) and the averaging $\left<\ldots\right>$ is applied over all
49orientations.
50
51Assuming $a = A/B < 1$, $b = B /B = 1$, and $c = C/B > 1$, the
52form factor is given by (Mittelbach and Porod, 1961 [#Mittelbach]_)
53
54.. math::
55
56    P(q, \alpha) = \int_0^1 \phi_Q\left(\mu \sqrt{1-\sigma^2},a\right)
57        \left[S(\mu c \sigma/2)\right]^2 d\sigma
58
59with
60
61.. math::
62
63    \phi_Q(\mu,a) &= \int_0^1
64        \left\{S\left[\frac{\mu}{2}\cos\left(\frac{\pi}{2}u\right)\right]
65               S\left[\frac{\mu a}{2}\sin\left(\frac{\pi}{2}u\right)\right]
66               \right\}^2 du \\
67    S(x) &= \frac{\sin x}{x} \\
68    \mu &= qB
69
70where substitution of $\sigma = cos\alpha$ and $\beta = \pi/2 \ u$ have been
71applied.
72
73NB: The 2nd virial coefficient of the parallelepiped is calculated based on
74the averaged effective radius, after appropriately sorting the three
75dimensions, to give an oblate or prolate particle, $(=\sqrt{AB/\pi})$ and
76length $(= C)$ values, and used as the effective radius for
77$S(q)$ when $P(q) \cdot S(q)$ is applied.
78
79For 2d data the orientation of the particle is required, described using
80angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details
81of the calculation and angular dispersions see :ref:orientation .
82
83.. Comment by Miguel Gonzalez:
84   The following text has been commented because I think there are two
85   mistakes. Psi is the rotational angle around C (but I cannot understand
86   what it means against the q plane) and psi=0 corresponds to a||x and b||y.
87
88   The angle $\Psi$ is the rotational angle around the $C$ axis against
89   the $q$ plane. For example, $\Psi = 0$ when the $B$ axis is parallel
90   to the $x$-axis of the detector.
91
92The angle $\Psi$ is the rotational angle around the $C$ axis.
93For $\theta = 0$ and $\phi = 0$, $\Psi = 0$ corresponds to the $B$ axis
94oriented parallel to the y-axis of the detector with $A$ along the x-axis.
95For other $\theta$, $\phi$ values, the parallelepiped has to be first rotated
96$\theta$ degrees in the $z-x$ plane and then $\phi$ degrees around the $z$ axis,
97before doing a final rotation of $\Psi$ degrees around the resulting $C$ axis
98of the particle to obtain the final orientation of the parallelepiped.
99
100.. _parallelepiped-orientation:
101
102.. figure:: img/parallelepiped_angle_definition.png
103
104    Definition of the angles for oriented parallelepiped, shown with $A<B<C$.
105
106.. figure:: img/parallelepiped_angle_projection.png
107
108    Examples of the angles for an oriented parallelepiped against the
109    detector plane.
110
111On introducing "Orientational Distribution" in the angles, "distribution of
112theta" and "distribution of phi" parameters will appear. These are actually
113rotations about axes $\delta_1$ and $\delta_2$ of the parallelepiped,
114perpendicular to the $a$ x $c$ and $b$ x $c$ faces. (When $\theta = \phi = 0$
115these are parallel to the $Y$ and $X$ axes of the instrument.) The third
116orientation distribution, in $\psi$, is about the $c$ axis of the particle,
117perpendicular to the $a$ x $b$ face. Some experimentation may be required to
118understand the 2d patterns fully as discussed in :ref:orientation .
119
120For a given orientation of the parallelepiped, the 2D form factor is
121calculated as
122
123.. math::
124
125    P(q_x, q_y) = \left[\frac{\sin(\tfrac{1}{2}qA\cos\alpha)}{(\tfrac{1}
126                   {2}qA\cos\alpha)}\right]^2
127                  \left[\frac{\sin(\tfrac{1}{2}qB\cos\beta)}{(\tfrac{1}
128                   {2}qB\cos\beta)}\right]^2
129                  \left[\frac{\sin(\tfrac{1}{2}qC\cos\gamma)}{(\tfrac{1}
130                   {2}qC\cos\gamma)}\right]^2
131
132with
133
134.. math::
135
136    \cos\alpha &= \hat A \cdot \hat q, \\
137    \cos\beta  &= \hat B \cdot \hat q, \\
138    \cos\gamma &= \hat C \cdot \hat q
139
140and the scattering intensity as:
141
142.. math::
143
144    I(q_x, q_y) = \frac{\text{scale}}{V} V^2 \Delta\rho^2 P(q_x, q_y)
145            + \text{background}
146
147.. Comment by Miguel Gonzalez:
148   This reflects the logic of the code, as in parallelepiped.c the call
149   to _pkernel returns $P(q_x, q_y)$ and then this is multiplied by
150   $V^2 * (\Delta \rho)^2$. And finally outside parallelepiped.c it will be
151   multiplied by scale, normalized by $V$ and the background added. But
152   mathematically it makes more sense to write
153   $I(q_x, q_y) = \text{scale} V \Delta\rho^2 P(q_x, q_y) + \text{background}$,
154   with scale being the volume fraction.
155
156
157Validation
158----------
159
160Validation of the code was done by comparing the output of the 1D calculation
161to the angular average of the output of a 2D calculation over all possible
162angles.
163
164
165References
166----------
167
168.. [#Mittelbach] P Mittelbach and G Porod, *Acta Physica Austriaca*,
169   14 (1961) 185-211
170.. [#] R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854
171
172Authorship and Verification
173----------------------------
174
175* **Author:** NIST IGOR/DANSE **Date:** pre 2010
177* **Last Reviewed by:**  Richard Heenan **Date:** April 06, 2017
178"""
179
180import numpy as np
181from numpy import pi, inf, sqrt, sin, cos
182
183name = "parallelepiped"
184title = "Rectangular parallelepiped with uniform scattering length density."
185description = """
186    I(q)= scale*V*(sld - sld_solvent)^2*P(q,alpha)+background
187        P(q,alpha) = integral from 0 to 1 of ...
188           phi(mu*sqrt(1-sigma^2),a) * S(mu*c*sigma/2)^2 * dsigma
189        with
190            phi(mu,a) = integral from 0 to 1 of ..
191            (S((mu/2)*cos(pi*u/2))*S((mu*a/2)*sin(pi*u/2)))^2 * du
192            S(x) = sin(x)/x
193            mu = q*B
194        V: Volume of the rectangular parallelepiped
195        alpha: angle between the long axis of the
196            parallelepiped and the q-vector for 1D
197"""
198category = "shape:parallelepiped"
199
200#             ["name", "units", default, [lower, upper], "type","description"],
201parameters = [["sld", "1e-6/Ang^2", 4, [-inf, inf], "sld",
202               "Parallelepiped scattering length density"],
203              ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "sld",
204               "Solvent scattering length density"],
205              ["length_a", "Ang", 35, [0, inf], "volume",
206               "Shorter side of the parallelepiped"],
207              ["length_b", "Ang", 75, [0, inf], "volume",
208               "Second side of the parallelepiped"],
209              ["length_c", "Ang", 400, [0, inf], "volume",
210               "Larger side of the parallelepiped"],
211              ["theta", "degrees", 60, [-360, 360], "orientation",
212               "c axis to beam angle"],
213              ["phi", "degrees", 60, [-360, 360], "orientation",
215              ["psi", "degrees", 60, [-360, 360], "orientation",
217             ]
218
219source = ["lib/gauss76.c", "parallelepiped.c"]
220
221def ER(length_a, length_b, length_c):
222    """
223    Return effective radius (ER) for P(q)*S(q)
224    """
225    # now that axes can be in any size order, need to sort a,b,c
226    # where a~b and c is either much smaller or much larger
227    abc = np.vstack((length_a, length_b, length_c))
228    abc = np.sort(abc, axis=0)
229    selector = (abc[1] - abc[0]) > (abc[2] - abc[1])
230    length = np.where(selector, abc[0], abc[2])
231    # surface average radius (rough approximation)
232    radius = sqrt(np.where(~selector, abc[0]*abc[1], abc[1]*abc[2]) / pi)
233
235    return 0.5 * (ddd) ** (1. / 3.)
236
237# VR defaults to 1.0
238
239
240def random():
241    length = 10**np.random.uniform(1, 4.7, size=3)
242    pars = dict(
243        length_a=length[0],
244        length_b=length[1],
245        length_c=length[2],
246    )
247    return pars
248
249
250# parameters for demo
251demo = dict(scale=1, background=0,
252            sld=6.3, sld_solvent=1.0,
253            length_a=35, length_b=75, length_c=400,
254            theta=45, phi=30, psi=15,
255            length_a_pd=0.1, length_a_pd_n=10,
256            length_b_pd=0.1, length_b_pd_n=1,
257            length_c_pd=0.1, length_c_pd_n=1,
258            theta_pd=10, theta_pd_n=1,
259            phi_pd=10, phi_pd_n=1,
260            psi_pd=10, psi_pd_n=10)
261# rkh 7/4/17 add random unit test for 2d, note make all params different,
262# 2d values not tested against other codes or models
263qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.)
264tests = [[{}, 0.2, 0.17758004974],
265         [{}, [0.2], [0.17758004974]],
266         [{'theta':10.0, 'phi':20.0}, (qx, qy), 0.0089517140475],
267         [{'theta':10.0, 'phi':20.0}, [(qx, qy)], [0.0089517140475]],
268        ]
269del qx, qy  # not necessary to delete, but cleaner
Note: See TracBrowser for help on using the repository browser.