source: sasmodels/sasmodels/models/ellipsoid.py @ 42356c8

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 42356c8 was 42356c8, checked in by Paul Kienzle <pkienzle@…>, 8 years ago

label all sld parameters

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