source: sasmodels/sasmodels/models/parallelepiped.py @ 830cf6b

ticket-1257-vesicle-productticket_1156ticket_822_more_unit_tests
Last change on this file since 830cf6b was a34b811, checked in by Paul Kienzle <pkienzle@…>, 5 years ago

use radius_effective/radius_effective_mode/radius_effective_modes consistently throughout the code

  • 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.. [#] L. Onsager, *Ann. New York Acad. Sci.*, 51 (1949) 627-659
183
184Source
185------
186
187`parallelepiped.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/parallelepiped.py>`_
188
189`parallelepiped.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/parallelepiped.c>`_
190
191Authorship and Verification
192----------------------------
193
194* **Author:** NIST IGOR/DANSE **Date:** pre 2010
195* **Last Modified by:**  Paul Kienzle **Date:** April 05, 2017
196* **Last Reviewed by:**  Miguel Gonzales and Paul Butler **Date:** May 24,
197  2018 - documentation updated
198* **Source added by :** Steve King **Date:** March 25, 2019
199"""
200
201import numpy as np
202from numpy import inf
203
204name = "parallelepiped"
205title = "Rectangular parallelepiped with uniform scattering length density."
206description = """
207    I(q)= scale*V*(sld - sld_solvent)^2*P(q,alpha)+background
208        P(q,alpha) = integral from 0 to 1 of ...
209           phi(mu*sqrt(1-sigma^2),a) * S(mu*c*sigma/2)^2 * dsigma
210        with
211            phi(mu,a) = integral from 0 to 1 of ..
212            (S((mu/2)*cos(pi*u/2))*S((mu*a/2)*sin(pi*u/2)))^2 * du
213            S(x) = sin(x)/x
214            mu = q*B
215        V: Volume of the rectangular parallelepiped
216        alpha: angle between the long axis of the
217            parallelepiped and the q-vector for 1D
218"""
219category = "shape:parallelepiped"
220
221#             ["name", "units", default, [lower, upper], "type","description"],
222parameters = [["sld", "1e-6/Ang^2", 4, [-inf, inf], "sld",
223               "Parallelepiped scattering length density"],
224              ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "sld",
225               "Solvent scattering length density"],
226              ["length_a", "Ang", 35, [0, inf], "volume",
227               "Shorter side of the parallelepiped"],
228              ["length_b", "Ang", 75, [0, inf], "volume",
229               "Second side of the parallelepiped"],
230              ["length_c", "Ang", 400, [0, inf], "volume",
231               "Larger side of the parallelepiped"],
232              ["theta", "degrees", 60, [-360, 360], "orientation",
233               "c axis to beam angle"],
234              ["phi", "degrees", 60, [-360, 360], "orientation",
235               "rotation about beam"],
236              ["psi", "degrees", 60, [-360, 360], "orientation",
237               "rotation about c axis"],
238             ]
239
240source = ["lib/gauss76.c", "parallelepiped.c"]
241have_Fq = True
242radius_effective_modes = [
243    "equivalent cylinder excluded volume", "equivalent volume sphere",
244    "half length_a", "half length_b", "half length_c",
245    "equivalent circular cross-section", "half ab diagonal", "half diagonal",
246    ]
247
248def random():
249    """Return a random parameter set for the model."""
250    length = 10**np.random.uniform(1, 4.7, size=3)
251    pars = dict(
252        length_a=length[0],
253        length_b=length[1],
254        length_c=length[2],
255    )
256    return pars
257
258
259# parameters for demo
260demo = dict(scale=1, background=0,
261            sld=6.3, sld_solvent=1.0,
262            length_a=35, length_b=75, length_c=400,
263            theta=45, phi=30, psi=15,
264            length_a_pd=0.1, length_a_pd_n=10,
265            length_b_pd=0.1, length_b_pd_n=1,
266            length_c_pd=0.1, length_c_pd_n=1,
267            theta_pd=10, theta_pd_n=1,
268            phi_pd=10, phi_pd_n=1,
269            psi_pd=10, psi_pd_n=10)
270# rkh 7/4/17 add random unit test for 2d, note make all params different,
271# 2d values not tested against other codes or models
272qx, qy = 0.2 * np.cos(np.pi/6.), 0.2 * np.sin(np.pi/6.)
273tests = [[{}, 0.2, 0.17758004974],
274         [{}, [0.2], [0.17758004974]],
275         [{'theta':10.0, 'phi':20.0}, (qx, qy), 0.0089517140475],
276         [{'theta':10.0, 'phi':20.0}, [(qx, qy)], [0.0089517140475]],
277        ]
278del qx, qy  # not necessary to delete, but cleaner
Note: See TracBrowser for help on using the repository browser.