source: sasmodels/sasmodels/models/parallelepiped.py @ 2d81cfe

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

lint

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