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

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

Final? edits to Parallelepiped (and core shell version) documentation

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