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

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

Merge branch 'ticket-776-orientation' into ticket-786

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