source: sasmodels/sasmodels/models/hollow_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, 6 years ago

update documentation for core_shell_parallelepiped, hollow_rectangular_prism and rectangular_prism

  • Property mode set to 100644
File size: 7.6 KB
RevLine 
[deb7ee0]1# rectangular_prism model
2# Note: model title and parameter table are inserted automatically
3r"""
4
[117090a]5This model provides the form factor, $P(q)$, for a hollow rectangular
6parallelepiped with a wall of thickness $\Delta$.
[393facf]7
[deb7ee0]8
9Definition
10----------
11
12The 1D scattering intensity for this model is calculated by forming
13the difference of the amplitudes of two massive parallelepipeds
14differing in their outermost dimensions in each direction by the
[ab2aea8]15same length increment $2\Delta$ (Nayuk, 2012).
[deb7ee0]16
17As in the case of the massive parallelepiped model (:ref:`rectangular-prism`),
18the scattering amplitude is computed for a particular orientation of the
19parallelepiped with respect to the scattering vector and then averaged over all
20possible orientations, giving
21
22.. math::
23  P(q) =  \frac{1}{V^2} \frac{2}{\pi} \times \, \int_0^{\frac{\pi}{2}} \,
24  \int_0^{\frac{\pi}{2}} A_{P\Delta}^2(q) \, \sin\theta \, d\theta \, d\phi
25
[117090a]26where $\theta$ is the angle between the $z$ axis and the longest axis
27of the parallelepiped, $\phi$ is the angle between the scattering vector
28(lying in the $xy$ plane) and the $y$ axis, and
[deb7ee0]29
30.. math::
[5111921]31  :nowrap:
32
[30b60d2]33  \begin{align*}
[deb7ee0]34  A_{P\Delta}(q) & =  A B C
35    \left[\frac{\sin \bigl( q \frac{C}{2} \cos\theta \bigr)}
36    {\left( q \frac{C}{2} \cos\theta \right)} \right]
37    \left[\frac{\sin \bigl( q \frac{A}{2} \sin\theta \sin\phi \bigr)}
38    {\left( q \frac{A}{2} \sin\theta \sin\phi \right)}\right]
39    \left[\frac{\sin \bigl( q \frac{B}{2} \sin\theta \cos\phi \bigr)}
40    {\left( q \frac{B}{2} \sin\theta \cos\phi \right)}\right] \\
41    & - 8
42    \left(\frac{A}{2}-\Delta\right) \left(\frac{B}{2}-\Delta\right) \left(\frac{C}{2}-\Delta\right)
43    \left[ \frac{\sin \bigl[ q \bigl(\frac{C}{2}-\Delta\bigr) \cos\theta \bigr]}
44    {q \bigl(\frac{C}{2}-\Delta\bigr) \cos\theta} \right]
45    \left[ \frac{\sin \bigl[ q \bigl(\frac{A}{2}-\Delta\bigr) \sin\theta \sin\phi \bigr]}
46    {q \bigl(\frac{A}{2}-\Delta\bigr) \sin\theta \sin\phi} \right]
47    \left[ \frac{\sin \bigl[ q \bigl(\frac{B}{2}-\Delta\bigr) \sin\theta \cos\phi \bigr]}
48    {q \bigl(\frac{B}{2}-\Delta\bigr) \sin\theta \cos\phi} \right]
[30b60d2]49  \end{align*}
[deb7ee0]50
[117090a]51where $A$, $B$ and $C$ are the external sides of the parallelepiped fulfilling
52$A \le B \le C$, and the volume $V$ of the parallelepiped is
[deb7ee0]53
54.. math::
55  V = A B C \, - \, (A - 2\Delta) (B - 2\Delta) (C - 2\Delta)
56
57The 1D scattering intensity is then calculated as
58
59.. math::
[ab2aea8]60  I(q) = \text{scale} \times V \times (\rho_\text{p} -
61  \rho_\text{solvent})^2 \times P(q) + \text{background}
[deb7ee0]62
[ab2aea8]63where $\rho_\text{p}$ is the scattering length of the parallelepiped,
64$\rho_\text{solvent}$ is the scattering length of the solvent,
[deb7ee0]65and (if the data are in absolute units) *scale* represents the volume fraction
66(which is unitless).
67
[393facf]68For 2d data the orientation of the particle is required, described using
69angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details
70of the calculation and angular dispersions see :ref:`orientation` .
71The angle $\Psi$ is the rotational angle around the long *C* axis. For example,
72$\Psi = 0$ when the *B* axis is parallel to the *x*-axis of the detector.
73
74For 2d, constraints must be applied during fitting to ensure that the inequality
75$A < B < C$ is not violated, and hence the correct definition of angles is preserved. The calculation will not report an error,
76but the results may be not correct.
77
78.. figure:: img/parallelepiped_angle_definition.png
79
80    Definition of the angles for oriented hollow rectangular prism.
81    Note that rotation $\theta$, initially in the $xz$ plane, is carried out first, then
82    rotation $\phi$ about the $z$ axis, finally rotation $\Psi$ is now around the axis of the prism.
83    The neutron or X-ray beam is along the $z$ axis.
84
85.. figure:: img/parallelepiped_angle_projection.png
86
87    Examples of the angles for oriented hollow rectangular prisms against the
88    detector plane.
[deb7ee0]89
90
91Validation
92----------
93
94Validation of the code was conducted by qualitatively comparing the output
95of the 1D model to the curves shown in (Nayuk, 2012).
96
[aa2edb2]97
98References
99----------
[deb7ee0]100
101R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854
102
103"""
104
105from numpy import pi, inf, sqrt
106
107name = "hollow_rectangular_prism"
108title = "Hollow rectangular parallelepiped with uniform scattering length density."
109description = """
[3d8283b]110    I(q)= scale*V*(sld - sld_solvent)^2*P(q,theta,phi)+background
[deb7ee0]111        P(q,theta,phi) = (2/pi/V^2) * double integral from 0 to pi/2 of ...
112           (AP1-AP2)^2(q)*sin(theta)*dtheta*dphi
113        AP1 = S(q*C*cos(theta)/2) * S(q*A*sin(theta)*sin(phi)/2) * S(q*B*sin(theta)*cos(phi)/2)
114        AP2 = S(q*C'*cos(theta)) * S(q*A'*sin(theta)*sin(phi)) * S(q*B'*sin(theta)*cos(phi))
115        C' = (C/2-thickness)
116        B' = (B/2-thickness)
117        A' = (A/2-thickness)
118        S(x) = sin(x)/x
119"""
120category = "shape:parallelepiped"
121
122#             ["name", "units", default, [lower, upper], "type","description"],
[42356c8]123parameters = [["sld", "1e-6/Ang^2", 6.3, [-inf, inf], "sld",
[deb7ee0]124               "Parallelepiped scattering length density"],
[42356c8]125              ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "sld",
[deb7ee0]126               "Solvent scattering length density"],
[a807206]127              ["length_a", "Ang", 35, [0, inf], "volume",
[deb7ee0]128               "Shorter side of the parallelepiped"],
129              ["b2a_ratio", "Ang", 1, [0, inf], "volume",
130               "Ratio sides b/a"],
131              ["c2a_ratio", "Ang", 1, [0, inf], "volume",
132               "Ratio sides c/a"],
133              ["thickness", "Ang", 1, [0, inf], "volume",
134               "Thickness of parallelepiped"],
[8de1477]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"],
[deb7ee0]141             ]
142
[40a87fa]143source = ["lib/gauss76.c", "hollow_rectangular_prism.c"]
[deb7ee0]144
[a807206]145def ER(length_a, b2a_ratio, c2a_ratio, thickness):
[deb7ee0]146    """
[40a87fa]147    Return equivalent radius (ER)
148    thickness parameter not used
[deb7ee0]149    """
[a807206]150    b_side = length_a * b2a_ratio
151    c_side = length_a * c2a_ratio
[deb7ee0]152
153    # surface average radius (rough approximation)
[a807206]154    surf_rad = sqrt(length_a * b_side / pi)
[deb7ee0]155
156    ddd = 0.75 * surf_rad * (2 * surf_rad * c_side + (c_side + surf_rad) * (c_side + pi * surf_rad))
157    return 0.5 * (ddd) ** (1. / 3.)
158
[a807206]159def VR(length_a, b2a_ratio, c2a_ratio, thickness):
[deb7ee0]160    """
[40a87fa]161    Return shell volume and total volume
[deb7ee0]162    """
[a807206]163    b_side = length_a * b2a_ratio
164    c_side = length_a * c2a_ratio
165    a_core = length_a - 2.0*thickness
[deb7ee0]166    b_core = b_side - 2.0*thickness
167    c_core = c_side - 2.0*thickness
168    vol_core = a_core * b_core * c_core
[a807206]169    vol_total = length_a * b_side * c_side
[deb7ee0]170    vol_shell = vol_total - vol_core
171    return vol_total, vol_shell
172
173
[31df0c9]174def random():
175    import numpy as np
176    a, b, c = 10**np.random.uniform(1, 4.7, size=3)
[8f04da4]177    # Thickness is limited to 1/2 the smallest dimension
178    # Use a distribution with a preference for thin shell or thin core
179    # Avoid core,shell radii < 1
180    min_dim = 0.5*min(a, b, c)
181    thickness = np.random.beta(0.5, 0.5)*(min_dim-2) + 1
182    #print(a, b, c, thickness, thickness/min_dim)
[31df0c9]183    pars = dict(
184        length_a=a,
185        b2a_ratio=b/a,
186        c2a_ratio=c/a,
187        thickness=thickness,
188    )
189    return pars
190
191
[deb7ee0]192# parameters for demo
193demo = dict(scale=1, background=0,
[ab2aea8]194            sld=6.3, sld_solvent=1.0,
[a807206]195            length_a=35, b2a_ratio=1, c2a_ratio=1, thickness=1,
196            length_a_pd=0.1, length_a_pd_n=10,
[deb7ee0]197            b2a_ratio_pd=0.1, b2a_ratio_pd_n=1,
198            c2a_ratio_pd=0.1, c2a_ratio_pd_n=1)
199
[6dd90c1]200tests = [[{}, 0.2, 0.76687283098],
201         [{}, [0.2], [0.76687283098]],
[deb7ee0]202        ]
Note: See TracBrowser for help on using the repository browser.