source: sasmodels/sasmodels/models/triaxial_ellipsoid.py @ c1e44e5

Last change on this file since c1e44e5 was c1e44e5, checked in by Paul Kienzle <pkienzle@…>, 5 years ago

Add local link to source files. Refs #1263.

  • Property mode set to 100644
File size: 6.9 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
118Authorship and Verification
119----------------------------
120
121* **Author:** NIST IGOR/DANSE **Date:** pre 2010
122* **Last Modified by:** Paul Kienzle (improved calculation) **Date:** April 4, 2017
123* **Last Reviewed by:** Paul Kienzle & Richard Heenan **Date:**  April 4, 2017
124"""
125
126import numpy as np
127from numpy import inf, sin, cos, pi
128
129name = "triaxial_ellipsoid"
130title = "Ellipsoid of uniform scattering length density with three independent axes."
131
132description = """
133   Triaxial ellipsoid - see main documentation.
134"""
135category = "shape:ellipsoid"
136
137#             ["name", "units", default, [lower, upper], "type","description"],
138parameters = [["sld", "1e-6/Ang^2", 4, [-inf, inf], "sld",
139               "Ellipsoid scattering length density"],
140              ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "sld",
141               "Solvent scattering length density"],
142              ["radius_equat_minor", "Ang", 20, [0, inf], "volume",
143               "Minor equatorial radius, Ra"],
144              ["radius_equat_major", "Ang", 400, [0, inf], "volume",
145               "Major equatorial radius, Rb"],
146              ["radius_polar", "Ang", 10, [0, inf], "volume",
147               "Polar radius, Rc"],
148              ["theta", "degrees", 60, [-360, 360], "orientation",
149               "polar axis to beam angle"],
150              ["phi", "degrees", 60, [-360, 360], "orientation",
151               "rotation about beam"],
152              ["psi", "degrees", 60, [-360, 360], "orientation",
153               "rotation about polar axis"],
154             ]
155
156source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "triaxial_ellipsoid.c"]
157have_Fq = True
158effective_radius_type = [
159    "equivalent biaxial ellipsoid average curvature",
160    "equivalent volume sphere", "min radius", "max radius",
161    ]
162
163def random():
164    """Return a random parameter set for the model."""
165    a, b, c = 10**np.random.uniform(1, 4.7, size=3)
166    pars = dict(
167        radius_equat_minor=a,
168        radius_equat_major=b,
169        radius_polar=c,
170    )
171    return pars
172
173
174demo = dict(scale=1, background=0,
175            sld=6, sld_solvent=1,
176            theta=30, phi=15, psi=5,
177            radius_equat_minor=25, radius_equat_major=36, radius_polar=50,
178            radius_equat_minor_pd=0, radius_equat_minor_pd_n=1,
179            radius_equat_major_pd=0, radius_equat_major_pd_n=1,
180            radius_polar_pd=.2, radius_polar_pd_n=30,
181            theta_pd=15, theta_pd_n=45,
182            phi_pd=15, phi_pd_n=1,
183            psi_pd=15, psi_pd_n=1)
184
185q = 0.1
186# april 6 2017, rkh add unit tests
187#     NOT compared with any other calc method, assume correct!
188# check 2d test after pull #890
189qx = q*cos(pi/6.0)
190qy = q*sin(pi/6.0)
191tests = [
192    [{}, 0.05, 24.8839548033],
193    [{'theta':80., 'phi':10.}, (qx, qy), 166.712060266],
194    ]
195del qx, qy  # not necessary to delete, but cleaner
Note: See TracBrowser for help on using the repository browser.