source: sasmodels/sasmodels/models/core_shell_parallelepiped.py @ 1f159bd

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

fcc/bcc paracrystal, core shell parallelepiped: suppress broken tests

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