source: sasmodels/sasmodels/models/rectangular_prism.py @ 393facf

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 393facf was 393facf, checked in by dirk, 5 years ago

update documentation for core_shell_parallelepiped, hollow_rectangular_prism and rectangular_prism

  • Property mode set to 100644
File size: 6.3 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
98P Mittelbach and G Porod, *Acta Physica Austriaca*, 14 (1961) 185-211
99
100R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854
101
102"""
103
104from numpy import pi, inf, sqrt
105
106name = "rectangular_prism"
107title = "Rectangular parallelepiped with uniform scattering length density."
108description = """
109    I(q)= scale*V*(sld - sld_solvent)^2*P(q,theta,phi)+background
110        P(q,theta,phi) = (2/pi) * double integral from 0 to pi/2 of ...
111           AP^2(q)*sin(theta)*dtheta*dphi
112        AP = S(q*C*cos(theta)/2) * S(q*A*sin(theta)*sin(phi)/2) * S(q*B*sin(theta)*cos(phi)/2)
113        S(x) = sin(x)/x
114"""
115category = "shape:parallelepiped"
116
117#             ["name", "units", default, [lower, upper], "type","description"],
118parameters = [["sld", "1e-6/Ang^2", 6.3, [-inf, inf], "sld",
119               "Parallelepiped scattering length density"],
120              ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "sld",
121               "Solvent scattering length density"],
122              ["length_a", "Ang", 35, [0, inf], "volume",
123               "Shorter side of the parallelepiped"],
124              ["b2a_ratio", "", 1, [0, inf], "volume",
125               "Ratio sides b/a"],
126              ["c2a_ratio", "", 1, [0, inf], "volume",
127               "Ratio sides c/a"],
128              ["theta", "degrees", 0, [-360, 360], "orientation",
129               "c axis to beam angle"],
130              ["phi", "degrees", 0, [-360, 360], "orientation",
131               "rotation about beam"],
132              ["psi", "degrees", 0, [-360, 360], "orientation",
133               "rotation about c axis"],
134             ]
135
136source = ["lib/gauss76.c", "rectangular_prism.c"]
137
138def ER(length_a, b2a_ratio, c2a_ratio):
139    """
140        Return equivalent radius (ER)
141    """
142    b_side = length_a * b2a_ratio
143    c_side = length_a * c2a_ratio
144
145    # surface average radius (rough approximation)
146    surf_rad = sqrt(length_a * b_side / pi)
147
148    ddd = 0.75 * surf_rad * (2 * surf_rad * c_side + (c_side + surf_rad) * (c_side + pi * surf_rad))
149    return 0.5 * (ddd) ** (1. / 3.)
150
151def random():
152    import numpy as np
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.