source: sasmodels/sasmodels/models/hollow_rectangular_prism.py @ 30b60d2

core_shell_microgelscostrafo411magnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 30b60d2 was 30b60d2, checked in by Paul Kienzle <pkienzle@…>, 5 years ago

allow build of both latex and html

  • 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 hollow rectangular
6parallelepiped with a wall of thickness $\Delta$.
7It computes only the 1D scattering, not the 2D.
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
15same length increment $2\Delta$ (Nayuk, 2012).
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
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
29
30.. math::
31  :nowrap:
32
33  \begin{align*}
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]
49  \end{align*}
50
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
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::
60  I(q) = \text{scale} \times V \times (\rho_\text{p} -
61  \rho_\text{solvent})^2 \times P(q) + \text{background}
62
63where $\rho_\text{p}$ is the scattering length of the parallelepiped,
64$\rho_\text{solvent}$ is the scattering length of the solvent,
65and (if the data are in absolute units) *scale* represents the volume fraction
66(which is unitless).
67
68**The 2D scattering intensity is not computed by this model.**
69
70
71Validation
72----------
73
74Validation of the code was conducted by qualitatively comparing the output
75of the 1D model to the curves shown in (Nayuk, 2012).
76
77
78References
79----------
80
81R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854
82
83"""
84
85from numpy import pi, inf, sqrt
86
87name = "hollow_rectangular_prism"
88title = "Hollow rectangular parallelepiped with uniform scattering length density."
89description = """
90    I(q)= scale*V*(sld - sld_solvent)^2*P(q,theta,phi)+background
91        P(q,theta,phi) = (2/pi/V^2) * double integral from 0 to pi/2 of ...
92           (AP1-AP2)^2(q)*sin(theta)*dtheta*dphi
93        AP1 = S(q*C*cos(theta)/2) * S(q*A*sin(theta)*sin(phi)/2) * S(q*B*sin(theta)*cos(phi)/2)
94        AP2 = S(q*C'*cos(theta)) * S(q*A'*sin(theta)*sin(phi)) * S(q*B'*sin(theta)*cos(phi))
95        C' = (C/2-thickness)
96        B' = (B/2-thickness)
97        A' = (A/2-thickness)
98        S(x) = sin(x)/x
99"""
100category = "shape:parallelepiped"
101
102#             ["name", "units", default, [lower, upper], "type","description"],
103parameters = [["sld", "1e-6/Ang^2", 6.3, [-inf, inf], "sld",
104               "Parallelepiped scattering length density"],
105              ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "sld",
106               "Solvent scattering length density"],
107              ["length_a", "Ang", 35, [0, inf], "volume",
108               "Shorter side of the parallelepiped"],
109              ["b2a_ratio", "Ang", 1, [0, inf], "volume",
110               "Ratio sides b/a"],
111              ["c2a_ratio", "Ang", 1, [0, inf], "volume",
112               "Ratio sides c/a"],
113              ["thickness", "Ang", 1, [0, inf], "volume",
114               "Thickness of parallelepiped"],
115             ]
116
117source = ["lib/gauss76.c", "hollow_rectangular_prism.c"]
118
119def ER(length_a, b2a_ratio, c2a_ratio, thickness):
120    """
121    Return equivalent radius (ER)
122    thickness parameter not used
123    """
124    b_side = length_a * b2a_ratio
125    c_side = length_a * c2a_ratio
126
127    # surface average radius (rough approximation)
128    surf_rad = sqrt(length_a * b_side / pi)
129
130    ddd = 0.75 * surf_rad * (2 * surf_rad * c_side + (c_side + surf_rad) * (c_side + pi * surf_rad))
131    return 0.5 * (ddd) ** (1. / 3.)
132
133def VR(length_a, b2a_ratio, c2a_ratio, thickness):
134    """
135    Return shell volume and total volume
136    """
137    b_side = length_a * b2a_ratio
138    c_side = length_a * c2a_ratio
139    a_core = length_a - 2.0*thickness
140    b_core = b_side - 2.0*thickness
141    c_core = c_side - 2.0*thickness
142    vol_core = a_core * b_core * c_core
143    vol_total = length_a * b_side * c_side
144    vol_shell = vol_total - vol_core
145    return vol_total, vol_shell
146
147
148def random():
149    import numpy as np
150    a, b, c = 10**np.random.uniform(1, 4.7, size=3)
151    # Thickness is limited to 1/2 the smallest dimension
152    # Use a distribution with a preference for thin shell or thin core
153    # Avoid core,shell radii < 1
154    min_dim = 0.5*min(a, b, c)
155    thickness = np.random.beta(0.5, 0.5)*(min_dim-2) + 1
156    #print(a, b, c, thickness, thickness/min_dim)
157    pars = dict(
158        length_a=a,
159        b2a_ratio=b/a,
160        c2a_ratio=c/a,
161        thickness=thickness,
162    )
163    return pars
164
165
166# parameters for demo
167demo = dict(scale=1, background=0,
168            sld=6.3, sld_solvent=1.0,
169            length_a=35, b2a_ratio=1, c2a_ratio=1, thickness=1,
170            length_a_pd=0.1, length_a_pd_n=10,
171            b2a_ratio_pd=0.1, b2a_ratio_pd_n=1,
172            c2a_ratio_pd=0.1, c2a_ratio_pd_n=1)
173
174tests = [[{}, 0.2, 0.76687283098],
175         [{}, [0.2], [0.76687283098]],
176        ]
177
Note: See TracBrowser for help on using the repository browser.