source: sasmodels/sasmodels/models/cylinder.py @ 31df0c9

core_shell_microgelscostrafo411magnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 31df0c9 was 31df0c9, checked in by Paul Kienzle <pkienzle@…>, 7 years ago

tuned random model generation for more models

  • Property mode set to 100644
File size: 6.8 KB
Line 
1# cylinder model
2# Note: model title and parameter table are inserted automatically
3r"""
4
5For information about polarised and magnetic scattering, see
6the :ref:`magnetism` documentation.
7
8Definition
9----------
10
11The output of the 2D scattering intensity function for oriented cylinders is
12given by (Guinier, 1955)
13
14.. math::
15
16    P(q,\alpha) = \frac{\text{scale}}{V} F^2(q,\alpha).sin(\alpha) + \text{background}
17
18where
19
20.. math::
21
22    F(q,\alpha) = 2 (\Delta \rho) V
23           \frac{\sin \left(\tfrac12 qL\cos\alpha \right)}
24                {\tfrac12 qL \cos \alpha}
25           \frac{J_1 \left(q R \sin \alpha\right)}{q R \sin \alpha}
26
27and $\alpha$ is the angle between the axis of the cylinder and $\vec q$, $V =\pi R^2L$
28is the volume of the cylinder, $L$ is the length of the cylinder, $R$ is the
29radius of the cylinder, and $\Delta\rho$ (contrast) is the scattering length
30density difference between the scatterer and the solvent. $J_1$ is the
31first order Bessel function.
32
33For randomly oriented particles:
34
35.. math::
36
37    F^2(q)=\int_{0}^{\pi/2}{F^2(q,\alpha)\sin(\alpha)d\alpha}=\int_{0}^{1}{F^2(q,u)du}
38
39
40Numerical integration is simplified by a change of variable to $u = cos(\alpha)$ with
41$sin(\alpha)=\sqrt{1-u^2}$.
42
43The output of the 1D scattering intensity function for randomly oriented
44cylinders is thus given by
45
46.. math::
47
48    P(q) = \frac{\text{scale}}{V}
49        \int_0^{\pi/2} F^2(q,\alpha) \sin \alpha\ d\alpha + \text{background}
50
51
52NB: The 2nd virial coefficient of the cylinder is calculated based on the
53radius and length values, and used as the effective radius for $S(q)$
54when $P(q) \cdot S(q)$ is applied.
55
56For oriented cylinders, we define the direction of the
57axis of the cylinder using two angles $\theta$ (note this is not the
58same as the scattering angle used in q) and $\phi$. Those angles
59are defined in :numref:`cylinder-angle-definition` .
60
61.. _cylinder-angle-definition:
62
63.. figure:: img/cylinder_angle_definition.png
64
65    Definition of the $\theta$ and $\phi$ orientation angles for a cylinder relative
66    to the beam line coordinates, plus an indication of their orientation distributions
67    which are described as rotations about each of the perpendicular axes $\delta_1$ and $\delta_2$
68    in the frame of the cylinder itself, which when $\theta = \phi = 0$ are parallel to the $Y$ and $X$ axes.
69
70.. figure:: img/cylinder_angle_projection.png
71
72    Examples for oriented cylinders.
73
74The $\theta$ and $\phi$ parameters to orient the cylinder only appear in the model when fitting 2d data.
75On introducing "Orientational Distribution" in the angles, "distribution of theta" and "distribution of phi" parameters will
76appear. These are actually rotations about the axes $\delta_1$ and $\delta_2$ of the cylinder, which when $\theta = \phi = 0$ are parallel
77to the $Y$ and $X$ axes of the instrument respectively. Some experimentation may be required to understand the 2d patterns fully.
78(Earlier implementations had numerical integration issues in some circumstances when orientation distributions passed through 90 degrees, such
79situations, with very broad distributions, should still be approached with care.)
80
81Validation
82----------
83
84Validation of the code was done by comparing the output of the 1D model
85to the output of the software provided by the NIST (Kline, 2006).
86The implementation of the intensity for fully oriented cylinders was done
87by averaging over a uniform distribution of orientations using
88
89.. math::
90
91    P(q) = \int_0^{\pi/2} d\phi
92        \int_0^\pi p(\theta) P_0(q,\theta) \sin \theta\ d\theta
93
94
95where $p(\theta,\phi) = 1$ is the probability distribution for the orientation
96and $P_0(q,\theta)$ is the scattering intensity for the fully oriented
97system, and then comparing to the 1D result.
98
99References
100----------
101
102J. S. Pedersen, Adv. Colloid Interface Sci. 70, 171-210 (1997).
103G. Fournet, Bull. Soc. Fr. Mineral. Cristallogr. 74, 39-113 (1951).
104"""
105
106import numpy as np  # type: ignore
107from numpy import pi, inf  # type: ignore
108
109name = "cylinder"
110title = "Right circular cylinder with uniform scattering length density."
111description = """
112     f(q,alpha) = 2*(sld - sld_solvent)*V*sin(qLcos(alpha)/2))
113                /[qLcos(alpha)/2]*J1(qRsin(alpha))/[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", "1e-6/Ang^2", 4, [-inf, inf], "sld",
130               "Cylinder scattering length density"],
131              ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "sld",
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, [-360, 360], "orientation",
138               "cylinder axis to beam angle"],
139              ["phi", "degrees",   60, [-360, 360], "orientation",
140               "rotation about beam"],
141             ]
142
143source = ["lib/polevl.c", "lib/sas_J1.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
152def random():
153    import numpy as np
154    V = 10**np.random.uniform(5, 12)
155    length = 10**np.random.uniform(-2, 2)*V**0.333
156    radius = np.sqrt(V/length/np.pi)
157    pars = dict(
158        #scale=1,
159        #background=0,
160        length=length,
161        radius=radius,
162    )
163    return pars
164
165
166# parameters for demo
167demo = dict(scale=1, background=0,
168            sld=6, sld_solvent=1,
169            radius=20, length=300,
170            theta=60, phi=60,
171            radius_pd=.2, radius_pd_n=9,
172            length_pd=.2, length_pd_n=10,
173            theta_pd=10, theta_pd_n=5,
174            phi_pd=10, phi_pd_n=5)
175
176qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5)
177# After redefinition of angles, find new tests values.  Was 10 10 in old coords
178tests = [[{}, 0.2, 0.042761386790780453],
179        [{}, [0.2], [0.042761386790780453]],
180#  new coords
181        [{'theta':80.1534480601659, 'phi':10.1510817110481}, (qx, qy), 0.03514647218513852],
182        [{'theta':80.1534480601659, 'phi':10.1510817110481}, [(qx, qy)], [0.03514647218513852]],
183# old coords   [{'theta':10.0, 'phi':10.0}, (qx, qy), 0.03514647218513852],
184#              [{'theta':10.0, 'phi':10.0}, [(qx, qy)], [0.03514647218513852]],
185        ]
186del qx, qy  # not necessary to delete, but cleaner
187# ADDED by:  RKH  ON: 18Mar2016 renamed sld's etc
Note: See TracBrowser for help on using the repository browser.