source: sasmodels/sasmodels/models/parallelepiped.py @ 99658f6

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 99658f6 was 99658f6, checked in by grethevj, 10 months ago

updated ER functions including cylinder excluded volume, to match 4.x

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