source: sasmodels/sasmodels/models/core_shell_parallelepiped.py @ 5bc6d21

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

more standardization of text

moved validation to its own section (is that the standard?) and moved
parallelipiped frong stuff around so start with Definition as per
standard.

  • Property mode set to 100644
File size: 9.9 KB
RevLine 
[44bd2be]1r"""
[5810f00]2Definition
3----------
4
[44bd2be]5Calculates the form factor for a rectangular solid with a core-shell structure.
[8f04da4]6The thickness and the scattering length density of the shell or
[393facf]7"rim" can be different on each (pair) of faces.
[cb0dc22]8
[500128b]9The form factor is normalized by the particle volume $V$ such that
[44bd2be]10
[500128b]11.. math::
12
[dbf1a60]13    I(q) = \frac{\text{scale}}{V} \langle P(q,\alpha,\beta) \rangle
[367886f]14    + \text{background}
[44bd2be]15
[500128b]16where $\langle \ldots \rangle$ is an average over all possible orientations
[dbf1a60]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.
[44bd2be]19
[500128b]20The core of the solid is defined by the dimensions $A$, $B$, $C$ such that
21$A < B < C$.
[44bd2be]22
[5bc373b]23.. figure:: img/parallelepiped_geometry.jpg
24
[331870d]25   Core of the core shell parallelepiped with the corresponding definition
[5bc373b]26   of sides.
27
[44bd2be]28
[500128b]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$
[dbf1a60]31$(=t_C)$ faces. The projection in the $AB$ plane is
[44bd2be]32
[5bc373b]33.. figure:: img/core_shell_parallelepiped_projection.jpg
34
[331870d]35   AB cut through the core-shell parallelipiped showing the cross secion of
36   four of the six shell slabs. As can be seen, this model leaves **"gaps"**
[dbf1a60]37   at the corners of the solid.
[44bd2be]38
[dbf1a60]39
40The total volume of the solid is thus given as
[44bd2be]41
42.. math::
43
44    V = ABC + 2t_ABC + 2t_BAC + 2t_CAB
45
[5810f00]46The intensity calculated follows the :ref:`parallelepiped` model, with the
47core-shell intensity being calculated as the square of the sum of the
[dbf1a60]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.
[44bd2be]54
[5810f00]55.. math::
[4493288]56
[dbf1a60]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}
[367886f]59
60and
61
62.. math::
63
[dbf1a60]64    F(q,\alpha,\beta)
[4493288]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})
[5bc373b]68        \left[S(Q_A, A+2t_A) - S(Q_A, A)\right] S(Q_B, B) S(Q_C, C) \\
[4493288]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]
[393facf]73
74with
[5810f00]75
[393facf]76.. math::
[5810f00]77
[331870d]78    S(Q_X, L) = L \frac{\sin (\tfrac{1}{2} Q_X L)}{\tfrac{1}{2} Q_X L}
[4493288]79
80and
81
82.. math::
[5810f00]83
[5bc373b]84    Q_A &= q \sin\alpha \sin\beta \\
85    Q_B &= q \sin\alpha \cos\beta \\
86    Q_C &= q \cos\alpha
[4493288]87
88
89where $\rho_\text{core}$, $\rho_\text{A}$, $\rho_\text{B}$ and $\rho_\text{C}$
[331870d]90are the scattering lengths of the parallelepiped core, and the rectangular
[4493288]91slabs of thickness $t_A$, $t_B$ and $t_C$, respectively. $\rho_\text{solvent}$
92is the scattering length of the solvent.
[44bd2be]93
[dbf1a60]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
[44bd2be]103FITTING NOTES
[4493288]104~~~~~~~~~~~~~
105
[92dfe0c]106If the scale is set equal to the particle volume fraction, $\phi$, the returned
[4493288]107value is the scattered intensity per unit volume, $I(q) = \phi P(q)$. However,
108**no interparticle interference effects are included in this calculation.**
[44bd2be]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})$
[4493288]117and length $(C+2t_C)$ values, after appropriately sorting the three dimensions
[331870d]118to give an oblate or prolate particle, to give an effective radius
[5bc373b]119for $S(q)$ when $P(q) * S(q)$ is applied.
[44bd2be]120
[904cd9c]121For 2d data the orientation of the particle is required, described using
[331870d]122angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below. For further
[4493288]123details of the calculation and angular dispersions see :ref:`orientation`.
[904cd9c]124The angle $\Psi$ is the rotational angle around the *long_c* axis. For example,
[eda8b30]125$\Psi = 0$ when the *short_b* axis is parallel to the *x*-axis of the detector.
[44bd2be]126
[dbf1a60]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.
[393facf]131
[15a90c1]132.. figure:: img/parallelepiped_angle_definition.png
[44bd2be]133
134    Definition of the angles for oriented core-shell parallelepipeds.
[2d81cfe]135    Note that rotation $\theta$, initially in the $xz$ plane, is carried
136    out first, then rotation $\phi$ about the $z$ axis, finally rotation
[331870d]137    $\Psi$ is now around the axis of the particle. The neutron or X-ray
[2d81cfe]138    beam is along the $z$ axis.
[44bd2be]139
[1916c52]140.. figure:: img/parallelepiped_angle_projection.png
[44bd2be]141
142    Examples of the angles for oriented core-shell parallelepipeds against the
143    detector plane.
144
[5bc6d21]145
146Validation
147----------
148
149Cross-checked against hollow rectangular prism and rectangular prism for equal
150thickness overlapping sides, and by Monte Carlo sampling of points within the
151shape for non-uniform, non-overlapping sides.
152
153
[aa2edb2]154References
155----------
[44bd2be]156
[5810f00]157.. [#] P Mittelbach and G Porod, *Acta Physica Austriaca*, 14 (1961) 185-211
158    Equations (1), (13-14). (in German)
159.. [#] D Singh (2009). *Small angle scattering studies of self assembly in
[fc0b7aa]160   lipid mixtures*, Johns Hopkins University Thesis (2009) 223-225. `Available
[5810f00]161   from Proquest <http://search.proquest.com/docview/304915826?accountid
162   =26379>`_
163
164Authorship and Verification
165----------------------------
[44bd2be]166
[5810f00]167* **Author:** NIST IGOR/DANSE **Date:** pre 2010
[331870d]168* **Converted to sasmodels by:** Miguel Gonzalez **Date:** February 26, 2016
[97be877]169* **Last Modified by:** Paul Kienzle **Date:** October 17, 2017
[44bd2be]170"""
171
172import numpy as np
[14207bb]173from numpy import pi, inf, sqrt, cos, sin
[44bd2be]174
175name = "core_shell_parallelepiped"
176title = "Rectangular solid with a core-shell structure."
177description = """
[8f04da4]178     P(q)=
[44bd2be]179"""
180category = "shape:parallelepiped"
181
182#             ["name", "units", default, [lower, upper], "type","description"],
[42356c8]183parameters = [["sld_core", "1e-6/Ang^2", 1, [-inf, inf], "sld",
[44bd2be]184               "Parallelepiped core scattering length density"],
[42356c8]185              ["sld_a", "1e-6/Ang^2", 2, [-inf, inf], "sld",
[44bd2be]186               "Parallelepiped A rim scattering length density"],
[42356c8]187              ["sld_b", "1e-6/Ang^2", 4, [-inf, inf], "sld",
[44bd2be]188               "Parallelepiped B rim scattering length density"],
[42356c8]189              ["sld_c", "1e-6/Ang^2", 2, [-inf, inf], "sld",
[44bd2be]190               "Parallelepiped C rim scattering length density"],
[42356c8]191              ["sld_solvent", "1e-6/Ang^2", 6, [-inf, inf], "sld",
[44bd2be]192               "Solvent scattering length density"],
[2222134]193              ["length_a", "Ang", 35, [0, inf], "volume",
[44bd2be]194               "Shorter side of the parallelepiped"],
[2222134]195              ["length_b", "Ang", 75, [0, inf], "volume",
[44bd2be]196               "Second side of the parallelepiped"],
[2222134]197              ["length_c", "Ang", 400, [0, inf], "volume",
[44bd2be]198               "Larger side of the parallelepiped"],
[2222134]199              ["thick_rim_a", "Ang", 10, [0, inf], "volume",
[44bd2be]200               "Thickness of A rim"],
[2222134]201              ["thick_rim_b", "Ang", 10, [0, inf], "volume",
[44bd2be]202               "Thickness of B rim"],
[2222134]203              ["thick_rim_c", "Ang", 10, [0, inf], "volume",
[44bd2be]204               "Thickness of C rim"],
[9b79f29]205              ["theta", "degrees", 0, [-360, 360], "orientation",
206               "c axis to beam angle"],
207              ["phi", "degrees", 0, [-360, 360], "orientation",
208               "rotation about beam"],
209              ["psi", "degrees", 0, [-360, 360], "orientation",
210               "rotation about c axis"],
[44bd2be]211             ]
212
[43b7eea]213source = ["lib/gauss76.c", "core_shell_parallelepiped.c"]
[44bd2be]214
215
[2222134]216def ER(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c):
[44bd2be]217    """
218        Return equivalent radius (ER)
219    """
[10ee838]220    from .parallelepiped import ER as ER_p
[44bd2be]221
[10ee838]222    a = length_a + 2*thick_rim_a
223    b = length_b + 2*thick_rim_b
224    c = length_c + 2*thick_rim_c
225    return ER_p(a, b, c)
[44bd2be]226
227# VR defaults to 1.0
228
[8f04da4]229def random():
230    outer = 10**np.random.uniform(1, 4.7, size=3)
231    thick = np.random.beta(0.5, 0.5, size=3)*(outer-2) + 1
232    length = outer - thick
233    pars = dict(
234        length_a=length[0],
235        length_b=length[1],
236        length_c=length[2],
237        thick_rim_a=thick[0],
238        thick_rim_b=thick[1],
239        thick_rim_c=thick[2],
240    )
241    return pars
242
[44bd2be]243# parameters for demo
244demo = dict(scale=1, background=0.0,
[14838a3]245            sld_core=1, sld_a=2, sld_b=4, sld_c=2, sld_solvent=6,
[2222134]246            length_a=35, length_b=75, length_c=400,
247            thick_rim_a=10, thick_rim_b=10, thick_rim_c=10,
[44bd2be]248            theta=0, phi=0, psi=0,
[2222134]249            length_a_pd=0.1, length_a_pd_n=1,
250            length_b_pd=0.1, length_b_pd_n=1,
251            length_c_pd=0.1, length_c_pd_n=1,
252            thick_rim_a_pd=0.1, thick_rim_a_pd_n=1,
253            thick_rim_b_pd=0.1, thick_rim_b_pd_n=1,
254            thick_rim_c_pd=0.1, thick_rim_c_pd_n=1,
[44bd2be]255            theta_pd=10, theta_pd_n=1,
256            phi_pd=10, phi_pd_n=1,
[14838a3]257            psi_pd=10, psi_pd_n=1)
[44bd2be]258
[4493288]259# rkh 7/4/17 add random unit test for 2d, note make all params different,
260# 2d values not tested against other codes or models
[fa70e04]261if 0:  # pak: model rewrite; need to update tests
262    qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.)
263    tests = [[{}, 0.2, 0.533149288477],
[2d81cfe]264             [{}, [0.2], [0.533149288477]],
265             [{'theta':10.0, 'phi':20.0}, (qx, qy), 0.0853299803222],
266             [{'theta':10.0, 'phi':20.0}, [(qx, qy)], [0.0853299803222]],
[fa70e04]267            ]
268    del qx, qy  # not necessary to delete, but cleaner
Note: See TracBrowser for help on using the repository browser.