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

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since b343226 was b343226, checked in by Paul Kienzle <pkienzle@…>, 6 years ago

trivial doc formatting

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