source: sasmodels/sasmodels/models/core_shell_parallelepiped.py @ 8f04da4

core_shell_microgelscostrafo411magnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 8f04da4 was 8f04da4, checked in by Paul Kienzle <pkienzle@…>, 7 years ago

tuned random model generation for more models

  • Property mode set to 100644
File size: 8.6 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
[cb0dc22]7"rim" can be different on each (pair) of faces. However at this time
[1f65db5]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.
[cb0dc22]10
11.. note::
[1916c52]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.
[44bd2be]14
[500128b]15The form factor is normalized by the particle volume $V$ such that
[44bd2be]16
[500128b]17.. math::
18
19    I(q) = \text{scale}\frac{\langle f^2 \rangle}{V} + \text{background}
[44bd2be]20
[500128b]21where $\langle \ldots \rangle$ is an average over all possible orientations
22of the rectangular solid.
[44bd2be]23
24
25The function calculated is the form factor of the rectangular solid below.
[500128b]26The core of the solid is defined by the dimensions $A$, $B$, $C$ such that
27$A < B < C$.
[44bd2be]28
[2f0c07d]29.. image:: img/core_shell_parallelepiped_geometry.jpg
[44bd2be]30
[500128b]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
[44bd2be]34
[1916c52]35.. image:: img/core_shell_parallelepiped_projection.jpg
[44bd2be]36
37The volume of the solid is
38
39.. math::
40
41    V = ABC + 2t_ABC + 2t_BAC + 2t_CAB
42
[cb0dc22]43**meaning that there are "gaps" at the corners of the solid.**  Again note that
[8f04da4]44$t_C = 0$ currently.
[44bd2be]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
48amplitudes of the core and shell, in the same manner as a core-shell model.
[44bd2be]49
[5810f00]50.. math::
51
52    F_{a}(Q,\alpha,\beta)=
[1916c52]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}
[14207bb]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]
[5810f00]57
58.. note::
59
[1916c52]60    Why does t_B not appear in the above equation?
[5810f00]61    For the calculation of the form factor to be valid, the sides of the solid
[1916c52]62    MUST (perhaps not any more?) be chosen such that** $A < B < C$.
[5810f00]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.
[44bd2be]65
66FITTING NOTES
[92dfe0c]67If the scale is set equal to the particle volume fraction, $\phi$, the returned
[500128b]68value is the scattered intensity per unit volume, $I(q) = \phi P(q)$.
[5810f00]69However, **no interparticle interference effects are included in this
70calculation.**
[44bd2be]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
[500128b]76$A < B < C$ is not violated. The calculation will not report an error,
[44bd2be]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})$
[15a90c1]83and length $(C+2t_C)$ values, after appropriately
[8f04da4]84sorting the three dimensions to give an oblate or prolate particle, to give an
[15a90c1]85effective radius, for $S(Q)$ when $P(Q) * S(Q)$ is applied.
[44bd2be]86
87To provide easy access to the orientation of the parallelepiped, we define the
[500128b]88axis of the cylinder using three angles $\theta$, $\phi$ and $\Psi$.
[2f0c07d]89(see :ref:`cylinder orientation <cylinder-angle-definition>`).
[500128b]90The angle $\Psi$ is the rotational angle around the *long_c* axis against the
91$q$ plane. For example, $\Psi = 0$ when the *short_b* axis is parallel to the
[44bd2be]92*x*-axis of the detector.
93
[15a90c1]94.. figure:: img/parallelepiped_angle_definition.png
[44bd2be]95
96    Definition of the angles for oriented core-shell parallelepipeds.
97
[1916c52]98.. figure:: img/parallelepiped_angle_projection.png
[44bd2be]99
100    Examples of the angles for oriented core-shell parallelepipeds against the
101    detector plane.
102
[aa2edb2]103References
104----------
[44bd2be]105
[5810f00]106.. [#] P Mittelbach and G Porod, *Acta Physica Austriaca*, 14 (1961) 185-211
107    Equations (1), (13-14). (in German)
108.. [#] D Singh (2009). *Small angle scattering studies of self assembly in
109   lipid mixtures*, John's Hopkins University Thesis (2009) 223-225. `Available
110   from Proquest <http://search.proquest.com/docview/304915826?accountid
111   =26379>`_
112
113Authorship and Verification
114----------------------------
[44bd2be]115
[5810f00]116* **Author:** NIST IGOR/DANSE **Date:** pre 2010
[cb0dc22]117* **Converted to sasmodels by:** Miguel Gonzales **Date:** February 26, 2016
118* **Last Modified by:** Wojciech Potrzebowski **Date:** January 11, 2017
119* **Currently Under review by:** Paul Butler
[44bd2be]120"""
121
122import numpy as np
[14207bb]123from numpy import pi, inf, sqrt, cos, sin
[44bd2be]124
125name = "core_shell_parallelepiped"
126title = "Rectangular solid with a core-shell structure."
127description = """
[8f04da4]128     P(q)=
[44bd2be]129"""
130category = "shape:parallelepiped"
131
132#             ["name", "units", default, [lower, upper], "type","description"],
[42356c8]133parameters = [["sld_core", "1e-6/Ang^2", 1, [-inf, inf], "sld",
[44bd2be]134               "Parallelepiped core scattering length density"],
[42356c8]135              ["sld_a", "1e-6/Ang^2", 2, [-inf, inf], "sld",
[44bd2be]136               "Parallelepiped A rim scattering length density"],
[42356c8]137              ["sld_b", "1e-6/Ang^2", 4, [-inf, inf], "sld",
[44bd2be]138               "Parallelepiped B rim scattering length density"],
[42356c8]139              ["sld_c", "1e-6/Ang^2", 2, [-inf, inf], "sld",
[44bd2be]140               "Parallelepiped C rim scattering length density"],
[42356c8]141              ["sld_solvent", "1e-6/Ang^2", 6, [-inf, inf], "sld",
[44bd2be]142               "Solvent scattering length density"],
[2222134]143              ["length_a", "Ang", 35, [0, inf], "volume",
[44bd2be]144               "Shorter side of the parallelepiped"],
[2222134]145              ["length_b", "Ang", 75, [0, inf], "volume",
[44bd2be]146               "Second side of the parallelepiped"],
[2222134]147              ["length_c", "Ang", 400, [0, inf], "volume",
[44bd2be]148               "Larger side of the parallelepiped"],
[2222134]149              ["thick_rim_a", "Ang", 10, [0, inf], "volume",
[44bd2be]150               "Thickness of A rim"],
[2222134]151              ["thick_rim_b", "Ang", 10, [0, inf], "volume",
[44bd2be]152               "Thickness of B rim"],
[2222134]153              ["thick_rim_c", "Ang", 10, [0, inf], "volume",
[44bd2be]154               "Thickness of C rim"],
[9b79f29]155              ["theta", "degrees", 0, [-360, 360], "orientation",
156               "c axis to beam angle"],
157              ["phi", "degrees", 0, [-360, 360], "orientation",
158               "rotation about beam"],
159              ["psi", "degrees", 0, [-360, 360], "orientation",
160               "rotation about c axis"],
[44bd2be]161             ]
162
[43b7eea]163source = ["lib/gauss76.c", "core_shell_parallelepiped.c"]
[44bd2be]164
165
[2222134]166def ER(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c):
[44bd2be]167    """
168        Return equivalent radius (ER)
169    """
170
171    # surface average radius (rough approximation)
[2222134]172    surf_rad = sqrt((length_a + 2.0*thick_rim_a) * (length_b + 2.0*thick_rim_b) / pi)
[44bd2be]173
[2222134]174    height = length_c + 2.0*thick_rim_c
[44bd2be]175
176    ddd = 0.75 * surf_rad * (2 * surf_rad * height + (height + surf_rad) * (height + pi * surf_rad))
177    return 0.5 * (ddd) ** (1. / 3.)
178
179# VR defaults to 1.0
180
[8f04da4]181def random():
182    import numpy as np
183    outer = 10**np.random.uniform(1, 4.7, size=3)
184    thick = np.random.beta(0.5, 0.5, size=3)*(outer-2) + 1
185    length = outer - thick
186    pars = dict(
187        length_a=length[0],
188        length_b=length[1],
189        length_c=length[2],
190        thick_rim_a=thick[0],
191        thick_rim_b=thick[1],
192        thick_rim_c=thick[2],
193    )
194    return pars
195
[44bd2be]196# parameters for demo
197demo = dict(scale=1, background=0.0,
[14838a3]198            sld_core=1, sld_a=2, sld_b=4, sld_c=2, sld_solvent=6,
[2222134]199            length_a=35, length_b=75, length_c=400,
200            thick_rim_a=10, thick_rim_b=10, thick_rim_c=10,
[44bd2be]201            theta=0, phi=0, psi=0,
[2222134]202            length_a_pd=0.1, length_a_pd_n=1,
203            length_b_pd=0.1, length_b_pd_n=1,
204            length_c_pd=0.1, length_c_pd_n=1,
205            thick_rim_a_pd=0.1, thick_rim_a_pd_n=1,
206            thick_rim_b_pd=0.1, thick_rim_b_pd_n=1,
207            thick_rim_c_pd=0.1, thick_rim_c_pd_n=1,
[44bd2be]208            theta_pd=10, theta_pd_n=1,
209            phi_pd=10, phi_pd_n=1,
[14838a3]210            psi_pd=10, psi_pd_n=1)
[44bd2be]211
[14207bb]212# rkh 7/4/17 add random unit test for 2d, note make all params different, 2d values not tested against other codes or models
213qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.)
[6dd90c1]214tests = [[{}, 0.2, 0.533149288477],
215         [{}, [0.2], [0.533149288477]],
[14207bb]216         [{'theta':10.0, 'phi':20.0}, (qx, qy), 0.0853299803222],
217         [{'theta':10.0, 'phi':20.0}, [(qx, qy)], [0.0853299803222]],
[44bd2be]218        ]
219del qx, qy  # not necessary to delete, but cleaner
Note: See TracBrowser for help on using the repository browser.