source: sasmodels/sasmodels/models/core_shell_parallelepiped.py @ 0507e09

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 0507e09 was 0507e09, checked in by smk78, 3 months ago

Added link to source code to each model. Closes #883

  • Property mode set to 100644
File size: 11.5 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 "rim" can be
7different on each (pair) of faces. The three dimensions of the core of the
8parallelepiped (strictly here a cuboid) may be given in *any* size order as
9long as the particles are randomly oriented (i.e. take on all possible
10orientations see notes on 2D below). To avoid multiple fit solutions,
11especially with Monte-Carlo fit methods, it may be advisable to restrict their
12ranges. There may be a number of closely similar "best fits", so some trial and
13error, or fixing of some dimensions at expected values, may help.
14
15The form factor is normalized by the particle volume $V$ such that
16
17.. math::
18
19    I(q) = \frac{\text{scale}}{V} \langle P(q,\alpha,\beta) \rangle
20    + \text{background}
21
22where $\langle \ldots \rangle$ is an average over all possible orientations
23of the rectangular solid, and the usual $\Delta \rho^2 \ V^2$ term cannot be
24pulled out of the form factor term due to the multiple slds in the model.
25
26The core of the solid is defined by the dimensions $A$, $B$, $C$ here shown
27such that $A < B < C$.
28
29.. figure:: img/parallelepiped_geometry.jpg
30
31   Core of the core shell parallelepiped with the corresponding definition
32   of sides.
33
34
35There are rectangular "slabs" of thickness $t_A$ that add to the $A$ dimension
36(on the $BC$ faces). There are similar slabs on the $AC$ $(=t_B)$ and $AB$
37$(=t_C)$ faces. The projection in the $AB$ plane is
38
39.. figure:: img/core_shell_parallelepiped_projection.jpg
40
41   AB cut through the core-shell parallelipiped showing the cross secion of
42   four of the six shell slabs. As can be seen, this model leaves **"gaps"**
43   at the corners of the solid.
44
45
46The total volume of the solid is thus given as
47
48.. math::
49
50    V = ABC + 2t_ABC + 2t_BAC + 2t_CAB
51
52The intensity calculated follows the :ref:`parallelepiped` model, with the
53core-shell intensity being calculated as the square of the sum of the
54amplitudes of the core and the slabs on the edges. The scattering amplitude is
55computed for a particular orientation of the core-shell parallelepiped with
56respect to the scattering vector and then averaged over all possible
57orientations, where $\alpha$ is the angle between the $z$ axis and the $C$ axis
58of the parallelepiped, and $\beta$ is the angle between the projection of the
59particle in the $xy$ detector plane and the $y$ axis.
60
61.. math::
62
63    P(q)=\frac {\int_{0}^{\pi/2}\int_{0}^{\pi/2}F^2(q,\alpha,\beta) \ sin\alpha
64    \ d\alpha \ d\beta} {\int_{0}^{\pi/2} \ sin\alpha \ d\alpha \ d\beta}
65
66and
67
68.. math::
69
70    F(q,\alpha,\beta)
71    &= (\rho_\text{core}-\rho_\text{solvent})
72       S(Q_A, A) S(Q_B, B) S(Q_C, C) \\
73    &+ (\rho_\text{A}-\rho_\text{solvent})
74        \left[S(Q_A, A+2t_A) - S(Q_A, A)\right] S(Q_B, B) S(Q_C, C) \\
75    &+ (\rho_\text{B}-\rho_\text{solvent})
76        S(Q_A, A) \left[S(Q_B, B+2t_B) - S(Q_B, B)\right] S(Q_C, C) \\
77    &+ (\rho_\text{C}-\rho_\text{solvent})
78        S(Q_A, A) S(Q_B, B) \left[S(Q_C, C+2t_C) - S(Q_C, C)\right]
79
80with
81
82.. math::
83
84    S(Q_X, L) = L \frac{\sin (\tfrac{1}{2} Q_X L)}{\tfrac{1}{2} Q_X L}
85
86and
87
88.. math::
89
90    Q_A &= q \sin\alpha \sin\beta \\
91    Q_B &= q \sin\alpha \cos\beta \\
92    Q_C &= q \cos\alpha
93
94
95where $\rho_\text{core}$, $\rho_\text{A}$, $\rho_\text{B}$ and $\rho_\text{C}$
96are the scattering lengths of the parallelepiped core, and the rectangular
97slabs of thickness $t_A$, $t_B$ and $t_C$, respectively. $\rho_\text{solvent}$
98is the scattering length of the solvent.
99
100.. note::
101
102   the code actually implements two substitutions: $d(cos\alpha)$ is
103   substituted for -$sin\alpha \ d\alpha$ (note that in the
104   :ref:`parallelepiped` code this is explicitly implemented with
105   $\sigma = cos\alpha$), and $\beta$ is set to $\beta = u \pi/2$ so that
106   $du = \pi/2 \ d\beta$.  Thus both integrals go from 0 to 1 rather than 0
107   to $\pi/2$.
108
109FITTING NOTES
110~~~~~~~~~~~~~
111
112#. There are many parameters in this model. Hold as many fixed as possible with
113   known values, or you will certainly end up at a solution that is unphysical.
114
115#. The 2nd virial coefficient of the core_shell_parallelepiped is calculated
116   based on the the averaged effective radius $(=\sqrt{(A+2t_A)(B+2t_B)/\pi})$
117   and length $(C+2t_C)$ values, after appropriately sorting the three
118   dimensions to give an oblate or prolate particle, to give an effective radius
119   for $S(q)$ when $P(q) * S(q)$ is applied.
120
121#. For 2d data the orientation of the particle is required, described using
122   angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, where $\theta$
123   and $\phi$ define the orientation of the director in the laboratry reference
124   frame of the beam direction ($z$) and detector plane ($x-y$ plane), while
125   the angle $\Psi$ is effectively the rotational angle around the particle
126   $C$ axis. For $\theta = 0$ and $\phi = 0$, $\Psi = 0$ corresponds to the
127   $B$ axis oriented parallel to the y-axis of the detector with $A$ along
128   the x-axis. For other $\theta$, $\phi$ values, the order of rotations
129   matters. In particular, the parallelepiped must first be rotated $\theta$
130   degrees in the $x-z$ plane before rotating $\phi$ degrees around the $z$
131   axis (in the $x-y$ plane). Applying orientational distribution to the
132   particle orientation (i.e  `jitter` to one or more of these angles) can get
133   more confusing as `jitter` is defined **NOT** with respect to the laboratory
134   frame but the particle reference frame. It is thus highly recmmended to
135   read :ref:`orientation` for further details of the calculation and angular
136   dispersions.
137
138.. note:: For 2d, constraints must be applied during fitting to ensure that the
139   order of sides chosen is not altered, and hence that the correct definition
140   of angles is preserved. For the default choice shown here, that means
141   ensuring that the inequality $A < B < C$ is not violated,  The calculation
142   will not report an error, but the results may be not correct.
143
144.. figure:: img/parallelepiped_angle_definition.png
145
146    Definition of the angles for oriented core-shell parallelepipeds.
147    Note that rotation $\theta$, initially in the $x-z$ plane, is carried
148    out first, then rotation $\phi$ about the $z$ axis, finally rotation
149    $\Psi$ is now around the $C$ axis of the particle. The neutron or X-ray
150    beam is along the $z$ axis and the detecotr defines the $x-y$ plane.
151
152.. figure:: img/parallelepiped_angle_projection.png
153
154    Examples of the angles for oriented core-shell parallelepipeds against the
155    detector plane.
156
157
158Validation
159----------
160
161Cross-checked against hollow rectangular prism and rectangular prism for equal
162thickness overlapping sides, and by Monte Carlo sampling of points within the
163shape for non-uniform, non-overlapping sides.
164
165
166References
167----------
168
169.. [#] P Mittelbach and G Porod, *Acta Physica Austriaca*, 14 (1961) 185-211
170    Equations (1), (13-14). (in German)
171.. [#] D Singh (2009). *Small angle scattering studies of self assembly in
172   lipid mixtures*, Johns Hopkins University Thesis (2009) 223-225. `Available
173   from Proquest <http://search.proquest.com/docview/304915826?accountid
174   =26379>`_
175.. [#] L. Onsager, *Ann. New York Acad. Sci.*, 51 (1949) 627-659
176
177Source
178------
179
180`core_shell_parallelepiped.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/core_shell_parallelepiped.py>`_
181
182`core_shell_parallelepiped.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/core_shell_parallelepiped.c>`_
183
184Authorship and Verification
185----------------------------
186
187* **Author:** NIST IGOR/DANSE **Date:** pre 2010
188* **Converted to sasmodels by:** Miguel Gonzalez **Date:** February 26, 2016
189* **Last Modified by:** Paul Kienzle **Date:** October 17, 2017
190* **Last Reviewed by:** Paul Butler **Date:** May 24, 2018 - documentation
191  updated
192* **Source added by :** Steve King **Date:** March 25, 2019
193"""
194
195import numpy as np
196from numpy import inf
197
198name = "core_shell_parallelepiped"
199title = "Rectangular solid with a core-shell structure."
200description = """
201     P(q)=
202"""
203category = "shape:parallelepiped"
204
205#             ["name", "units", default, [lower, upper], "type","description"],
206parameters = [["sld_core", "1e-6/Ang^2", 1, [-inf, inf], "sld",
207               "Parallelepiped core scattering length density"],
208              ["sld_a", "1e-6/Ang^2", 2, [-inf, inf], "sld",
209               "Parallelepiped A rim scattering length density"],
210              ["sld_b", "1e-6/Ang^2", 4, [-inf, inf], "sld",
211               "Parallelepiped B rim scattering length density"],
212              ["sld_c", "1e-6/Ang^2", 2, [-inf, inf], "sld",
213               "Parallelepiped C rim scattering length density"],
214              ["sld_solvent", "1e-6/Ang^2", 6, [-inf, inf], "sld",
215               "Solvent scattering length density"],
216              ["length_a", "Ang", 35, [0, inf], "volume",
217               "Shorter side of the parallelepiped"],
218              ["length_b", "Ang", 75, [0, inf], "volume",
219               "Second side of the parallelepiped"],
220              ["length_c", "Ang", 400, [0, inf], "volume",
221               "Larger side of the parallelepiped"],
222              ["thick_rim_a", "Ang", 10, [0, inf], "volume",
223               "Thickness of A rim"],
224              ["thick_rim_b", "Ang", 10, [0, inf], "volume",
225               "Thickness of B rim"],
226              ["thick_rim_c", "Ang", 10, [0, inf], "volume",
227               "Thickness of C rim"],
228              ["theta", "degrees", 0, [-360, 360], "orientation",
229               "c axis to beam angle"],
230              ["phi", "degrees", 0, [-360, 360], "orientation",
231               "rotation about beam"],
232              ["psi", "degrees", 0, [-360, 360], "orientation",
233               "rotation about c axis"],
234             ]
235
236source = ["lib/gauss76.c", "core_shell_parallelepiped.c"]
237have_Fq = True
238effective_radius_type = [
239    "equivalent cylinder excluded volume",
240    "equivalent volume sphere",
241    "half outer length_a", "half outer length_b", "half outer length_c",
242    "equivalent circular cross-section",
243    "half outer ab diagonal", "half outer diagonal",
244    ]
245
246def random():
247    """Return a random parameter set for the model."""
248    outer = 10**np.random.uniform(1, 4.7, size=3)
249    thick = np.random.beta(0.5, 0.5, size=3)*(outer-2) + 1
250    length = outer - thick
251    pars = dict(
252        length_a=length[0],
253        length_b=length[1],
254        length_c=length[2],
255        thick_rim_a=thick[0],
256        thick_rim_b=thick[1],
257        thick_rim_c=thick[2],
258    )
259    return pars
260
261# parameters for demo
262demo = dict(scale=1, background=0.0,
263            sld_core=1, sld_a=2, sld_b=4, sld_c=2, sld_solvent=6,
264            length_a=35, length_b=75, length_c=400,
265            thick_rim_a=10, thick_rim_b=10, thick_rim_c=10,
266            theta=0, phi=0, psi=0,
267            length_a_pd=0.1, length_a_pd_n=1,
268            length_b_pd=0.1, length_b_pd_n=1,
269            length_c_pd=0.1, length_c_pd_n=1,
270            thick_rim_a_pd=0.1, thick_rim_a_pd_n=1,
271            thick_rim_b_pd=0.1, thick_rim_b_pd_n=1,
272            thick_rim_c_pd=0.1, thick_rim_c_pd_n=1,
273            theta_pd=10, theta_pd_n=1,
274            phi_pd=10, phi_pd_n=1,
275            psi_pd=10, psi_pd_n=1)
276
277# rkh 7/4/17 add random unit test for 2d, note make all params different,
278# 2d values not tested against other codes or models
279if 0:  # pak: model rewrite; need to update tests
280    qx, qy = 0.2 * np.cos(np.pi/6.), 0.2 * np.sin(np.pi/6.)
281    tests = [[{}, 0.2, 0.533149288477],
282             [{}, [0.2], [0.533149288477]],
283             [{'theta':10.0, 'phi':20.0}, (qx, qy), 0.0853299803222],
284             [{'theta':10.0, 'phi':20.0}, [(qx, qy)], [0.0853299803222]],
285            ]
286    del qx, qy  # not necessary to delete, but cleaner
Note: See TracBrowser for help on using the repository browser.