source: sasmodels/sasmodels/models/core_shell_parallelepiped.py @ dbf1a60

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since dbf1a60 was dbf1a60, checked in by butler, 5 years ago

Add comments to c code and clean documentatin

Added comments to c code in both parallelepiped and core shell
parallelepiped noting the change of integration varialbes in the
computation. Cleaned up and final corrections to the core shell
documentation and did some cleaning of parallelipiped. In particular
tried to bring a bit more consistency between the docs.

addresses #896

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