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

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

lint

  • Property mode set to 100644
File size: 7.3 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[1] Finnigan, J.A., Jacobs, D.J., 1971.
117*Light scattering by ellipsoidal particles in solution*,
118J. Phys. D: Appl. Phys. 4, 72-77. doi:10.1088/0022-3727/4/1/310
119
120Authorship and Verification
121----------------------------
122
123* **Author:** NIST IGOR/DANSE **Date:** pre 2010
124* **Last Modified by:** Paul Kienzle (improved calculation) **Date:** April 4, 2017
125* **Last Reviewed by:** Paul Kienzle & Richard Heenan **Date:**  April 4, 2017
126"""
127
128import numpy as np
129from numpy import inf, sin, cos, pi
130
131name = "triaxial_ellipsoid"
132title = "Ellipsoid of uniform scattering length density with three independent axes."
133
134description = """
135   Triaxial ellipsoid - see main documentation.
136"""
137category = "shape:ellipsoid"
138
139#             ["name", "units", default, [lower, upper], "type","description"],
140parameters = [["sld", "1e-6/Ang^2", 4, [-inf, inf], "sld",
141               "Ellipsoid scattering length density"],
142              ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "sld",
143               "Solvent scattering length density"],
144              ["radius_equat_minor", "Ang", 20, [0, inf], "volume",
145               "Minor equatorial radius, Ra"],
146              ["radius_equat_major", "Ang", 400, [0, inf], "volume",
147               "Major equatorial radius, Rb"],
148              ["radius_polar", "Ang", 10, [0, inf], "volume",
149               "Polar radius, Rc"],
150              ["theta", "degrees", 60, [-360, 360], "orientation",
151               "polar axis to beam angle"],
152              ["phi", "degrees", 60, [-360, 360], "orientation",
153               "rotation about beam"],
154              ["psi", "degrees", 60, [-360, 360], "orientation",
155               "rotation about polar axis"],
156             ]
157
158source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "triaxial_ellipsoid.c"]
159
160def ER(radius_equat_minor, radius_equat_major, radius_polar):
161    """
162    Returns the effective radius used in the S*P calculation
163    """
164    from .ellipsoid import ER as ellipsoid_ER
165
166    # now that radii can be in any size order, radii need sorting a,b,c
167    # where a~b and c is either much smaller or much larger
168    radii = np.vstack((radius_equat_major, radius_equat_minor, radius_polar))
169    radii = np.sort(radii, axis=0)
170    selector = (radii[1] - radii[0]) > (radii[2] - radii[1])
171    polar = np.where(selector, radii[0], radii[2])
172    equatorial = np.sqrt(np.where(~selector, radii[0]*radii[1], radii[1]*radii[2]))
173    return ellipsoid_ER(polar, equatorial)
174
175def random():
176    a, b, c = 10**np.random.uniform(1, 4.7, size=3)
177    pars = dict(
178        radius_equat_minor=a,
179        radius_equat_major=b,
180        radius_polar=c,
181    )
182    return pars
183
184
185demo = dict(scale=1, background=0,
186            sld=6, sld_solvent=1,
187            theta=30, phi=15, psi=5,
188            radius_equat_minor=25, radius_equat_major=36, radius_polar=50,
189            radius_equat_minor_pd=0, radius_equat_minor_pd_n=1,
190            radius_equat_major_pd=0, radius_equat_major_pd_n=1,
191            radius_polar_pd=.2, radius_polar_pd_n=30,
192            theta_pd=15, theta_pd_n=45,
193            phi_pd=15, phi_pd_n=1,
194            psi_pd=15, psi_pd_n=1)
195
196q = 0.1
197# april 6 2017, rkh add unit tests
198#     NOT compared with any other calc method, assume correct!
199# check 2d test after pull #890
200qx = q*cos(pi/6.0)
201qy = q*sin(pi/6.0)
202tests = [
203    [{}, 0.05, 24.8839548033],
204    [{'theta':80., 'phi':10.}, (qx, qy), 166.712060266],
205    ]
206del qx, qy  # not necessary to delete, but cleaner
Note: See TracBrowser for help on using the repository browser.