source: sasmodels/sasmodels/models/rectangular_prism.py @ 0507e09

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 0507e09 was 0507e09, checked in by smk78, 7 months ago

Added link to source code to each model. Closes #883

  • Property mode set to 100644
File size: 6.7 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
102Source
103------
104
105`rectangular_prism.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/rectangular_prism.py>`_
106
107`rectangular_prism.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/rectangular_prism.c>`_
108
109Authorship and Verification
110----------------------------
111
112* **Author:**
113* **Last Modified by:**
114* **Last Reviewed by:**
115* **Source added by :** Steve King **Date:** March 25, 2019
116"""
117
118import numpy as np
119from numpy import inf
120
121name = "rectangular_prism"
122title = "Rectangular parallelepiped with uniform scattering length density."
123description = """
124    I(q)= scale*V*(sld - sld_solvent)^2*P(q,theta,phi)+background
125        P(q,theta,phi) = (2/pi) * double integral from 0 to pi/2 of ...
126           AP^2(q)*sin(theta)*dtheta*dphi
127        AP = S(q*C*cos(theta)/2) * S(q*A*sin(theta)*sin(phi)/2) * S(q*B*sin(theta)*cos(phi)/2)
128        S(x) = sin(x)/x
129"""
130category = "shape:parallelepiped"
131
132#             ["name", "units", default, [lower, upper], "type","description"],
133parameters = [["sld", "1e-6/Ang^2", 6.3, [-inf, inf], "sld",
134               "Parallelepiped scattering length density"],
135              ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "sld",
136               "Solvent scattering length density"],
137              ["length_a", "Ang", 35, [0, inf], "volume",
138               "Shorter side of the parallelepiped"],
139              ["b2a_ratio", "", 1, [0, inf], "volume",
140               "Ratio sides b/a"],
141              ["c2a_ratio", "", 1, [0, inf], "volume",
142               "Ratio sides c/a"],
143              ["theta", "degrees", 0, [-360, 360], "orientation",
144               "c axis to beam angle"],
145              ["phi", "degrees", 0, [-360, 360], "orientation",
146               "rotation about beam"],
147              ["psi", "degrees", 0, [-360, 360], "orientation",
148               "rotation about c axis"],
149             ]
150
151source = ["lib/gauss76.c", "rectangular_prism.c"]
152have_Fq = True
153effective_radius_type = [
154    "equivalent cylinder excluded volume", "equivalent volume sphere",
155    "half length_a", "half length_b", "half length_c",
156    "equivalent circular cross-section", "half ab diagonal", "half diagonal",
157    ]
158
159def random():
160    """Return a random parameter set for the model."""
161    a, b, c = 10**np.random.uniform(1, 4.7, size=3)
162    pars = dict(
163        length_a=a,
164        b2a_ratio=b/a,
165        c2a_ratio=c/a,
166    )
167    return pars
168
169# parameters for demo
170demo = dict(scale=1, background=0,
171            sld=6.3, sld_solvent=1.0,
172            length_a=35, b2a_ratio=1, c2a_ratio=1,
173            length_a_pd=0.1, length_a_pd_n=10,
174            b2a_ratio_pd=0.1, b2a_ratio_pd_n=1,
175            c2a_ratio_pd=0.1, c2a_ratio_pd_n=1)
176
177tests = [[{}, 0.2, 0.375248406825],
178         [{}, [0.2], [0.375248406825]],
179        ]
Note: See TracBrowser for help on using the repository browser.