source: sasmodels/sasmodels/models/parallelepiped.py @ ef07e95

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since ef07e95 was ef07e95, checked in by Paul Kienzle <pkienzle@…>, 19 months ago

use consistent formatting for 'Last Reviewed'

  • Property mode set to 100644
File size: 9.9 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:** NIST IGOR/DANSE **Date:** pre 2010
170* **Last Modified by:**  Paul Kienzle **Date:** April 05, 2017
171* **Last Reviewed by:**  Richard Heenan **Date:** April 06, 2017
172"""
173
174import numpy as np
175from numpy import pi, inf, sqrt, sin, cos
176
177name = "parallelepiped"
178title = "Rectangular parallelepiped with uniform scattering length density."
179description = """
180    I(q)= scale*V*(sld - sld_solvent)^2*P(q,alpha)+background
181        P(q,alpha) = integral from 0 to 1 of ...
182           phi(mu*sqrt(1-sigma^2),a) * S(mu*c*sigma/2)^2 * dsigma
183        with
184            phi(mu,a) = integral from 0 to 1 of ..
185            (S((mu/2)*cos(pi*u/2))*S((mu*a/2)*sin(pi*u/2)))^2 * du
186            S(x) = sin(x)/x
187            mu = q*B
188        V: Volume of the rectangular parallelepiped
189        alpha: angle between the long axis of the
190            parallelepiped and the q-vector for 1D
191"""
192category = "shape:parallelepiped"
193
194#             ["name", "units", default, [lower, upper], "type","description"],
195parameters = [["sld", "1e-6/Ang^2", 4, [-inf, inf], "sld",
196               "Parallelepiped scattering length density"],
197              ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "sld",
198               "Solvent scattering length density"],
199              ["length_a", "Ang", 35, [0, inf], "volume",
200               "Shorter side of the parallelepiped"],
201              ["length_b", "Ang", 75, [0, inf], "volume",
202               "Second side of the parallelepiped"],
203              ["length_c", "Ang", 400, [0, inf], "volume",
204               "Larger side of the parallelepiped"],
205              ["theta", "degrees", 60, [-360, 360], "orientation",
206               "c axis to beam angle"],
207              ["phi", "degrees", 60, [-360, 360], "orientation",
208               "rotation about beam"],
209              ["psi", "degrees", 60, [-360, 360], "orientation",
210               "rotation about c axis"],
211             ]
212
213source = ["lib/gauss76.c", "parallelepiped.c"]
214
215def ER(length_a, length_b, length_c):
216    """
217    Return effective radius (ER) for P(q)*S(q)
218    """
219    # now that axes can be in any size order, need to sort a,b,c
220    # where a~b and c is either much smaller or much larger
221    abc = np.vstack((length_a, length_b, length_c))
222    abc = np.sort(abc, axis=0)
223    selector = (abc[1] - abc[0]) > (abc[2] - abc[1])
224    length = np.where(selector, abc[0], abc[2])
225    # surface average radius (rough approximation)
226    radius = sqrt(np.where(~selector, abc[0]*abc[1], abc[1]*abc[2]) / pi)
227
228    ddd = 0.75 * radius * (2*radius*length + (length + radius)*(length + pi*radius))
229    return 0.5 * (ddd) ** (1. / 3.)
230
231# VR defaults to 1.0
232
233
234def random():
235    length = 10**np.random.uniform(1, 4.7, size=3)
236    pars = dict(
237        length_a=length[0],
238        length_b=length[1],
239        length_c=length[2],
240    )
241    return pars
242
243
244# parameters for demo
245demo = dict(scale=1, background=0,
246            sld=6.3, sld_solvent=1.0,
247            length_a=35, length_b=75, length_c=400,
248            theta=45, phi=30, psi=15,
249            length_a_pd=0.1, length_a_pd_n=10,
250            length_b_pd=0.1, length_b_pd_n=1,
251            length_c_pd=0.1, length_c_pd_n=1,
252            theta_pd=10, theta_pd_n=1,
253            phi_pd=10, phi_pd_n=1,
254            psi_pd=10, psi_pd_n=10)
255# rkh 7/4/17 add random unit test for 2d, note make all params different,
256# 2d values not tested against other codes or models
257qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.)
258tests = [[{}, 0.2, 0.17758004974],
259         [{}, [0.2], [0.17758004974]],
260         [{'theta':10.0, 'phi':20.0}, (qx, qy), 0.0089517140475],
261         [{'theta':10.0, 'phi':20.0}, [(qx, qy)], [0.0089517140475]],
262        ]
263del qx, qy  # not necessary to delete, but cleaner
Note: See TracBrowser for help on using the repository browser.