source: sasmodels/sasmodels/models/parallelepiped.py @ 9802ab3

core_shell_microgelscostrafo411magnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 9802ab3 was 9802ab3, checked in by richardh, 7 years ago

changes to docs, anticipating new orientation angle integrals

  • 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
22.. note::
23
24   The edge of the solid used to have to satisfy the condition that $A < B < C$.
25   After some improvements to the effective radius calculation, used with
26   an S(Q), it is beleived that this is no longer the case.
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) \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}$,
45$P(q, \alpha)$ is the form factor corresponding to a parallelepiped oriented
46at an angle $\alpha$ (angle between the long axis C and $\vec q$),
47and the averaging $\left<\ldots\right>$ is applied over all orientations.
48
49Assuming $a = A/B < 1$, $b = B /B = 1$, and $c = C/B > 1$, the
50form factor is given by (Mittelbach and Porod, 1961)
51
52.. math::
53
54    P(q, \alpha) = \int_0^1 \phi_Q\left(\mu \sqrt{1-\sigma^2},a\right)
55        \left[S(\mu c \sigma/2)\right]^2 d\sigma
56
57with
58
59.. math::
60
61    \phi_Q(\mu,a) &= \int_0^1
62        \left\{S\left[\frac{\mu}{2}\cos\left(\frac{\pi}{2}u\right)\right]
63               S\left[\frac{\mu a}{2}\sin\left(\frac{\pi}{2}u\right)\right]
64               \right\}^2 du
65
66    S(x) &= \frac{\sin x}{x}
67
68    \mu &= qB
69
70The scattering intensity per unit volume is returned in units of |cm^-1|.
71
72NB: The 2nd virial coefficient of the parallelepiped is calculated based on
73the averaged effective radius, after appropriately sorting the three
74dimensions, to give an oblate or prolate particle, $(=\sqrt{AB/\pi})$ and
75length $(= C)$ values, and used as the effective radius for
76$S(q)$ when $P(q) \cdot S(q)$ is applied.
77
78To provide easy access to the orientation of the parallelepiped, we define
79three angles $\theta$, $\phi$ and $\Psi$. The definition of $\theta$ and
80$\phi$ is the same as for the cylinder model (see also figures below).
81
82.. Comment by Miguel Gonzalez:
83   The following text has been commented because I think there are two
84   mistakes. Psi is the rotational angle around C (but I cannot understand
85   what it means against the q plane) and psi=0 corresponds to a||x and b||y.
86
87   The angle $\Psi$ is the rotational angle around the $C$ axis against
88   the $q$ plane. For example, $\Psi = 0$ when the $B$ axis is parallel
89   to the $x$-axis of the detector.
90
91The angle $\Psi$ is the rotational angle around the $C$ axis.
92For $\theta = 0$ and $\phi = 0$, $\Psi = 0$ corresponds to the $B$ axis
93oriented parallel to the y-axis of the detector with $A$ along the z-axis.
94For other $\theta$, $\phi$ values, the parallelepiped has to be first rotated
95$\theta$ degrees around $z$ and $\phi$ degrees around $y$,
96before doing a final rotation of $\Psi$ degrees around the resulting $C$ to
97obtain the final orientation of the parallelepiped.
98For example, for $\theta = 0$ and $\phi = 90$, we have that $\Psi = 0$
99corresponds to $A$ along $x$ and $B$ along $y$,
100while for $\theta = 90$ and $\phi = 0$, $\Psi = 0$ corresponds to
101$A$ along $z$ and $B$ along $x$.
102
103.. _parallelepiped-orientation:
104
105.. figure:: img/parallelepiped_angle_definition.png
106
107    Definition of the angles for oriented parallelepiped, shown with $A<B<C$.
108
109.. figure:: img/parallelepiped_angle_projection.png
110
111    Examples of the angles for an oriented parallelepiped against the
112    detector plane.
113
114On introducing "Orientational Distribution" in the angles, "distribution of theta" and "distribution of phi" parameters will
115appear. These are actually rotations about axes $\delta_1$ and $\delta_2$ of the parallelepiped, perpendicular to the $a$ x $c$ and $b$ x $c$ faces.
116(When $\theta = \phi = 0$ these are parallel to the $Y$ and $X$ axes of the instrument.) The third orientation distribution, in $\psi$, is
117about the $c$ axis of the particle, perpendicular to the $a$ x $b$ face. Some experimentation may be required to
118understand the 2d patterns fully. (Earlier implementations had numerical integration issues in some circumstances when orientation
119distributions passed through 90 degrees, such situations, with very broad distributions, should still be approached with care.)
120
121   
122For a given orientation of the parallelepiped, the 2D form factor is
123calculated as
124
125.. math::
126
127    P(q_x, q_y) = \left[\frac{\sin(\tfrac{1}{2}qA\cos\alpha)}{(\tfrac{1}{2}qA\cos\alpha)}\right]^2
128                  \left[\frac{\sin(\tfrac{1}{2}qB\cos\beta)}{(\tfrac{1}{2}qB\cos\beta)}\right]^2
129                  \left[\frac{\sin(\tfrac{1}{2}qC\cos\gamma)}{(\tfrac{1}{2}qC\cos\gamma)}\right]^2
130
131with
132
133.. math::
134
135    \cos\alpha &= \hat A \cdot \hat q,
136
137    \cos\beta  &= \hat B \cdot \hat q,
138
139    \cos\gamma &= \hat C \cdot \hat q
140
141and the scattering intensity as:
142
143.. math::
144
145    I(q_x, q_y) = \frac{\text{scale}}{V} V^2 \Delta\rho^2 P(q_x, q_y)
146            + \text{background}
147
148.. Comment by Miguel Gonzalez:
149   This reflects the logic of the code, as in parallelepiped.c the call
150   to _pkernel returns $P(q_x, q_y)$ and then this is multiplied by
151   $V^2 * (\Delta \rho)^2$. And finally outside parallelepiped.c it will be
152   multiplied by scale, normalized by $V$ and the background added. But
153   mathematically it makes more sense to write
154   $I(q_x, q_y) = \text{scale} V \Delta\rho^2 P(q_x, q_y) + \text{background}$,
155   with scale being the volume fraction.
156
157
158Validation
159----------
160
161Validation of the code was done by comparing the output of the 1D calculation
162to the angular average of the output of a 2D calculation over all possible
163angles.
164
165
166References
167----------
168
169P Mittelbach and G Porod, *Acta Physica Austriaca*, 14 (1961) 185-211
170
171R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854
172
173Authorship and Verification
174----------------------------
175
176* **Author:** This model is based on form factor calculations implemented
177    in a c-library provided by the NIST Center for Neutron Research (Kline, 2006).
178* **Last Modified by:**  Paul Kienzle **Date:** April 05, 2017
179* **Last Reviewed by:**  Richard Heenan **Date:** April 06, 2017
180
181"""
182
183import numpy as np
184from numpy import pi, inf, sqrt, sin, cos
185
186name = "parallelepiped"
187title = "Rectangular parallelepiped with uniform scattering length density."
188description = """
189    I(q)= scale*V*(sld - sld_solvent)^2*P(q,alpha)+background
190        P(q,alpha) = integral from 0 to 1 of ...
191           phi(mu*sqrt(1-sigma^2),a) * S(mu*c*sigma/2)^2 * dsigma
192        with
193            phi(mu,a) = integral from 0 to 1 of ..
194            (S((mu/2)*cos(pi*u/2))*S((mu*a/2)*sin(pi*u/2)))^2 * du
195            S(x) = sin(x)/x
196            mu = q*B
197        V: Volume of the rectangular parallelepiped
198        alpha: angle between the long axis of the
199            parallelepiped and the q-vector for 1D
200"""
201category = "shape:parallelepiped"
202
203#             ["name", "units", default, [lower, upper], "type","description"],
204parameters = [["sld", "1e-6/Ang^2", 4, [-inf, inf], "sld",
205               "Parallelepiped scattering length density"],
206              ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "sld",
207               "Solvent scattering length density"],
208              ["length_a", "Ang", 35, [0, inf], "volume",
209               "Shorter side of the parallelepiped"],
210              ["length_b", "Ang", 75, [0, inf], "volume",
211               "Second side of the parallelepiped"],
212              ["length_c", "Ang", 400, [0, inf], "volume",
213               "Larger side of the parallelepiped"],
214              ["theta", "degrees", 60, [-360, 360], "orientation",
215               "c axis to beam angle"],
216              ["phi", "degrees", 60, [-360, 360], "orientation",
217               "rotation about beam"],
218              ["psi", "degrees", 60, [-360, 360], "orientation",
219               "rotation about c axis"],
220             ]
221
222source = ["lib/gauss76.c", "parallelepiped.c"]
223
224def ER(length_a, length_b, length_c):
225    """
226    Return effective radius (ER) for P(q)*S(q)
227    """
228    # now that axes can be in any size order, need to sort a,b,c where a~b and c is either much smaller
229    # or much larger
230    abc = np.vstack((length_a, length_b, length_c))
231    abc = np.sort(abc, axis=0)
232    selector = (abc[1] - abc[0]) > (abc[2] - abc[1])
233    length = np.where(selector, abc[0], abc[2])
234    # surface average radius (rough approximation)
235    radius = np.sqrt(np.where(~selector, abc[0]*abc[1], abc[1]*abc[2]) / pi)
236
237    ddd = 0.75 * radius * (2*radius*length + (length + radius)*(length + pi*radius))
238    return 0.5 * (ddd) ** (1. / 3.)
239
240# VR defaults to 1.0
241
242# parameters for demo
243demo = dict(scale=1, background=0,
244            sld=6.3, sld_solvent=1.0,
245            length_a=35, length_b=75, length_c=400,
246            theta=45, phi=30, psi=15,
247            length_a_pd=0.1, length_a_pd_n=10,
248            length_b_pd=0.1, length_b_pd_n=1,
249            length_c_pd=0.1, length_c_pd_n=1,
250            theta_pd=10, theta_pd_n=1,
251            phi_pd=10, phi_pd_n=1,
252            psi_pd=10, psi_pd_n=10)
253# rkh 7/4/17 add random unit test for 2d, note make all params different, 2d values not tested against other codes or models
254qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.)
255tests = [[{}, 0.2, 0.17758004974],
256         [{}, [0.2], [0.17758004974]],
257         [{'theta':10.0, 'phi':20.0}, (qx, qy), 0.0089517140475],
258         [{'theta':10.0, 'phi':20.0}, [(qx, qy)], [0.0089517140475]],
259        ]
260del qx, qy  # not necessary to delete, but cleaner
Note: See TracBrowser for help on using the repository browser.