source: sasmodels/sasmodels/models/ellipsoid.py @ c88f983

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since c88f983 was c88f983, checked in by Torin Cooper-Bennun <torin.cooper-bennun@…>, 6 years ago

Merge branch 'master' into beta_approx

  • Property mode set to 100644
File size: 7.4 KB
Line 
1# ellipsoid model
2# Note: model title and parameter table are inserted automatically
3r"""
4The form factor is normalized by the particle volume
5
6Definition
7----------
8
9The output of the 2D scattering intensity function for oriented ellipsoids
10is given by (Feigin, 1987)
11
12.. math::
13
14    P(q,\alpha) = \frac{\text{scale}}{V} F^2(q,\alpha) + \text{background}
15
16where
17
18.. math::
19
20    F(q,\alpha) = \Delta \rho V \frac{3(\sin qr  - qr \cos qr)}{(qr)^3}
21
22for
23
24.. math::
25
26    r = \left[ R_e^2 \sin^2 \alpha + R_p^2 \cos^2 \alpha \right]^{1/2}
27
28
29$\alpha$ is the angle between the axis of the ellipsoid and $\vec q$,
30$V = (4/3)\pi R_pR_e^2$ is the volume of the ellipsoid, $R_p$ is the polar
31radius along the rotational axis of the ellipsoid, $R_e$ is the equatorial
32radius perpendicular to the rotational axis of the ellipsoid and
33$\Delta \rho$ (contrast) is the scattering length density difference between
34the scatterer and the solvent.
35
36For randomly oriented particles use the orientational average,
37
38.. math::
39
40   \langle F^2(q) \rangle = \int_{0}^{\pi/2}{F^2(q,\alpha)\sin(\alpha)\,d\alpha}
41
42
43computed via substitution of $u=\sin(\alpha)$, $du=\cos(\alpha)\,d\alpha$ as
44
45.. math::
46
47    \langle F^2(q) \rangle = \int_0^1{F^2(q, u)\,du}
48
49with
50
51.. math::
52
53    r = R_e \left[ 1 + u^2\left(R_p^2/R_e^2 - 1\right)\right]^{1/2}
54
55For 2d data from oriented ellipsoids the direction of the rotation axis of
56the ellipsoid is defined using two angles $\theta$ and $\phi$ as for the
57:ref:`cylinder orientation figure <cylinder-angle-definition>`.
58For the ellipsoid, $\theta$ is the angle between the rotational axis
59and the $z$ -axis in the $xz$ plane followed by a rotation by $\phi$
60in the $xy$ plane, for further details of the calculation and angular
61dispersions see :ref:`orientation` .
62
63NB: The 2nd virial coefficient of the solid ellipsoid is calculated based
64on the $R_p$ and $R_e$ values, and used as the effective radius for
65$S(q)$ when $P(q) \cdot S(q)$ is applied.
66
67
68The $\theta$ and $\phi$ parameters are not used for the 1D output.
69
70
71
72Validation
73----------
74
75Validation of the code was done by comparing the output of the 1D model
76to the output of the software provided by the NIST (Kline, 2006).
77
78The implementation of the intensity for fully oriented ellipsoids was
79validated by averaging the 2D output using a uniform distribution
80$p(\theta,\phi) = 1.0$ and comparing with the output of the 1D calculation.
81
82
83.. _ellipsoid-comparison-2d:
84
85.. figure:: img/ellipsoid_comparison_2d.jpg
86
87    Comparison of the intensity for uniformly distributed ellipsoids
88    calculated from our 2D model and the intensity from the NIST SANS
89    analysis software. The parameters used were: *scale* = 1.0,
90    *radius_polar* = 20 |Ang|, *radius_equatorial* = 400 |Ang|,
91    *contrast* = 3e-6 |Ang^-2|, and *background* = 0.0 |cm^-1|.
92
93The discrepancy above $q$ = 0.3 |cm^-1| is due to the way the form factors
94are calculated in the c-library provided by NIST. A numerical integration
95has to be performed to obtain $P(q)$ for randomly oriented particles.
96The NIST software performs that integration with a 76-point Gaussian
97quadrature rule, which will become imprecise at high $q$ where the amplitude
98varies quickly as a function of $q$. The SasView result shown has been
99obtained by summing over 501 equidistant points. Our result was found
100to be stable over the range of $q$ shown for a number of points higher
101than 500.
102
103Model was also tested against the triaxial ellipsoid model with equal major
104and minor equatorial radii.  It is also consistent with the cyclinder model
105with polar radius equal to length and equatorial radius equal to radius.
106
107References
108----------
109
110L A Feigin and D I Svergun.
111*Structure Analysis by Small-Angle X-Ray and Neutron Scattering*,
112Plenum Press, New York, 1987.
113
114A. Isihara. J. Chem. Phys. 18(1950) 1446-1449
115
116Authorship and Verification
117----------------------------
118
119* **Author:** NIST IGOR/DANSE **Date:** pre 2010
120* **Converted to sasmodels by:** Helen Park **Date:** July 9, 2014
121* **Last Modified by:** Paul Kienzle **Date:** March 22, 2017
122"""
123from __future__ import division
124
125import numpy as np
126from numpy import inf, sin, cos, pi
127
128try:
129    from numpy import cbrt
130except ImportError:
131    def cbrt(x): return x ** (1.0/3.0)
132
133name = "ellipsoid"
134title = "Ellipsoid of revolution with uniform scattering length density."
135
136description = """\
137P(q.alpha)= scale*f(q)^2 + background, where f(q)= 3*(sld
138        - sld_solvent)*V*[sin(q*r(Rp,Re,alpha))
139        -q*r*cos(qr(Rp,Re,alpha))]
140        /[qr(Rp,Re,alpha)]^3"
141
142     r(Rp,Re,alpha)= [Re^(2)*(sin(alpha))^2
143        + Rp^(2)*(cos(alpha))^2]^(1/2)
144
145        sld: SLD of the ellipsoid
146        sld_solvent: SLD of the solvent
147        V: volume of the ellipsoid
148        Rp: polar radius of the ellipsoid
149        Re: equatorial radius of the ellipsoid
150"""
151category = "shape:ellipsoid"
152
153#             ["name", "units", default, [lower, upper], "type","description"],
154parameters = [["sld", "1e-6/Ang^2", 4, [-inf, inf], "sld",
155               "Ellipsoid scattering length density"],
156              ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "sld",
157               "Solvent scattering length density"],
158              ["radius_polar", "Ang", 20, [0, inf], "volume",
159               "Polar radius"],
160              ["radius_equatorial", "Ang", 400, [0, inf], "volume",
161               "Equatorial radius"],
162              ["theta", "degrees", 60, [-360, 360], "orientation",
163               "ellipsoid axis to beam angle"],
164              ["phi", "degrees", 60, [-360, 360], "orientation",
165               "rotation about beam"],
166             ]
167
168
169source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "ellipsoid.c"]
170have_Fq = True
171
172def ER(radius_polar, radius_equatorial):
173    # see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449
174    ee = np.empty_like(radius_polar)
175    idx = radius_polar > radius_equatorial
176    ee[idx] = (radius_polar[idx] ** 2 - radius_equatorial[idx] ** 2) / radius_polar[idx] ** 2
177    idx = radius_polar < radius_equatorial
178    ee[idx] = (radius_equatorial[idx] ** 2 - radius_polar[idx] ** 2) / radius_equatorial[idx] ** 2
179    valid = (radius_polar * radius_equatorial != 0) & (radius_polar != radius_equatorial)
180    bd = 1.0 - ee[valid]
181    e1 = np.sqrt(ee[valid])
182    b1 = 1.0 + np.arcsin(e1) / (e1 * np.sqrt(bd))
183    bL = (1.0 + e1) / (1.0 - e1)
184    b2 = 1.0 + bd / 2 / e1 * np.log(bL)
185    delta = 0.75 * b1 * b2
186    ddd = 2.0 * (delta + 1.0) * (radius_polar * radius_equatorial**2)[valid]
187
188    r = np.zeros_like(radius_polar)
189    r[valid] = 0.5 * cbrt(ddd)
190    idx = radius_polar == radius_equatorial
191    r[idx] = radius_polar[idx]
192    return r
193
194def random():
195    volume = 10**np.random.uniform(5, 12)
196    radius_polar = 10**np.random.uniform(1.3, 4)
197    radius_equatorial = np.sqrt(volume/radius_polar) # ignore 4/3 pi
198    pars = dict(
199        #background=0, sld=0, sld_solvent=1,
200        radius_polar=radius_polar,
201        radius_equatorial=radius_equatorial,
202    )
203    return pars
204
205demo = dict(scale=1, background=0,
206            sld=6, sld_solvent=1,
207            radius_polar=50, radius_equatorial=30,
208            theta=30, phi=15,
209            radius_polar_pd=.2, radius_polar_pd_n=15,
210            radius_equatorial_pd=.2, radius_equatorial_pd_n=15,
211            theta_pd=15, theta_pd_n=45,
212            phi_pd=15, phi_pd_n=1)
213q = 0.1
214# april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct!
215qx = q*cos(pi/6.0)
216qy = q*sin(pi/6.0)
217tests = [
218    [{}, 0.05, 54.8525847025],
219    [{'theta':80., 'phi':10.}, (qx, qy), 1.74134670026],
220]
221del qx, qy  # not necessary to delete, but cleaner
Note: See TracBrowser for help on using the repository browser.