source: sasmodels/sasmodels/models/rectangular_prism.py @ 99658f6

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 99658f6 was 99658f6, checked in by grethevj, 13 months ago

updated ER functions including cylinder excluded volume, to match 4.x

  • Property mode set to 100644
File size: 6.2 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
102L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).
103
104"""
105
106import numpy as np
107from numpy import pi, inf, sqrt
108
109name = "rectangular_prism"
110title = "Rectangular parallelepiped with uniform scattering length density."
111description = """
112    I(q)= scale*V*(sld - sld_solvent)^2*P(q,theta,phi)+background
113        P(q,theta,phi) = (2/pi) * double integral from 0 to pi/2 of ...
114           AP^2(q)*sin(theta)*dtheta*dphi
115        AP = S(q*C*cos(theta)/2) * S(q*A*sin(theta)*sin(phi)/2) * S(q*B*sin(theta)*cos(phi)/2)
116        S(x) = sin(x)/x
117"""
118category = "shape:parallelepiped"
119
120#             ["name", "units", default, [lower, upper], "type","description"],
121parameters = [["sld", "1e-6/Ang^2", 6.3, [-inf, inf], "sld",
122               "Parallelepiped scattering length density"],
123              ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "sld",
124               "Solvent scattering length density"],
125              ["length_a", "Ang", 35, [0, inf], "volume",
126               "Shorter side of the parallelepiped"],
127              ["b2a_ratio", "", 1, [0, inf], "volume",
128               "Ratio sides b/a"],
129              ["c2a_ratio", "", 1, [0, inf], "volume",
130               "Ratio sides c/a"],
131              ["theta", "degrees", 0, [-360, 360], "orientation",
132               "c axis to beam angle"],
133              ["phi", "degrees", 0, [-360, 360], "orientation",
134               "rotation about beam"],
135              ["psi", "degrees", 0, [-360, 360], "orientation",
136               "rotation about c axis"],
137             ]
138
139source = ["lib/gauss76.c", "rectangular_prism.c"]
140have_Fq = True
141effective_radius_type = [
142    "equivalent cylinder excluded volume", "equivalent volume sphere", 
143    "half length_a", "half length_b", "half length_c",
144    "equivalent circular cross-section", "half ab diagonal", "half diagonal",
145    ]
146
147def random():
148    a, b, c = 10**np.random.uniform(1, 4.7, size=3)
149    pars = dict(
150        length_a=a,
151        b2a_ratio=b/a,
152        c2a_ratio=c/a,
153    )
154    return pars
155
156# parameters for demo
157demo = dict(scale=1, background=0,
158            sld=6.3, sld_solvent=1.0,
159            length_a=35, b2a_ratio=1, c2a_ratio=1,
160            length_a_pd=0.1, length_a_pd_n=10,
161            b2a_ratio_pd=0.1, b2a_ratio_pd_n=1,
162            c2a_ratio_pd=0.1, c2a_ratio_pd_n=1)
163
164tests = [[{}, 0.2, 0.375248406825],
165         [{}, [0.2], [0.375248406825]],
166        ]
Note: See TracBrowser for help on using the repository browser.