source: sasmodels/sasmodels/models/triaxial_ellipsoid.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: 7.2 KB
Line 
1# triaxial ellipsoid model
2# Note: model title and parameter table are inserted automatically
3r"""
4Definition
5----------
6
7.. figure:: img/triaxial_ellipsoid_geometry.jpg
8
9    Ellipsoid with $R_a$ as *radius_equat_minor*, $R_b$ as *radius_equat_major*
10    and $R_c$ as *radius_polar*.
11
12Given an ellipsoid
13
14.. math::
15
16    \frac{X^2}{R_a^2} + \frac{Y^2}{R_b^2} + \frac{Z^2}{R_c^2} = 1
17
18the scattering for randomly oriented particles is defined by the average over
19all orientations $\Omega$ of:
20
21.. math::
22
23    P(q) = \text{scale}(\Delta\rho)^2\frac{V}{4 \pi}\int_\Omega\Phi^2(qr)\,d\Omega
24           + \text{background}
25
26where
27
28.. math::
29
30    \Phi(qr) &= 3 j_1(qr)/qr = 3 (\sin qr - qr \cos qr)/(qr)^3 \\
31    r^2 &= R_a^2e^2 + R_b^2f^2 + R_c^2g^2 \\
32    V &= \tfrac{4}{3} \pi R_a R_b R_c
33
34The $e$, $f$ and $g$ terms are the projections of the orientation vector on $X$,
35$Y$ and $Z$ respectively.  Keeping the orientation fixed at the canonical
36axes, we can integrate over the incident direction using polar angle
37$-\pi/2 \le \gamma \le \pi/2$ and equatorial angle $0 \le \phi \le 2\pi$
38(as defined in ref [1]),
39
40 .. math::
41
42     \langle\Phi^2\rangle = \int_0^{2\pi} \int_{-\pi/2}^{\pi/2} \Phi^2(qr)
43                                                \cos \gamma\,d\gamma d\phi
44
45with $e = \cos\gamma \sin\phi$, $f = \cos\gamma \cos\phi$ and $g = \sin\gamma$.
46A little algebra yields
47
48.. math::
49
50    r^2 = b^2(p_a \sin^2 \phi \cos^2 \gamma + 1 + p_c \sin^2 \gamma)
51
52for
53
54.. math::
55
56    p_a = \frac{a^2}{b^2} - 1 \text{ and } p_c = \frac{c^2}{b^2} - 1
57
58Due to symmetry, the ranges can be restricted to a single quadrant
59$0 \le \gamma \le \pi/2$ and $0 \le \phi \le \pi/2$, scaling the resulting
60integral by 8. The computation is done using the substitution $u = \sin\gamma$,
61$du = \cos\gamma\,d\gamma$, giving
62
63.. math::
64
65    \langle\Phi^2\rangle &= 8 \int_0^{\pi/2} \int_0^1 \Phi^2(qr) du d\phi \\
66    r^2 &= b^2(p_a \sin^2(\phi)(1 - u^2) + 1 + p_c u^2)
67
68Though for convenience we describe the three radii of the ellipsoid as equatorial
69and polar, they may be given in $any$ size order. To avoid multiple solutions, especially
70with Monte-Carlo fit methods, it may be advisable to restrict their ranges. For typical
71small angle diffraction situations there may be a number of closely similar "best fits",
72so some trial and error, or fixing of some radii at expected values, may help.
73
74To provide easy access to the orientation of the triaxial ellipsoid,
75we define the axis of the cylinder using the angles $\theta$, $\phi$
76and $\psi$. These angles are defined analogously to the elliptical_cylinder below, note that
77angle $\phi$ is now NOT the same as in the equations above.
78
79.. figure:: img/elliptical_cylinder_angle_definition.png
80
81    Definition of angles for oriented triaxial ellipsoid, where radii are for illustration here
82    $a < b << c$ and angle $\Psi$ is a rotation around the axis of the particle.
83
84For oriented ellipsoids the *theta*, *phi* and *psi* orientation parameters will appear when fitting 2D data,
85see the :ref:`elliptical-cylinder` model for further information.
86
87.. _triaxial-ellipsoid-angles:
88
89.. figure:: img/triaxial_ellipsoid_angle_projection.png
90
91    Some examples for an oriented triaxial ellipsoid.
92
93The radius-of-gyration for this system is  $R_g^2 = (R_a R_b R_c)^2/5$.
94
95The contrast $\Delta\rho$ is defined as SLD(ellipsoid) - SLD(solvent).  In the
96parameters, $R_a$ is the minor equatorial radius, $R_b$ is the major
97equatorial radius, and $R_c$ is the polar radius of the ellipsoid.
98
99NB: The 2nd virial coefficient of the triaxial solid ellipsoid is
100calculated after sorting the three radii to give the most appropriate
101prolate or oblate form, from the new polar radius $R_p = R_c$ and effective equatorial
102radius,  $R_e = \sqrt{R_a R_b}$, to then be used as the effective radius for
103$S(q)$ when $P(q) \cdot S(q)$ is applied.
104
105Validation
106----------
107
108Validation of our code was done by comparing the output of the
1091D calculation to the angular average of the output of 2D calculation
110over all possible angles.
111
112
113References
114----------
115
116.. [#] Finnigan, J.A., Jacobs, D.J., 1971. *Light scattering by ellipsoidal particles in solution*, J. Phys. D: Appl. Phys. 4, 72-77. doi:10.1088/0022-3727/4/1/310
117
118Source
119------
120
121`triaxial_ellipsoid.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/triaxial_ellipsoid.py>`_
122
123`triaxial_ellipsoid.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/triaxial_ellipsoid.c>`_
124
125Authorship and Verification
126----------------------------
127
128* **Author:** NIST IGOR/DANSE **Date:** pre 2010
129* **Last Modified by:** Paul Kienzle (improved calculation) **Date:** April 4, 2017
130* **Last Reviewed by:** Paul Kienzle & Richard Heenan **Date:**  April 4, 2017
131* **Source added by :** Steve King **Date:** March 25, 2019
132"""
133
134import numpy as np
135from numpy import inf, sin, cos, pi
136
137name = "triaxial_ellipsoid"
138title = "Ellipsoid of uniform scattering length density with three independent axes."
139
140description = """
141   Triaxial ellipsoid - see main documentation.
142"""
143category = "shape:ellipsoid"
144
145#             ["name", "units", default, [lower, upper], "type","description"],
146parameters = [["sld", "1e-6/Ang^2", 4, [-inf, inf], "sld",
147               "Ellipsoid scattering length density"],
148              ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "sld",
149               "Solvent scattering length density"],
150              ["radius_equat_minor", "Ang", 20, [0, inf], "volume",
151               "Minor equatorial radius, Ra"],
152              ["radius_equat_major", "Ang", 400, [0, inf], "volume",
153               "Major equatorial radius, Rb"],
154              ["radius_polar", "Ang", 10, [0, inf], "volume",
155               "Polar radius, Rc"],
156              ["theta", "degrees", 60, [-360, 360], "orientation",
157               "polar axis to beam angle"],
158              ["phi", "degrees", 60, [-360, 360], "orientation",
159               "rotation about beam"],
160              ["psi", "degrees", 60, [-360, 360], "orientation",
161               "rotation about polar axis"],
162             ]
163
164source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "triaxial_ellipsoid.c"]
165have_Fq = True
166effective_radius_type = [
167    "equivalent biaxial ellipsoid average curvature",
168    "equivalent volume sphere", "min radius", "max radius",
169    ]
170
171def random():
172    """Return a random parameter set for the model."""
173    a, b, c = 10**np.random.uniform(1, 4.7, size=3)
174    pars = dict(
175        radius_equat_minor=a,
176        radius_equat_major=b,
177        radius_polar=c,
178    )
179    return pars
180
181
182demo = dict(scale=1, background=0,
183            sld=6, sld_solvent=1,
184            theta=30, phi=15, psi=5,
185            radius_equat_minor=25, radius_equat_major=36, radius_polar=50,
186            radius_equat_minor_pd=0, radius_equat_minor_pd_n=1,
187            radius_equat_major_pd=0, radius_equat_major_pd_n=1,
188            radius_polar_pd=.2, radius_polar_pd_n=30,
189            theta_pd=15, theta_pd_n=45,
190            phi_pd=15, phi_pd_n=1,
191            psi_pd=15, psi_pd_n=1)
192
193q = 0.1
194# april 6 2017, rkh add unit tests
195#     NOT compared with any other calc method, assume correct!
196# check 2d test after pull #890
197qx = q*cos(pi/6.0)
198qy = q*sin(pi/6.0)
199tests = [
200    [{}, 0.05, 24.8839548033],
201    [{'theta':80., 'phi':10.}, (qx, qy), 166.712060266],
202    ]
203del qx, qy  # not necessary to delete, but cleaner
Note: See TracBrowser for help on using the repository browser.