source: sasmodels/sasmodels/models/triaxial_ellipsoid.py @ 1f65db5

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

docs edits for bcc_ fc_ sc_paracrystal

  • Property mode set to 100644
File size: 6.5 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 analogously to the elliptical_cylinder below
68
69.. figure:: img/elliptical_cylinder_angle_definition.png
70
71    Definition of angles for oriented triaxial ellipsoid, where radii shown here are $a < b << c$
72    and angle $\Psi$ is a rotation around the axis of the particle.
73
74The angle $\psi$ is the rotational angle around its own $c$ axis
75against the $q$ plane. For example, $\psi = 0$ when the
76$a$ axis is parallel to the $x$ axis of the detector.
77
78.. _triaxial-ellipsoid-angles:
79
80.. figure:: img/triaxial_ellipsoid_angle_projection.png
81
82    Some example angles for oriented ellipsoid.
83
84The radius-of-gyration for this system is  $R_g^2 = (R_a R_b R_c)^2/5$.
85
86The contrast $\Delta\rho$ is defined as SLD(ellipsoid) - SLD(solvent).  In the
87parameters, $R_a$ is the minor equatorial radius, $R_b$ is the major
88equatorial radius, and $R_c$ is the polar radius of the ellipsoid.
89
90NB: The 2nd virial coefficient of the triaxial solid ellipsoid is
91calculated based on the polar radius $R_p = R_c$ and equatorial
92radius $R_e = \sqrt{R_a R_b}$, and used as the effective radius for
93$S(q)$ when $P(q) \cdot S(q)$ is applied.
94
95Validation
96----------
97
98Validation of our code was done by comparing the output of the
991D calculation to the angular average of the output of 2D calculation
100over all possible angles.
101
102
103References
104----------
105
106[1] Finnigan, J.A., Jacobs, D.J., 1971.
107*Light scattering by ellipsoidal particles in solution*,
108J. Phys. D: Appl. Phys. 4, 72-77. doi:10.1088/0022-3727/4/1/310
109
110Authorship and Verification
111----------------------------
112
113* **Author:** NIST IGOR/DANSE **Date:** pre 2010
114* **Last Modified by:** Paul Kienzle (improved calculation) **Date:** April 4, 2017
115* **Last Reviewed by:** Paul Kienzle &Richard Heenan **Date:**  April 4, 2017
116
117"""
118
119from numpy import inf, sin, cos, pi
120
121name = "triaxial_ellipsoid"
122title = "Ellipsoid of uniform scattering length density with three independent axes."
123
124description = """\
125Note: During fitting ensure that the inequality ra<rb<rc is not
126        violated. Otherwise the calculation will
127        not be correct.
128"""
129category = "shape:ellipsoid"
130
131#             ["name", "units", default, [lower, upper], "type","description"],
132parameters = [["sld", "1e-6/Ang^2", 4, [-inf, inf], "sld",
133               "Ellipsoid scattering length density"],
134              ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "sld",
135               "Solvent scattering length density"],
136              ["radius_equat_minor", "Ang", 20, [0, inf], "volume",
137               "Minor equatorial radius, Ra"],
138              ["radius_equat_major", "Ang", 400, [0, inf], "volume",
139               "Major equatorial radius, Rb"],
140              ["radius_polar", "Ang", 10, [0, inf], "volume",
141               "Polar radius, Rc"],
142              ["theta", "degrees", 60, [-inf, inf], "orientation",
143               "In plane angle"],
144              ["phi", "degrees", 60, [-inf, inf], "orientation",
145               "Out of plane angle"],
146              ["psi", "degrees", 60, [-inf, inf], "orientation",
147               "Out of plane angle"],
148             ]
149
150source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "triaxial_ellipsoid.c"]
151
152def ER(radius_equat_minor, radius_equat_major, radius_polar):
153    """
154    Returns the effective radius used in the S*P calculation
155    """
156    import numpy as np
157    from .ellipsoid import ER as ellipsoid_ER
158
159    # now that radii can be in any size order, radii need sorting a,b,c where a~b and c is either much smaller
160    # or much larger
161    radii = np.vstack((radius_equat_major, radius_equat_minor, radius_polar))
162    radii = np.sort(radii, axis=0)
163    selector = (radii[1] - radii[0]) > (radii[2] - radii[1])
164    polar = np.where(selector, radii[0], radii[2])
165    equatorial = np.sqrt(np.where(~selector, radii[0]*radii[1], radii[1]*radii[2]))
166    return ellipsoid_ER(polar, equatorial)
167
168demo = dict(scale=1, background=0,
169            sld=6, sld_solvent=1,
170            theta=30, phi=15, psi=5,
171            radius_equat_minor=25, radius_equat_major=36, radius_polar=50,
172            radius_equat_minor_pd=0, radius_equat_minor_pd_n=1,
173            radius_equat_major_pd=0, radius_equat_major_pd_n=1,
174            radius_polar_pd=.2, radius_polar_pd_n=30,
175            theta_pd=15, theta_pd_n=45,
176            phi_pd=15, phi_pd_n=1,
177            psi_pd=15, psi_pd_n=1)
178
179q = 0.1
180# april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct!
181# add 2d test after pull #890
182qx = q*cos(pi/6.0)
183qy = q*sin(pi/6.0)
184tests = [[{}, 0.05, 24.8839548033],
185#        [{'theta':80., 'phi':10.}, (qx, qy), 9999. ],
186        ]
187del qx, qy  # not necessary to delete, but cleaner
Note: See TracBrowser for help on using the repository browser.