source: sasmodels/sasmodels/models/ellipsoid.py @ 041bc75

core_shell_microgelscostrafo411magnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 041bc75 was 416f5c7, checked in by richardh, 7 years ago

fixes for numref warnings in docu, new equations core_shell_bicelle core_shell_ellipsoid

  • Property mode set to 100644
File size: 5.9 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) = \frac{3 \Delta \rho V (\sin[qr(R_p,R_e,\alpha)]
21                - \cos[qr(R_p,R_e,\alpha)])}
22                {[qr(R_p,R_e,\alpha)]^3}
23
24and
25
26.. math::
27
28    r(R_p,R_e,\alpha) = \left[ R_e^2 \sin^2 \alpha
29        + R_p^2 \cos^2 \alpha \right]^{1/2}
30
31
32$\alpha$ is the angle between the axis of the ellipsoid and $\vec q$,
33$V = (4/3)\pi R_pR_e^2$ is the volume of the ellipsoid , $R_p$ is the polar radius along the
34rotational axis of the ellipsoid, $R_e$ is the equatorial radius perpendicular
35to the rotational axis of the ellipsoid and $\Delta \rho$ (contrast) is the
36scattering length density difference between the scatterer and the solvent.
37
38For randomly oriented particles:
39
40.. math::
41
42   F^2(q)=\int_{0}^{\pi/2}{F^2(q,\alpha)\sin(\alpha)d\alpha}
43
44
45To provide easy access to the orientation of the ellipsoid, we define
46the rotation axis of the ellipsoid using two angles $\theta$ and $\phi$.
47These angles are defined in the
48:ref:`cylinder orientation figure <cylinder-angle-definition>`.
49For the ellipsoid, $\theta$ is the angle between the rotational axis
50and the $z$ -axis.
51
52NB: The 2nd virial coefficient of the solid ellipsoid is calculated based
53on the $R_p$ and $R_e$ values, and used as the effective radius for
54$S(q)$ when $P(q) \cdot S(q)$ is applied.
55
56
57The $\theta$ and $\phi$ parameters are not used for the 1D output.
58
59.. _ellipsoid-geometry:
60
61.. figure:: img/ellipsoid_angle_projection.jpg
62
63    The angles for oriented ellipsoid, shown here as oblate, $a$ = $R_p$ and $b$ = $R_e$
64
65Validation
66----------
67
68Validation of the code was done by comparing the output of the 1D model
69to the output of the software provided by the NIST (Kline, 2006).
70
71The implementation of the intensity for fully oriented ellipsoids was
72validated by averaging the 2D output using a uniform distribution
73$p(\theta,\phi) = 1.0$ and comparing with the output of the 1D calculation.
74
75
76.. _ellipsoid-comparison-2d:
77
78.. figure:: img/ellipsoid_comparison_2d.jpg
79
80    Comparison of the intensity for uniformly distributed ellipsoids
81    calculated from our 2D model and the intensity from the NIST SANS
82    analysis software. The parameters used were: *scale* = 1.0,
83    *radius_polar* = 20 |Ang|, *radius_equatorial* = 400 |Ang|,
84    *contrast* = 3e-6 |Ang^-2|, and *background* = 0.0 |cm^-1|.
85
86The discrepancy above $q$ = 0.3 |cm^-1| is due to the way the form factors
87are calculated in the c-library provided by NIST. A numerical integration
88has to be performed to obtain $P(q)$ for randomly oriented particles.
89The NIST software performs that integration with a 76-point Gaussian
90quadrature rule, which will become imprecise at high $q$ where the amplitude
91varies quickly as a function of $q$. The SasView result shown has been
92obtained by summing over 501 equidistant points. Our result was found
93to be stable over the range of $q$ shown for a number of points higher
94than 500.
95
96References
97----------
98
99L A Feigin and D I Svergun.
100*Structure Analysis by Small-Angle X-Ray and Neutron Scattering*,
101Plenum Press, New York, 1987.
102"""
103
104from numpy import inf
105
106name = "ellipsoid"
107title = "Ellipsoid of revolution with uniform scattering length density."
108
109description = """\
110P(q.alpha)= scale*f(q)^2 + background, where f(q)= 3*(sld
111        - sld_solvent)*V*[sin(q*r(Rp,Re,alpha))
112        -q*r*cos(qr(Rp,Re,alpha))]
113        /[qr(Rp,Re,alpha)]^3"
114
115     r(Rp,Re,alpha)= [Re^(2)*(sin(alpha))^2
116        + Rp^(2)*(cos(alpha))^2]^(1/2)
117
118        sld: SLD of the ellipsoid
119        sld_solvent: SLD of the solvent
120        V: volume of the ellipsoid
121        Rp: polar radius of the ellipsoid
122        Re: equatorial radius of the ellipsoid
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_polar", "Ang", 20, [0, inf], "volume",
132               "Polar radius"],
133              ["radius_equatorial", "Ang", 400, [0, inf], "volume",
134               "Equatorial radius"],
135              ["theta", "degrees", 60, [-inf, inf], "orientation",
136               "In plane angle"],
137              ["phi", "degrees", 60, [-inf, inf], "orientation",
138               "Out of plane angle"],
139             ]
140
141source = ["lib/sph_j1c.c", "lib/gauss76.c", "ellipsoid.c"]
142
143def ER(radius_polar, radius_equatorial):
144    import numpy as np
145
146    ee = np.empty_like(radius_polar)
147    idx = radius_polar > radius_equatorial
148    ee[idx] = (radius_polar[idx] ** 2 - radius_equatorial[idx] ** 2) / radius_polar[idx] ** 2
149    idx = radius_polar < radius_equatorial
150    ee[idx] = (radius_equatorial[idx] ** 2 - radius_polar[idx] ** 2) / radius_equatorial[idx] ** 2
151    idx = radius_polar == radius_equatorial
152    ee[idx] = 2 * radius_polar[idx]
153    valid = (radius_polar * radius_equatorial != 0)
154    bd = 1.0 - ee[valid]
155    e1 = np.sqrt(ee[valid])
156    b1 = 1.0 + np.arcsin(e1) / (e1 * np.sqrt(bd))
157    bL = (1.0 + e1) / (1.0 - e1)
158    b2 = 1.0 + bd / 2 / e1 * np.log(bL)
159    delta = 0.75 * b1 * b2
160
161    ddd = np.zeros_like(radius_polar)
162    ddd[valid] = 2.0 * (delta + 1.0) * radius_polar * radius_equatorial ** 2
163    return 0.5 * ddd ** (1.0 / 3.0)
164
165
166demo = dict(scale=1, background=0,
167            sld=6, sld_solvent=1,
168            radius_polar=50, radius_equatorial=30,
169            theta=30, phi=15,
170            radius_polar_pd=.2, radius_polar_pd_n=15,
171            radius_equatorial_pd=.2, radius_equatorial_pd_n=15,
172            theta_pd=15, theta_pd_n=45,
173            phi_pd=15, phi_pd_n=1)
Note: See TracBrowser for help on using the repository browser.