source: sasmodels/sasmodels/models/core_shell_parallelepiped.py @ 97be877

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 97be877 was 97be877, checked in by Paul Kienzle <pkienzle@…>, 6 years ago

update authorship/verification for core-shell parallelepiped

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