source: sasmodels/sasmodels/models/rectangular_prism.py @ d57b06c

Last change on this file since d57b06c was d57b06c, checked in by Paul Kienzle <pkienzle@…>, 7 months ago

Merge remote-tracking branch 'origin/master' into ticket-1263-source-link

  • Property mode set to 100644
File size: 6.4 KB
Line 
1# rectangular_prism model
2# Note: model title and parameter table are inserted automatically
3r"""
4
5This model provides the form factor, $P(q)$, for a rectangular prism.
6
7Note that this model is almost totally equivalent to the existing
8:ref:`parallelepiped` model.
9The only difference is that the way the relevant
10parameters are defined here ($a$, $b/a$, $c/a$ instead of $a$, $b$, $c$)
11which allows use of polydispersity with this model while keeping the shape of
12the prism (e.g. setting $b/a = 1$ and $c/a = 1$ and applying polydispersity
13to *a* will generate a distribution of cubes of different sizes).
14
15
16Definition
17----------
18
19The 1D scattering intensity for this model was calculated by Mittelbach and
20Porod (Mittelbach, 1961), but the implementation here is closer to the
21equations given by Nayuk and Huber (Nayuk, 2012).
22Note also that the angle definitions used in the code and the present
23documentation correspond to those used in (Nayuk, 2012) (see Fig. 1 of
24that reference), with $\theta$ corresponding to $\alpha$ in that paper,
25and not to the usual convention used for example in the
26:ref:`parallelepiped` model.
27
28In this model the scattering from a massive parallelepiped with an
29orientation with respect to the scattering vector given by $\theta$
30and $\phi$
31
32.. math::
33
34    A_P\,(q) =
35      \frac{\sin \left( \tfrac{1}{2}qC \cos\theta \right) }{\tfrac{1}{2} qC \cos\theta}
36      \,\times\,
37      \frac{\sin \left( \tfrac{1}{2}qA \cos\theta \right) }{\tfrac{1}{2} qA \cos\theta}
38      \,\times\ ,
39      \frac{\sin \left( \tfrac{1}{2}qB \cos\theta \right) }{\tfrac{1}{2} qB \cos\theta}
40
41where $A$, $B$ and $C$ are the sides of the parallelepiped and must fulfill
42$A \le B \le C$, $\theta$ is the angle between the $z$ axis and the
43longest axis of the parallelepiped $C$, and $\phi$ is the angle between the
44scattering vector (lying in the $xy$ plane) and the $y$ axis.
45
46The normalized form factor in 1D is obtained averaging over all possible
47orientations
48
49.. math::
50    P(q) =  \frac{2}{\pi} \int_0^{\frac{\pi}{2}} \,
51    \int_0^{\frac{\pi}{2}} A_P^2(q) \, \sin\theta \, d\theta \, d\phi
52
53And the 1D scattering intensity is calculated as
54
55.. math::
56    I(q) = \text{scale} \times V \times (\rho_\text{p} -
57    \rho_\text{solvent})^2 \times P(q)
58
59where $V$ is the volume of the rectangular prism, $\rho_\text{p}$
60is the scattering length of the parallelepiped, $\rho_\text{solvent}$
61is the scattering length of the solvent, and (if the data are in absolute
62units) *scale* represents the volume fraction (which is unitless).
63
64For 2d data the orientation of the particle is required, described using
65angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details
66of the calculation and angular dispersions see :ref:`orientation` .
67The angle $\Psi$ is the rotational angle around the long *C* axis. For example,
68$\Psi = 0$ when the *B* axis is parallel to the *x*-axis of the detector.
69
70For 2d, constraints must be applied during fitting to ensure that the inequality
71$A < B < C$ is not violated, and hence the correct definition of angles is preserved. The calculation will not report an error,
72but the results may be not correct.
73
74.. figure:: img/parallelepiped_angle_definition.png
75
76    Definition of the angles for oriented core-shell parallelepipeds.
77    Note that rotation $\theta$, initially in the $xz$ plane, is carried out first, then
78    rotation $\phi$ about the $z$ axis, finally rotation $\Psi$ is now around the axis of the cylinder.
79    The neutron or X-ray beam is along the $z$ axis.
80
81.. figure:: img/parallelepiped_angle_projection.png
82
83    Examples of the angles for oriented rectangular prisms against the
84    detector plane.
85
86
87
88Validation
89----------
90
91Validation of the code was conducted by comparing the output of the 1D model
92to the output of the existing :ref:`parallelepiped` model.
93
94
95References
96----------
97
98.. [#] P Mittelbach and G Porod, *Acta Physica Austriaca*, 14 (1961) 185-211
99.. [#] R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854
100.. [#] L. Onsager, *Ann. New York Acad. Sci.*, 51 (1949) 627-659
101
102Authorship and Verification
103----------------------------
104
105* **Author:**
106* **Last Modified by:**
107* **Last Reviewed by:**
108"""
109
110import numpy as np
111from numpy import inf
112
113name = "rectangular_prism"
114title = "Rectangular parallelepiped with uniform scattering length density."
115description = """
116    I(q)= scale*V*(sld - sld_solvent)^2*P(q,theta,phi)+background
117        P(q,theta,phi) = (2/pi) * double integral from 0 to pi/2 of ...
118           AP^2(q)*sin(theta)*dtheta*dphi
119        AP = S(q*C*cos(theta)/2) * S(q*A*sin(theta)*sin(phi)/2) * S(q*B*sin(theta)*cos(phi)/2)
120        S(x) = sin(x)/x
121"""
122category = "shape:parallelepiped"
123
124#             ["name", "units", default, [lower, upper], "type","description"],
125parameters = [["sld", "1e-6/Ang^2", 6.3, [-inf, inf], "sld",
126               "Parallelepiped scattering length density"],
127              ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "sld",
128               "Solvent scattering length density"],
129              ["length_a", "Ang", 35, [0, inf], "volume",
130               "Shorter side of the parallelepiped"],
131              ["b2a_ratio", "", 1, [0, inf], "volume",
132               "Ratio sides b/a"],
133              ["c2a_ratio", "", 1, [0, inf], "volume",
134               "Ratio sides c/a"],
135              ["theta", "degrees", 0, [-360, 360], "orientation",
136               "c axis to beam angle"],
137              ["phi", "degrees", 0, [-360, 360], "orientation",
138               "rotation about beam"],
139              ["psi", "degrees", 0, [-360, 360], "orientation",
140               "rotation about c axis"],
141             ]
142
143source = ["lib/gauss76.c", "rectangular_prism.c"]
144have_Fq = True
145radius_effective_modes = [
146    "equivalent cylinder excluded volume", "equivalent volume sphere",
147    "half length_a", "half length_b", "half length_c",
148    "equivalent circular cross-section", "half ab diagonal", "half diagonal",
149    ]
150
151def random():
152    """Return a random parameter set for the model."""
153    a, b, c = 10**np.random.uniform(1, 4.7, size=3)
154    pars = dict(
155        length_a=a,
156        b2a_ratio=b/a,
157        c2a_ratio=c/a,
158    )
159    return pars
160
161# parameters for demo
162demo = dict(scale=1, background=0,
163            sld=6.3, sld_solvent=1.0,
164            length_a=35, b2a_ratio=1, c2a_ratio=1,
165            length_a_pd=0.1, length_a_pd_n=10,
166            b2a_ratio_pd=0.1, b2a_ratio_pd_n=1,
167            c2a_ratio_pd=0.1, c2a_ratio_pd_n=1)
168
169tests = [[{}, 0.2, 0.375248406825],
170         [{}, [0.2], [0.375248406825]],
171        ]
Note: See TracBrowser for help on using the repository browser.