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

ticket-1257-vesicle-productticket_1156ticket_822_more_unit_tests
Last change on this file since a34b811 was a34b811, checked in by Paul Kienzle <pkienzle@…>, 6 months ago

use radius_effective/radius_effective_mode/radius_effective_modes consistently throughout the code

  • Property mode set to 100644
File size: 6.8 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
70Validation
71----------
72
73Validation of the code was done by comparing the output of the 1D model
74to the output of the software provided by the NIST (Kline, 2006).
75
76The implementation of the intensity for fully oriented ellipsoids was
77validated by averaging the 2D output using a uniform distribution
78$p(\theta,\phi) = 1.0$ and comparing with the output of the 1D calculation.
79
80
81.. _ellipsoid-comparison-2d:
82
83.. figure:: img/ellipsoid_comparison_2d.jpg
84
85    Comparison of the intensity for uniformly distributed ellipsoids
86    calculated from our 2D model and the intensity from the NIST SANS
87    analysis software. The parameters used were: *scale* = 1.0,
88    *radius_polar* = 20 |Ang|, *radius_equatorial* = 400 |Ang|,
89    *contrast* = 3e-6 |Ang^-2|, and *background* = 0.0 |cm^-1|.
90
91The discrepancy above $q$ = 0.3 |cm^-1| is due to the way the form factors
92are calculated in the c-library provided by NIST. A numerical integration
93has to be performed to obtain $P(q)$ for randomly oriented particles.
94The NIST software performs that integration with a 76-point Gaussian
95quadrature rule, which will become imprecise at high $q$ where the amplitude
96varies quickly as a function of $q$. The SasView result shown has been
97obtained by summing over 501 equidistant points. Our result was found
98to be stable over the range of $q$ shown for a number of points higher
99than 500.
100
101Model was also tested against the triaxial ellipsoid model with equal major
102and minor equatorial radii.  It is also consistent with the cyclinder model
103with polar radius equal to length and equatorial radius equal to radius.
104
105References
106----------
107
108.. [#] L A Feigin and D I Svergun. *Structure Analysis by Small-Angle X-Ray and Neutron Scattering*, Plenum Press, New York, 1987
109.. [#] A. Isihara. *J. Chem. Phys.*, 18 (1950) 1446-1449
110
111Source
112------
113
114`ellipsoid.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/ellipsoid.py>`_
115
116`ellipsoid.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/ellipsoid.c>`_
117
118Authorship and Verification
119----------------------------
120
121* **Author:** NIST IGOR/DANSE **Date:** pre 2010
122* **Converted to sasmodels by:** Helen Park **Date:** July 9, 2014
123* **Last Modified by:** Paul Kienzle **Date:** March 22, 2017
124* **Source added by :** Steve King **Date:** March 25, 2019
125"""
126from __future__ import division
127
128import numpy as np
129from numpy import inf, sin, cos, pi
130
131name = "ellipsoid"
132title = "Ellipsoid of revolution with uniform scattering length density."
133
134description = """\
135P(q.alpha)= scale*f(q)^2 + background, where f(q)= 3*(sld
136        - sld_solvent)*V*[sin(q*r(Rp,Re,alpha))
137        -q*r*cos(qr(Rp,Re,alpha))]
138        /[qr(Rp,Re,alpha)]^3"
139
140     r(Rp,Re,alpha)= [Re^(2)*(sin(alpha))^2
141        + Rp^(2)*(cos(alpha))^2]^(1/2)
142
143        sld: SLD of the ellipsoid
144        sld_solvent: SLD of the solvent
145        V: volume of the ellipsoid
146        Rp: polar radius of the ellipsoid
147        Re: equatorial radius of the ellipsoid
148"""
149category = "shape:ellipsoid"
150
151#             ["name", "units", default, [lower, upper], "type","description"],
152parameters = [["sld", "1e-6/Ang^2", 4, [-inf, inf], "sld",
153               "Ellipsoid scattering length density"],
154              ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "sld",
155               "Solvent scattering length density"],
156              ["radius_polar", "Ang", 20, [0, inf], "volume",
157               "Polar radius"],
158              ["radius_equatorial", "Ang", 400, [0, inf], "volume",
159               "Equatorial radius"],
160              ["theta", "degrees", 60, [-360, 360], "orientation",
161               "ellipsoid axis to beam angle"],
162              ["phi", "degrees", 60, [-360, 360], "orientation",
163               "rotation about beam"],
164             ]
165
166
167source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "ellipsoid.c"]
168have_Fq = True
169radius_effective_modes = [
170    "average curvature", "equivalent volume sphere", "min radius", "max radius",
171    ]
172
173def random():
174    """Return a random parameter set for the model."""
175    volume = 10**np.random.uniform(5, 12)
176    radius_polar = 10**np.random.uniform(1.3, 4)
177    radius_equatorial = np.sqrt(volume/radius_polar) # ignore 4/3 pi
178    pars = dict(
179        #background=0, sld=0, sld_solvent=1,
180        radius_polar=radius_polar,
181        radius_equatorial=radius_equatorial,
182    )
183    return pars
184
185demo = dict(scale=1, background=0,
186            sld=6, sld_solvent=1,
187            radius_polar=50, radius_equatorial=30,
188            theta=30, phi=15,
189            radius_polar_pd=.2, radius_polar_pd_n=15,
190            radius_equatorial_pd=.2, radius_equatorial_pd_n=15,
191            theta_pd=15, theta_pd_n=45,
192            phi_pd=15, phi_pd_n=1)
193q = 0.1
194# april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct!
195qx = q*cos(pi/6.0)
196qy = q*sin(pi/6.0)
197tests = [
198    [{}, 0.05, 54.8525847025],
199    [{'theta':80., 'phi':10.}, (qx, qy), 1.74134670026],
200]
201del qx, qy  # not necessary to delete, but cleaner
Note: See TracBrowser for help on using the repository browser.