source: sasmodels/sasmodels/models/triaxial_ellipsoid.py @ 0b56f38

core_shell_microgelscostrafo411magnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 0b56f38 was 0b56f38, checked in by richardh, 7 years ago

unit tests added, needs opencl testing

  • Property mode set to 100644
File size: 6.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 all orientations $\Omega$ of:
19
20.. math::
21
22    P(q) = \text{scale}(\Delta\rho)^2\frac{V}{4 \pi}\int_\Omega \Phi^2(qr) d\Omega + \text{background}
23
24where
25
26.. math::
27
28    \Phi(qr) &= 3 j_1(qr)/qr = 3 (\sin qr - qr \cos qr)/(qr)^3 \\
29    r^2 &= R_a^2e^2 + R_b^2f^2 + R_c^2g^2 \\
30    V &= \tfrac{4}{3} \pi R_a R_b R_c
31
32The $e$, $f$ and $g$ terms are the projections of the orientation vector on $X$,
33$Y$ and $Z$ respectively.  Keeping the orientation fixed at the canonical
34axes, we can integrate over the incident direction using polar angle
35$-\pi/2 \le \gamma \le \pi/2$ and equatorial angle $0 \le \phi \le 2\pi$
36(as defined in ref [1]),
37
38 .. math::
39
40     \langle\Phi^2\rangle = \int_0^{2\pi} \int_{-\pi/2}^{\pi/2} \Phi^2(qr) \cos \gamma\,d\gamma d\phi
41
42with $e = \cos\gamma \sin\phi$, $f = \cos\gamma \cos\phi$ and $g = \sin\gamma$.
43A little algebra yields
44
45.. math::
46
47    r^2 = b^2(p_a \sin^2 \phi \cos^2 \gamma + 1 + p_c \sin^2 \gamma)
48
49for
50
51.. math::
52
53    p_a = \frac{a^2}{b^2} - 1 \text{ and } p_c = \frac{c^2}{b^2} - 1
54
55Due to symmetry, the ranges can be restricted to a single quadrant
56$0 \le \gamma \le \pi/2$ and $0 \le \phi \le \pi/2$, scaling the resulting
57integral by 8. The computation is done using the substitution $u = \sin\gamma$,
58$du = \cos\gamma\,d\gamma$, giving
59
60.. math::
61
62    \langle\Phi^2\rangle &= 8 \int_0^{\pi/2} \int_0^1 \Phi^2(qr) du d\phi \\
63    r^2 &= b^2(p_a \sin^2(\phi)(1 - u^2) + 1 + p_c u^2)
64
65To provide easy access to the orientation of the triaxial ellipsoid,
66we define the axis of the cylinder using the angles $\theta$, $\phi$
67and $\psi$. These angles are defined on
68:numref:`triaxial-ellipsoid-angles` .
69The angle $\psi$ is the rotational angle around its own $c$ axis
70against the $q$ plane. For example, $\psi = 0$ when the
71$a$ axis is parallel to the $x$ axis of the detector.
72
73.. _triaxial-ellipsoid-angles:
74
75.. figure:: img/triaxial_ellipsoid_angle_projection.jpg
76
77    The angles for oriented ellipsoid.
78
79The radius-of-gyration for this system is  $R_g^2 = (R_a R_b R_c)^2/5$.
80
81The contrast $\Delta\rho$ is defined as SLD(ellipsoid) - SLD(solvent).  In the
82parameters, $R_a$ is the minor equatorial radius, $R_b$ is the major
83equatorial radius, and $R_c$ is the polar radius of the ellipsoid.
84
85NB: The 2nd virial coefficient of the triaxial solid ellipsoid is
86calculated based on the polar radius $R_p = R_c$ and equatorial
87radius $R_e = \sqrt{R_a R_b}$, and used as the effective radius for
88$S(q)$ when $P(q) \cdot S(q)$ is applied.
89
90Validation
91----------
92
93Validation of our code was done by comparing the output of the
941D calculation to the angular average of the output of 2D calculation
95over all possible angles.
96
97
98References
99----------
100
101[1] Finnigan, J.A., Jacobs, D.J., 1971.
102*Light scattering by ellipsoidal particles in solution*,
103J. Phys. D: Appl. Phys. 4, 72-77. doi:10.1088/0022-3727/4/1/310
104
105Authorship and Verification
106----------------------------
107
108* **Author:** NIST IGOR/DANSE **Date:** pre 2010
109* **Last Modified by:** Paul Kienzle (improved calculation) **Date:** April 4, 2017
110* **Last Reviewed by:** Paul Kienzle &Richard Heenan **Date:**  April 4, 2017
111
112"""
113
114from numpy import inf, sin, cos, pi
115
116name = "triaxial_ellipsoid"
117title = "Ellipsoid of uniform scattering length density with three independent axes."
118
119description = """\
120Note: During fitting ensure that the inequality ra<rb<rc is not
121        violated. Otherwise the calculation will
122        not be correct.
123"""
124category = "shape:ellipsoid"
125
126#             ["name", "units", default, [lower, upper], "type","description"],
127parameters = [["sld", "1e-6/Ang^2", 4, [-inf, inf], "sld",
128               "Ellipsoid scattering length density"],
129              ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "sld",
130               "Solvent scattering length density"],
131              ["radius_equat_minor", "Ang", 20, [0, inf], "volume",
132               "Minor equatorial radius, Ra"],
133              ["radius_equat_major", "Ang", 400, [0, inf], "volume",
134               "Major equatorial radius, Rb"],
135              ["radius_polar", "Ang", 10, [0, inf], "volume",
136               "Polar radius, Rc"],
137              ["theta", "degrees", 60, [-inf, inf], "orientation",
138               "In plane angle"],
139              ["phi", "degrees", 60, [-inf, inf], "orientation",
140               "Out of plane angle"],
141              ["psi", "degrees", 60, [-inf, inf], "orientation",
142               "Out of plane angle"],
143             ]
144
145source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "triaxial_ellipsoid.c"]
146
147def ER(radius_equat_minor, radius_equat_major, radius_polar):
148    """
149        Returns the effective radius used in the S*P calculation
150    """
151    import numpy as np
152    from .ellipsoid import ER as ellipsoid_ER
153     # now that radii can be in any size order, radii need sorting a,b,c where a~b and c is either much smaller or much larger
154     # also need some unit tests!
155   
156    return ellipsoid_ER(radius_polar, np.sqrt(radius_equat_minor * radius_equat_major))
157
158demo = dict(scale=1, background=0,
159            sld=6, sld_solvent=1,
160            theta=30, phi=15, psi=5,
161            radius_equat_minor=25, radius_equat_major=36, radius_polar=50,
162            radius_equat_minor_pd=0, radius_equat_minor_pd_n=1,
163            radius_equat_major_pd=0, radius_equat_major_pd_n=1,
164            radius_polar_pd=.2, radius_polar_pd_n=30,
165            theta_pd=15, theta_pd_n=45,
166            phi_pd=15, phi_pd_n=1,
167            psi_pd=15, psi_pd_n=1)
168q = 0.1
169# april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct!
170# add 2d test after pull #890
171qx = q*cos(pi/6.0)
172qy = q*sin(pi/6.0)
173tests = [[{}, 0.05, 24.8839548033],
174#        [{'theta':80., 'phi':10.}, (qx, qy), 9999. ],
175        ]
176del qx, qy  # not necessary to delete, but cleaner
Note: See TracBrowser for help on using the repository browser.