source: sasmodels/sasmodels/models/parallelepiped.py @ 96153e4

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 96153e4 was 96153e4, checked in by butler, 14 months ago

more corrections and normalizations

Addreses #896. Brings parallelepiped and core-shell parallelepiped more
into line with each other and corrects angle definitions refering to
orietnational distribution docs for details. Still need to sort out
with Paul Kienzle the proper order the angles are applied.

  • Property mode set to 100644
File size: 10.1 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 parallelepiped
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
74NB: The 2nd virial coefficient of the parallelepiped is calculated based on
75the averaged effective radius, after appropriately sorting the three
76dimensions, to give an oblate or prolate particle, $(=\sqrt{AB/\pi})$ and
77length $(= C)$ values, and used as the effective radius for
78$S(q)$ when $P(q) \cdot S(q)$ is applied.
79
80For 2d data the orientation of the particle is required, described using
81angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details
82of the calculation and angular dispersions see :ref:`orientation` .
83
84The angle $\Psi$ is the rotational angle around the $C$ axis.
85For $\theta = 0$ and $\phi = 0$, $\Psi = 0$ corresponds to the $B$ axis
86oriented parallel to the y-axis of the detector with $A$ along the x-axis.
87For other $\theta$, $\phi$ values, the parallelepiped has to be first rotated
88$\theta$ degrees in the $z-x$ plane and then $\phi$ degrees around the $z$ axis,
89before doing a final rotation of $\Psi$ degrees around the resulting $C$ axis
90of the particle to obtain the final orientation of the parallelepiped.
91
92.. note:: For 2d, constraints must be applied during fitting to ensure that the
93   inequality $A < B < C$ is not violated, and hence the correct definition
94   of angles is preserved. The calculation will not report an error,
95   but the results may be not correct.
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}
123                   {2}qA\cos\alpha)}\right]^2
124                  \left[\frac{\sin(\tfrac{1}{2}qB\cos\beta)}{(\tfrac{1}
125                   {2}qB\cos\beta)}\right]^2
126                  \left[\frac{\sin(\tfrac{1}{2}qC\cos\gamma)}{(\tfrac{1}
127                   {2}qC\cos\gamma)}\right]^2
128
129with
130
131.. math::
132
133    \cos\alpha &= \hat A \cdot \hat q, \\
134    \cos\beta  &= \hat B \cdot \hat q, \\
135    \cos\gamma &= \hat C \cdot \hat q
136
137and the scattering intensity as:
138
139.. math::
140
141    I(q_x, q_y) = \frac{\text{scale}}{V} V^2 \Delta\rho^2 P(q_x, q_y)
142            + \text{background}
143
144.. Comment by Miguel Gonzalez:
145   This reflects the logic of the code, as in parallelepiped.c the call
146   to _pkernel returns $P(q_x, q_y)$ and then this is multiplied by
147   $V^2 * (\Delta \rho)^2$. And finally outside parallelepiped.c it will be
148   multiplied by scale, normalized by $V$ and the background added. But
149   mathematically it makes more sense to write
150   $I(q_x, q_y) = \text{scale} V \Delta\rho^2 P(q_x, q_y) + \text{background}$,
151   with scale being the volume fraction.
152
153
154Validation
155----------
156
157Validation of the code was done by comparing the output of the 1D calculation
158to the angular average of the output of a 2D calculation over all possible
159angles.
160
161
162References
163----------
164
165.. [#Mittelbach] P Mittelbach and G Porod, *Acta Physica Austriaca*,
166   14 (1961) 185-211
167.. [#] R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854
168
169Authorship and Verification
170----------------------------
171
172* **Author:** NIST IGOR/DANSE **Date:** pre 2010
173* **Last Modified by:**  Paul Kienzle **Date:** April 05, 2017
174* **Last Reviewed by:**  Richard Heenan **Date:** April 06, 2017
175"""
176
177import numpy as np
178from numpy import pi, inf, sqrt, sin, cos
179
180name = "parallelepiped"
181title = "Rectangular parallelepiped with uniform scattering length density."
182description = """
183    I(q)= scale*V*(sld - sld_solvent)^2*P(q,alpha)+background
184        P(q,alpha) = integral from 0 to 1 of ...
185           phi(mu*sqrt(1-sigma^2),a) * S(mu*c*sigma/2)^2 * dsigma
186        with
187            phi(mu,a) = integral from 0 to 1 of ..
188            (S((mu/2)*cos(pi*u/2))*S((mu*a/2)*sin(pi*u/2)))^2 * du
189            S(x) = sin(x)/x
190            mu = q*B
191        V: Volume of the rectangular parallelepiped
192        alpha: angle between the long axis of the
193            parallelepiped and the q-vector for 1D
194"""
195category = "shape:parallelepiped"
196
197#             ["name", "units", default, [lower, upper], "type","description"],
198parameters = [["sld", "1e-6/Ang^2", 4, [-inf, inf], "sld",
199               "Parallelepiped scattering length density"],
200              ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "sld",
201               "Solvent scattering length density"],
202              ["length_a", "Ang", 35, [0, inf], "volume",
203               "Shorter side of the parallelepiped"],
204              ["length_b", "Ang", 75, [0, inf], "volume",
205               "Second side of the parallelepiped"],
206              ["length_c", "Ang", 400, [0, inf], "volume",
207               "Larger side of the parallelepiped"],
208              ["theta", "degrees", 60, [-360, 360], "orientation",
209               "c axis to beam angle"],
210              ["phi", "degrees", 60, [-360, 360], "orientation",
211               "rotation about beam"],
212              ["psi", "degrees", 60, [-360, 360], "orientation",
213               "rotation about c axis"],
214             ]
215
216source = ["lib/gauss76.c", "parallelepiped.c"]
217
218def ER(length_a, length_b, length_c):
219    """
220    Return effective radius (ER) for P(q)*S(q)
221    """
222    # now that axes can be in any size order, need to sort a,b,c
223    # where a~b and c is either much smaller or much larger
224    abc = np.vstack((length_a, length_b, length_c))
225    abc = np.sort(abc, axis=0)
226    selector = (abc[1] - abc[0]) > (abc[2] - abc[1])
227    length = np.where(selector, abc[0], abc[2])
228    # surface average radius (rough approximation)
229    radius = sqrt(np.where(~selector, abc[0]*abc[1], abc[1]*abc[2]) / pi)
230
231    ddd = 0.75 * radius * (2*radius*length + (length + radius)*(length + pi*radius))
232    return 0.5 * (ddd) ** (1. / 3.)
233
234# VR defaults to 1.0
235
236
237def random():
238    length = 10**np.random.uniform(1, 4.7, size=3)
239    pars = dict(
240        length_a=length[0],
241        length_b=length[1],
242        length_c=length[2],
243    )
244    return pars
245
246
247# parameters for demo
248demo = dict(scale=1, background=0,
249            sld=6.3, sld_solvent=1.0,
250            length_a=35, length_b=75, length_c=400,
251            theta=45, phi=30, psi=15,
252            length_a_pd=0.1, length_a_pd_n=10,
253            length_b_pd=0.1, length_b_pd_n=1,
254            length_c_pd=0.1, length_c_pd_n=1,
255            theta_pd=10, theta_pd_n=1,
256            phi_pd=10, phi_pd_n=1,
257            psi_pd=10, psi_pd_n=10)
258# rkh 7/4/17 add random unit test for 2d, note make all params different,
259# 2d values not tested against other codes or models
260qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.)
261tests = [[{}, 0.2, 0.17758004974],
262         [{}, [0.2], [0.17758004974]],
263         [{'theta':10.0, 'phi':20.0}, (qx, qy), 0.0089517140475],
264         [{'theta':10.0, 'phi':20.0}, [(qx, qy)], [0.0089517140475]],
265        ]
266del qx, qy  # not necessary to delete, but cleaner
Note: See TracBrowser for help on using the repository browser.