source: sasmodels/sasmodels/models/cylinder.py @ 50e1e40

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

use lib functions to simplify barbell, capped_cylinder, cylinder, ellipsoid, triaxial_ellipsoid; fix barbell calculation

  • Property mode set to 100644
File size: 5.9 KB
Line 
1# cylinder 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 cylinders is
10given by (Guinier, 1955)
11
12.. math::
13
14    P(q,\alpha) = \frac{\text{scale}}{V} F^2(q) + \text{background}
15
16where
17
18.. math::
19
20    F(q) = 2 (\Delta \rho) V
21           \frac{\sin \left(q\tfrac12 L\cos\alpha \right)}
22                {q\tfrac12 L \cos \alpha}
23           \frac{J_1 \left(q R \sin \alpha\right)}{q R \sin \alpha}
24
25and $\alpha$ is the angle between the axis of the cylinder and $\vec q$, $V$
26is the volume of the cylinder, $L$ is the length of the cylinder, $R$ is the
27radius of the cylinder, and $\Delta\rho$ (contrast) is the scattering length
28density difference between the scatterer and the solvent. $J_1$ is the
29first order Bessel function.
30
31To provide easy access to the orientation of the cylinder, we define the
32axis of the cylinder using two angles $\theta$ and $\phi$. Those angles
33are defined in :num:`figure #cylinder-orientation`.
34
35.. _cylinder-orientation:
36
37.. figure:: img/orientation.jpg
38
39    Definition of the angles for oriented cylinders.
40
41.. figure:: img/orientation2.jpg
42
43    Examples of the angles for oriented cylinders against the detector plane.
44
45NB: The 2nd virial coefficient of the cylinder is calculated based on the
46radius and length values, and used as the effective radius for $S(q)$
47when $P(q) \cdot S(q)$ is applied.
48
49The output of the 1D scattering intensity function for randomly oriented
50cylinders is then given by
51
52.. math::
53
54    P(q) = \frac{\text{scale}}{V}
55        \int_0^{\pi/2} F^2(q,\alpha) \sin \alpha\ d\alpha + \text{background}
56
57The $\theta$ and $\phi$ parameters are not used for the 1D output.
58
59Validation
60----------
61
62Validation of our code was done by comparing the output of the 1D model
63to the output of the software provided by the NIST (Kline, 2006).
64:num:`Figure #cylinder-compare` shows a comparison of
65the 1D output of our model and the output of the NIST software.
66
67.. _cylinder-compare:
68
69.. figure:: img/cylinder_compare.jpg
70
71    Comparison of the SasView scattering intensity for a cylinder with the
72    output of the NIST SANS analysis software.
73    The parameters were set to: *scale* = 1.0, *radius* = 20 |Ang|,
74    *length* = 400 |Ang|, *contrast* = 3e-6 |Ang^-2|, and
75    *background* = 0.01 |cm^-1|.
76
77In general, averaging over a distribution of orientations is done by
78evaluating the following
79
80.. math::
81
82    P(q) = \int_0^{\pi/2} d\phi
83        \int_0^\pi p(\theta, \phi) P_0(q,\alpha) \sin \theta\ d\theta
84
85
86where $p(\theta,\phi)$ is the probability distribution for the orientation
87and $P_0(q,\alpha)$ is the scattering intensity for the fully oriented
88system. Since we have no other software to compare the implementation of
89the intensity for fully oriented cylinders, we can compare the result of
90averaging our 2D output using a uniform distribution $p(\theta, \phi) = 1.0$.
91:num:`Figure #cylinder-crosscheck` shows the result of
92such a cross-check.
93
94.. _cylinder-crosscheck:
95
96.. figure:: img/cylinder_crosscheck.jpg
97
98    Comparison of the intensity for uniformly distributed cylinders
99    calculated from our 2D model and the intensity from the NIST SANS
100    analysis software.
101    The parameters used were: *scale* = 1.0, *radius* = 20 |Ang|,
102    *length* = 400 |Ang|, *contrast* = 3e-6 |Ang^-2|, and
103    *background* = 0.0 |cm^-1|.
104"""
105
106import numpy as np
107from numpy import pi, inf
108
109name = "cylinder"
110title = "Right circular cylinder with uniform scattering length density."
111description = """
112     f(q,alpha) = 2*(sld - solvent_sld)*V*sin(qLcos(alpha/2))
113                /[qLcos(alpha/2)]*J1(qRsin(alpha/2))/[qRsin(alpha)]
114
115            P(q,alpha)= scale/V*f(q,alpha)^(2)+background
116            V: Volume of the cylinder
117            R: Radius of the cylinder
118            L: Length of the cylinder
119            J1: The bessel function
120            alpha: angle between the axis of the
121            cylinder and the q-vector for 1D
122            :the ouput is P(q)=scale/V*integral
123            from pi/2 to zero of...
124            f(q,alpha)^(2)*sin(alpha)*dalpha + background
125"""
126category = "shape:cylinder"
127
128#             [ "name", "units", default, [lower, upper], "type", "description"],
129parameters = [["sld", "4e-6/Ang^2", 4, [-inf, inf], "",
130               "Cylinder scattering length density"],
131              ["solvent_sld", "1e-6/Ang^2", 1, [-inf, inf], "",
132               "Solvent scattering length density"],
133              ["radius", "Ang", 20, [0, inf], "volume",
134               "Cylinder radius"],
135              ["length", "Ang", 400, [0, inf], "volume",
136               "Cylinder length"],
137              ["theta", "degrees", 60, [-inf, inf], "orientation",
138               "In plane angle"],
139              ["phi", "degrees", 60, [-inf, inf], "orientation",
140               "Out of plane angle"],
141             ]
142
143source = ["lib/J1c.c", "lib/gauss76.c", "cylinder.c"]
144
145def ER(radius, length):
146    """
147        Return equivalent radius (ER)
148    """
149    ddd = 0.75 * radius * (2 * radius * length + (length + radius) * (length + pi * radius))
150    return 0.5 * (ddd) ** (1. / 3.)
151
152# parameters for demo
153demo = dict(scale=1, background=0,
154            sld=6, solvent_sld=1,
155            radius=20, length=300,
156            theta=60, phi=60,
157            radius_pd=.2, radius_pd_n=9,
158            length_pd=.2, length_pd_n=10,
159            theta_pd=10, theta_pd_n=5,
160            phi_pd=10, phi_pd_n=5)
161
162# For testing against the old sasview models, include the converted parameter
163# names and the target sasview model name.
164oldname = 'CylinderModel'
165oldpars = dict(theta='cyl_theta', phi='cyl_phi', sld='sldCyl', solvent_sld='sldSolv')
166
167
168qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5)
169tests = [[{}, 0.2, 0.041761386790780453],
170         [{}, [0.2], [0.041761386790780453]],
171         [{'theta':10.0, 'phi':10.0}, (qx, qy), 0.03414647218513852],
172         [{'theta':10.0, 'phi':10.0}, [(qx, qy)], [0.03414647218513852]],
173        ]
174del qx, qy  # not necessary to delete, but cleaner
Note: See TracBrowser for help on using the repository browser.