source: sasmodels/sasmodels/models/core_shell_parallelepiped.py @ 367886f

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 367886f was 367886f, checked in by butler, 6 years ago

More edits

Add angle variables, Put in terms of P(q), and add the double integral
equation over which the F(q) given will be integrated.

  • Property mode set to 100644
File size: 9.4 KB
Line 
1r"""
2Definition
3----------
4
5Calculates the form factor for a rectangular solid with a core-shell structure.
6The thickness and the scattering length density of the shell or
7"rim" can be different on each (pair) of faces.
8
9The form factor is normalized by the particle volume $V$ such that
10
11.. math::
12
13    I(q) = \text{scale}\frac{\langle P(q,\alpha,\beta) \rangle}{V}
14    + \text{background}
15
16where $\langle \ldots \rangle$ is an average over all possible orientations
17of the rectangular solid.
18
19The function calculated is the form factor of the rectangular solid below.
20The core of the solid is defined by the dimensions $A$, $B$, $C$ such that
21$A < B < C$.
22
23.. figure:: img/parallelepiped_geometry.jpg
24
25   Core of the core shell Parallelepiped with the corresponding definition
26   of sides.
27
28
29There are rectangular "slabs" of thickness $t_A$ that add to the $A$ dimension
30(on the $BC$ faces). There are similar slabs on the $AC$ $(=t_B)$ and $AB$
31$(=t_C)$ faces. The projection in the $AB$ plane is then
32
33.. figure:: img/core_shell_parallelepiped_projection.jpg
34
35   AB cut through the core-shell parllelipiped showing the cross secion of
36   four of the six shell slabs
37
38The volume of the solid is
39
40.. math::
41
42    V = ABC + 2t_ABC + 2t_BAC + 2t_CAB
43
44**meaning that there are "gaps" at the corners of the solid.**
45
46The intensity calculated follows the :ref:`parallelepiped` model, with the
47core-shell intensity being calculated as the square of the sum of the
48amplitudes of the core and the slabs on the edges.
49
50the scattering amplitude is computed for a particular orientation of the
51core-shell parallelepiped with respect to the scattering vector and then
52averaged over all possible orientations, where $\alpha$ is the angle between
53the $z$ axis and the $C$ axis of the parallelepiped, $\beta$ is
54the angle between projection of the particle in the $xy$ detector plane
55and the $y$ axis.
56
57.. math::
58
59    P(q)=\int_{0}^{\pi/2}\int_{0}^{\pi/2}F^2(q,\alpha,\beta) \ cos\alpha
60    \ d\alpha \ d\beta
61
62and
63
64.. math::
65
66    F(q)
67    &= (\rho_\text{core}-\rho_\text{solvent})
68       S(Q_A, A) S(Q_B, B) S(Q_C, C) \\
69    &+ (\rho_\text{A}-\rho_\text{solvent})
70        \left[S(Q_A, A+2t_A) - S(Q_A, A)\right] S(Q_B, B) S(Q_C, C) \\
71    &+ (\rho_\text{B}-\rho_\text{solvent})
72        S(Q_A, A) \left[S(Q_B, B+2t_B) - S(Q_B, B)\right] S(Q_C, C) \\
73    &+ (\rho_\text{C}-\rho_\text{solvent})
74        S(Q_A, A) S(Q_B, B) \left[S(Q_C, C+2t_C) - S(Q_C, C)\right]
75
76with
77
78.. math::
79
80    S(Q, L) = L \frac{\sin \tfrac{1}{2} Q L}{\tfrac{1}{2} Q L}
81
82and
83
84.. math::
85
86    Q_A &= q \sin\alpha \sin\beta \\
87    Q_B &= q \sin\alpha \cos\beta \\
88    Q_C &= q \cos\alpha
89
90
91where $\rho_\text{core}$, $\rho_\text{A}$, $\rho_\text{B}$ and $\rho_\text{C}$
92are the scattering length of the parallelepiped core, and the rectangular
93slabs of thickness $t_A$, $t_B$ and $t_C$, respectively. $\rho_\text{solvent}$
94is the scattering length of the solvent.
95
96FITTING NOTES
97~~~~~~~~~~~~~
98
99If the scale is set equal to the particle volume fraction, $\phi$, the returned
100value is the scattered intensity per unit volume, $I(q) = \phi P(q)$. However,
101**no interparticle interference effects are included in this calculation.**
102
103There are many parameters in this model. Hold as many fixed as possible with
104known values, or you will certainly end up at a solution that is unphysical.
105
106The returned value is in units of |cm^-1|, on absolute scale.
107
108NB: The 2nd virial coefficient of the core_shell_parallelepiped is calculated
109based on the the averaged effective radius $(=\sqrt{(A+2t_A)(B+2t_B)/\pi})$
110and length $(C+2t_C)$ values, after appropriately sorting the three dimensions
111to give an oblate or prolate particle, to give an effective radius,
112for $S(q)$ when $P(q) * S(q)$ is applied.
113
114For 2d data the orientation of the particle is required, described using
115angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further
116details of the calculation and angular dispersions see :ref:`orientation`.
117The angle $\Psi$ is the rotational angle around the *long_c* axis. For example,
118$\Psi = 0$ when the *short_b* axis is parallel to the *x*-axis of the detector.
119
120For 2d, constraints must be applied during fitting to ensure that the
121inequality $A < B < C$ is not violated, and hence the correct definition
122of angles is preserved. The calculation will not report an error,
123but the results may be not correct.
124
125.. figure:: img/parallelepiped_angle_definition.png
126
127    Definition of the angles for oriented core-shell parallelepipeds.
128    Note that rotation $\theta$, initially in the $xz$ plane, is carried
129    out first, then rotation $\phi$ about the $z$ axis, finally rotation
130    $\Psi$ is now around the axis of the cylinder. The neutron or X-ray
131    beam is along the $z$ axis.
132
133.. figure:: img/parallelepiped_angle_projection.png
134
135    Examples of the angles for oriented core-shell parallelepipeds against the
136    detector plane.
137
138References
139----------
140
141.. [#] P Mittelbach and G Porod, *Acta Physica Austriaca*, 14 (1961) 185-211
142    Equations (1), (13-14). (in German)
143.. [#] D Singh (2009). *Small angle scattering studies of self assembly in
144   lipid mixtures*, Johns Hopkins University Thesis (2009) 223-225. `Available
145   from Proquest <http://search.proquest.com/docview/304915826?accountid
146   =26379>`_
147
148Authorship and Verification
149----------------------------
150
151* **Author:** NIST IGOR/DANSE **Date:** pre 2010
152* **Converted to sasmodels by:** Miguel Gonzales **Date:** February 26, 2016
153* **Last Modified by:** Paul Kienzle **Date:** October 17, 2017
154* Cross-checked against hollow rectangular prism and rectangular prism for
155  equal thickness overlapping sides, and by Monte Carlo sampling of points
156  within the shape for non-uniform, non-overlapping sides.
157"""
158
159import numpy as np
160from numpy import pi, inf, sqrt, cos, sin
161
162name = "core_shell_parallelepiped"
163title = "Rectangular solid with a core-shell structure."
164description = """
165     P(q)=
166"""
167category = "shape:parallelepiped"
168
169#             ["name", "units", default, [lower, upper], "type","description"],
170parameters = [["sld_core", "1e-6/Ang^2", 1, [-inf, inf], "sld",
171               "Parallelepiped core scattering length density"],
172              ["sld_a", "1e-6/Ang^2", 2, [-inf, inf], "sld",
173               "Parallelepiped A rim scattering length density"],
174              ["sld_b", "1e-6/Ang^2", 4, [-inf, inf], "sld",
175               "Parallelepiped B rim scattering length density"],
176              ["sld_c", "1e-6/Ang^2", 2, [-inf, inf], "sld",
177               "Parallelepiped C rim scattering length density"],
178              ["sld_solvent", "1e-6/Ang^2", 6, [-inf, inf], "sld",
179               "Solvent scattering length density"],
180              ["length_a", "Ang", 35, [0, inf], "volume",
181               "Shorter side of the parallelepiped"],
182              ["length_b", "Ang", 75, [0, inf], "volume",
183               "Second side of the parallelepiped"],
184              ["length_c", "Ang", 400, [0, inf], "volume",
185               "Larger side of the parallelepiped"],
186              ["thick_rim_a", "Ang", 10, [0, inf], "volume",
187               "Thickness of A rim"],
188              ["thick_rim_b", "Ang", 10, [0, inf], "volume",
189               "Thickness of B rim"],
190              ["thick_rim_c", "Ang", 10, [0, inf], "volume",
191               "Thickness of C rim"],
192              ["theta", "degrees", 0, [-360, 360], "orientation",
193               "c axis to beam angle"],
194              ["phi", "degrees", 0, [-360, 360], "orientation",
195               "rotation about beam"],
196              ["psi", "degrees", 0, [-360, 360], "orientation",
197               "rotation about c axis"],
198             ]
199
200source = ["lib/gauss76.c", "core_shell_parallelepiped.c"]
201
202
203def ER(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c):
204    """
205        Return equivalent radius (ER)
206    """
207    from .parallelepiped import ER as ER_p
208
209    a = length_a + 2*thick_rim_a
210    b = length_b + 2*thick_rim_b
211    c = length_c + 2*thick_rim_c
212    return ER_p(a, b, c)
213
214# VR defaults to 1.0
215
216def random():
217    outer = 10**np.random.uniform(1, 4.7, size=3)
218    thick = np.random.beta(0.5, 0.5, size=3)*(outer-2) + 1
219    length = outer - thick
220    pars = dict(
221        length_a=length[0],
222        length_b=length[1],
223        length_c=length[2],
224        thick_rim_a=thick[0],
225        thick_rim_b=thick[1],
226        thick_rim_c=thick[2],
227    )
228    return pars
229
230# parameters for demo
231demo = dict(scale=1, background=0.0,
232            sld_core=1, sld_a=2, sld_b=4, sld_c=2, sld_solvent=6,
233            length_a=35, length_b=75, length_c=400,
234            thick_rim_a=10, thick_rim_b=10, thick_rim_c=10,
235            theta=0, phi=0, psi=0,
236            length_a_pd=0.1, length_a_pd_n=1,
237            length_b_pd=0.1, length_b_pd_n=1,
238            length_c_pd=0.1, length_c_pd_n=1,
239            thick_rim_a_pd=0.1, thick_rim_a_pd_n=1,
240            thick_rim_b_pd=0.1, thick_rim_b_pd_n=1,
241            thick_rim_c_pd=0.1, thick_rim_c_pd_n=1,
242            theta_pd=10, theta_pd_n=1,
243            phi_pd=10, phi_pd_n=1,
244            psi_pd=10, psi_pd_n=1)
245
246# rkh 7/4/17 add random unit test for 2d, note make all params different,
247# 2d values not tested against other codes or models
248if 0:  # pak: model rewrite; need to update tests
249    qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.)
250    tests = [[{}, 0.2, 0.533149288477],
251             [{}, [0.2], [0.533149288477]],
252             [{'theta':10.0, 'phi':20.0}, (qx, qy), 0.0853299803222],
253             [{'theta':10.0, 'phi':20.0}, [(qx, qy)], [0.0853299803222]],
254            ]
255    del qx, qy  # not necessary to delete, but cleaner
Note: See TracBrowser for help on using the repository browser.