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

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

core_shell_parallelepiped: simplify the calculation and update the docs

  • Property mode set to 100644
File size: 9.1 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
13    I(q) = \text{scale}\frac{\langle f^2 \rangle}{V} + \text{background}
[44bd2be]14
[500128b]15where $\langle \ldots \rangle$ is an average over all possible orientations
16of the rectangular solid.
[44bd2be]17
18The function calculated is the form factor of the rectangular solid below.
[500128b]19The core of the solid is defined by the dimensions $A$, $B$, $C$ such that
20$A < B < C$.
[44bd2be]21
[2f0c07d]22.. image:: img/core_shell_parallelepiped_geometry.jpg
[44bd2be]23
[500128b]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
[44bd2be]27
[1916c52]28.. image:: img/core_shell_parallelepiped_projection.jpg
[44bd2be]29
30The volume of the solid is
31
32.. math::
33
34    V = ABC + 2t_ABC + 2t_BAC + 2t_CAB
35
[393facf]36**meaning that there are "gaps" at the corners of the solid.**
[44bd2be]37
[5810f00]38The intensity calculated follows the :ref:`parallelepiped` model, with the
39core-shell intensity being calculated as the square of the sum of the
[393facf]40amplitudes of the core and the slabs on the edges.
41
[4493288]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.
[44bd2be]48
[5810f00]49.. math::
[4493288]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]
[393facf]60
61with
[5810f00]62
[393facf]63.. math::
[5810f00]64
[4493288]65    S(Q, L) = L \frac{\sin \tfrac{1}{2} Q L}{\tfrac{1}{2} Q L}
66
67and
68
69.. math::
[5810f00]70
[4493288]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.
[44bd2be]80
81FITTING NOTES
[4493288]82~~~~~~~~~~~~~
83
[92dfe0c]84If the scale is set equal to the particle volume fraction, $\phi$, the returned
[4493288]85value is the scattered intensity per unit volume, $I(q) = \phi P(q)$. However,
86**no interparticle interference effects are included in this calculation.**
[44bd2be]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})$
[4493288]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.
[44bd2be]98
[904cd9c]99For 2d data the orientation of the particle is required, described using
[4493288]100angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further
101details of the calculation and angular dispersions see :ref:`orientation`.
[904cd9c]102The angle $\Psi$ is the rotational angle around the *long_c* axis. For example,
[eda8b30]103$\Psi = 0$ when the *short_b* axis is parallel to the *x*-axis of the detector.
[44bd2be]104
[4493288]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,
[393facf]108but the results may be not correct.
109
[15a90c1]110.. figure:: img/parallelepiped_angle_definition.png
[44bd2be]111
112    Definition of the angles for oriented core-shell parallelepipeds.
[2d81cfe]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.
[44bd2be]117
[1916c52]118.. figure:: img/parallelepiped_angle_projection.png
[44bd2be]119
120    Examples of the angles for oriented core-shell parallelepipeds against the
121    detector plane.
122
[aa2edb2]123References
124----------
[44bd2be]125
[5810f00]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
[fc0b7aa]129   lipid mixtures*, Johns Hopkins University Thesis (2009) 223-225. `Available
[5810f00]130   from Proquest <http://search.proquest.com/docview/304915826?accountid
131   =26379>`_
132
133Authorship and Verification
134----------------------------
[44bd2be]135
[5810f00]136* **Author:** NIST IGOR/DANSE **Date:** pre 2010
[cb0dc22]137* **Converted to sasmodels by:** Miguel Gonzales **Date:** February 26, 2016
138* **Last Modified by:** Wojciech Potrzebowski **Date:** January 11, 2017
139* **Currently Under review by:** Paul Butler
[44bd2be]140"""
141
142import numpy as np
[14207bb]143from numpy import pi, inf, sqrt, cos, sin
[44bd2be]144
145name = "core_shell_parallelepiped"
146title = "Rectangular solid with a core-shell structure."
147description = """
[8f04da4]148     P(q)=
[44bd2be]149"""
150category = "shape:parallelepiped"
151
152#             ["name", "units", default, [lower, upper], "type","description"],
[42356c8]153parameters = [["sld_core", "1e-6/Ang^2", 1, [-inf, inf], "sld",
[44bd2be]154               "Parallelepiped core scattering length density"],
[42356c8]155              ["sld_a", "1e-6/Ang^2", 2, [-inf, inf], "sld",
[44bd2be]156               "Parallelepiped A rim scattering length density"],
[42356c8]157              ["sld_b", "1e-6/Ang^2", 4, [-inf, inf], "sld",
[44bd2be]158               "Parallelepiped B rim scattering length density"],
[42356c8]159              ["sld_c", "1e-6/Ang^2", 2, [-inf, inf], "sld",
[44bd2be]160               "Parallelepiped C rim scattering length density"],
[42356c8]161              ["sld_solvent", "1e-6/Ang^2", 6, [-inf, inf], "sld",
[44bd2be]162               "Solvent scattering length density"],
[2222134]163              ["length_a", "Ang", 35, [0, inf], "volume",
[44bd2be]164               "Shorter side of the parallelepiped"],
[2222134]165              ["length_b", "Ang", 75, [0, inf], "volume",
[44bd2be]166               "Second side of the parallelepiped"],
[2222134]167              ["length_c", "Ang", 400, [0, inf], "volume",
[44bd2be]168               "Larger side of the parallelepiped"],
[2222134]169              ["thick_rim_a", "Ang", 10, [0, inf], "volume",
[44bd2be]170               "Thickness of A rim"],
[2222134]171              ["thick_rim_b", "Ang", 10, [0, inf], "volume",
[44bd2be]172               "Thickness of B rim"],
[2222134]173              ["thick_rim_c", "Ang", 10, [0, inf], "volume",
[44bd2be]174               "Thickness of C rim"],
[9b79f29]175              ["theta", "degrees", 0, [-360, 360], "orientation",
176               "c axis to beam angle"],
177              ["phi", "degrees", 0, [-360, 360], "orientation",
178               "rotation about beam"],
179              ["psi", "degrees", 0, [-360, 360], "orientation",
180               "rotation about c axis"],
[44bd2be]181             ]
182
[43b7eea]183source = ["lib/gauss76.c", "core_shell_parallelepiped.c"]
[44bd2be]184
185
[2222134]186def ER(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c):
[44bd2be]187    """
188        Return equivalent radius (ER)
189    """
190
191    # surface average radius (rough approximation)
[4493288]192    surf_rad = sqrt((length_a + 2.0*thick_rim_a)
193                    * (length_b + 2.0*thick_rim_b) / pi)
[44bd2be]194
[2222134]195    height = length_c + 2.0*thick_rim_c
[44bd2be]196
[4493288]197    ddd = (0.75 * surf_rad * (2 * surf_rad * height
198           + (height + surf_rad) * (height + pi * surf_rad)))
[44bd2be]199    return 0.5 * (ddd) ** (1. / 3.)
200
201# VR defaults to 1.0
202
[8f04da4]203def random():
204    outer = 10**np.random.uniform(1, 4.7, size=3)
205    thick = np.random.beta(0.5, 0.5, size=3)*(outer-2) + 1
206    length = outer - thick
207    pars = dict(
208        length_a=length[0],
209        length_b=length[1],
210        length_c=length[2],
211        thick_rim_a=thick[0],
212        thick_rim_b=thick[1],
213        thick_rim_c=thick[2],
214    )
215    return pars
216
[44bd2be]217# parameters for demo
218demo = dict(scale=1, background=0.0,
[14838a3]219            sld_core=1, sld_a=2, sld_b=4, sld_c=2, sld_solvent=6,
[2222134]220            length_a=35, length_b=75, length_c=400,
221            thick_rim_a=10, thick_rim_b=10, thick_rim_c=10,
[44bd2be]222            theta=0, phi=0, psi=0,
[2222134]223            length_a_pd=0.1, length_a_pd_n=1,
224            length_b_pd=0.1, length_b_pd_n=1,
225            length_c_pd=0.1, length_c_pd_n=1,
226            thick_rim_a_pd=0.1, thick_rim_a_pd_n=1,
227            thick_rim_b_pd=0.1, thick_rim_b_pd_n=1,
228            thick_rim_c_pd=0.1, thick_rim_c_pd_n=1,
[44bd2be]229            theta_pd=10, theta_pd_n=1,
230            phi_pd=10, phi_pd_n=1,
[14838a3]231            psi_pd=10, psi_pd_n=1)
[44bd2be]232
[4493288]233# rkh 7/4/17 add random unit test for 2d, note make all params different,
234# 2d values not tested against other codes or models
[fa70e04]235if 0:  # pak: model rewrite; need to update tests
236    qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.)
237    tests = [[{}, 0.2, 0.533149288477],
[2d81cfe]238             [{}, [0.2], [0.533149288477]],
239             [{'theta':10.0, 'phi':20.0}, (qx, qy), 0.0853299803222],
240             [{'theta':10.0, 'phi':20.0}, [(qx, qy)], [0.0853299803222]],
[fa70e04]241            ]
242    del qx, qy  # not necessary to delete, but cleaner
Note: See TracBrowser for help on using the repository browser.