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

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

Update Verification section

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