source: sasmodels/sasmodels/models/cylinder.py @ 0507e09

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 0507e09 was 0507e09, checked in by smk78, 5 months ago

Added link to source code to each model. Closes #883

  • Property mode set to 100644
File size: 8.3 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 2d scattering from 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` , for further details see :ref:`orientation` .
60
61.. _cylinder-angle-definition:
62
63.. figure:: img/cylinder_angle_definition.png
64
65    Angles $\theta$ and $\phi$ orient the cylinder relative
66    to the beam line coordinates, where the beam is along the $z$ axis. Rotation $\theta$, initially
67    in the $xz$ plane, is carried out first, then rotation $\phi$ about the $z$ axis. Orientation distributions
68    are described as rotations about two perpendicular axes $\delta_1$ and $\delta_2$
69    in the frame of the cylinder itself, which when $\theta = \phi = 0$ are parallel to the $Y$ and $X$ axes.
70
71.. figure:: img/cylinder_angle_projection.png
72
73    Examples for oriented cylinders.
74
75The $\theta$ and $\phi$ parameters to orient the cylinder only appear in the model when fitting 2d data.
76
77Validation
78----------
79
80Validation of the code was done by comparing the output of the 1D model
81to the output of the software provided by the NIST (Kline, 2006).
82The implementation of the intensity for fully oriented cylinders was done
83by averaging over a uniform distribution of orientations using
84
85.. math::
86
87    P(q) = \int_0^{\pi/2} d\phi
88        \int_0^\pi p(\theta) P_0(q,\theta) \sin \theta\ d\theta
89
90
91where $p(\theta,\phi) = 1$ is the probability distribution for the orientation
92and $P_0(q,\theta)$ is the scattering intensity for the fully oriented
93system, and then comparing to the 1D result.
94
95References
96----------
97
98.. [#] J. Pedersen, *Adv. Colloid Interface Sci.*, 70 (1997) 171-210
99.. [#] G. Fournet, *Bull. Soc. Fr. Mineral. Cristallogr.*, 74 (1951) 39-113
100.. [#] L. Onsager, *Ann. New York Acad. Sci.*, 51 (1949) 627-659
101
102Source
103------
104
105`cylinder.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/cylinder.py>`_
106
107`cylinder.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/cylinder.c>`_
108
109Authorship and Verification
110----------------------------
111
112* **Author:**
113* **Last Modified by:**
114* **Last Reviewed by:**
115* **Source added by :** Steve King **Date:** March 25, 2019
116"""
117
118import numpy as np  # type: ignore
119from numpy import pi, inf  # type: ignore
120
121name = "cylinder"
122title = "Right circular cylinder with uniform scattering length density."
123description = """
124     f(q,alpha) = 2*(sld - sld_solvent)*V*sin(qLcos(alpha)/2))
125                /[qLcos(alpha)/2]*J1(qRsin(alpha))/[qRsin(alpha)]
126
127            P(q,alpha)= scale/V*f(q,alpha)^(2)+background
128            V: Volume of the cylinder
129            R: Radius of the cylinder
130            L: Length of the cylinder
131            J1: The bessel function
132            alpha: angle between the axis of the
133            cylinder and the q-vector for 1D
134            :the ouput is P(q)=scale/V*integral
135            from pi/2 to zero of...
136            f(q,alpha)^(2)*sin(alpha)*dalpha + background
137"""
138category = "shape:cylinder"
139
140#             [ "name", "units", default, [lower, upper], "type", "description"],
141parameters = [["sld", "1e-6/Ang^2", 4, [-inf, inf], "sld",
142               "Cylinder scattering length density"],
143              ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "sld",
144               "Solvent scattering length density"],
145              ["radius", "Ang", 20, [0, inf], "volume",
146               "Cylinder radius"],
147              ["length", "Ang", 400, [0, inf], "volume",
148               "Cylinder length"],
149              ["theta", "degrees", 60, [-360, 360], "orientation",
150               "cylinder axis to beam angle"],
151              ["phi", "degrees", 60, [-360, 360], "orientation",
152               "rotation about beam"],
153             ]
154
155source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "cylinder.c"]
156have_Fq = True
157effective_radius_type = [
158    "excluded volume", "equivalent volume sphere", "radius",
159    "half length", "half min dimension", "half max dimension", "half diagonal",
160    ]
161
162def random():
163    """Return a random parameter set for the model."""
164    volume = 10**np.random.uniform(5, 12)
165    length = 10**np.random.uniform(-2, 2)*volume**0.333
166    radius = np.sqrt(volume/length/np.pi)
167    pars = dict(
168        #scale=1,
169        #background=0,
170        length=length,
171        radius=radius,
172    )
173    return pars
174
175
176# parameters for demo
177demo = dict(scale=1, background=0,
178            sld=6, sld_solvent=1,
179            radius=20, length=300,
180            theta=60, phi=60,
181            radius_pd=.2, radius_pd_n=9,
182            length_pd=.2, length_pd_n=10,
183            theta_pd=10, theta_pd_n=5,
184            phi_pd=10, phi_pd_n=5)
185
186# pylint: disable=bad-whitespace, line-too-long
187qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5)
188# After redefinition of angles, find new tests values.  Was 10 10 in old coords
189tests = [
190    [{}, 0.2, 0.042761386790780453],
191    [{}, [0.2], [0.042761386790780453]],
192    #  new coords
193    [{'theta':80.1534480601659, 'phi':10.1510817110481}, (qx, qy), 0.03514647218513852],
194    [{'theta':80.1534480601659, 'phi':10.1510817110481}, [(qx, qy)], [0.03514647218513852]],
195    # old coords
196    #[{'theta':10.0, 'phi':10.0}, (qx, qy), 0.03514647218513852],
197    #[{'theta':10.0, 'phi':10.0}, [(qx, qy)], [0.03514647218513852]],
198]
199del qx, qy  # not necessary to delete, but cleaner
200
201# Default radius and length
202def calc_volume(radius, length):
203    """Return form volume for cylinder."""
204    return pi*radius**2*length
205def calc_r_effs(radius, length):
206    """Return effective radii for modes 0-7 of cylinder."""
207    return [
208        0.,
209        0.5*(0.75*radius*(2.0*radius*length
210                          + (radius + length)*(pi*radius + length)))**(1./3.),
211        (0.75*radius**2*length)**(1./3.),
212        radius,
213        length/2.,
214        min(radius, length/2.),
215        max(radius, length/2.),
216        np.sqrt(4*radius**2 + length**2)/2.,
217    ]
218r_effs = calc_r_effs(parameters[2][2], parameters[3][2])
219cyl_vol = calc_volume(parameters[2][2], parameters[3][2])
220tests.extend([
221    ({'radius_effective_mode': 0}, 0.1, None, None, r_effs[0], cyl_vol, 1.0),
222    ({'radius_effective_mode': 1}, 0.1, None, None, r_effs[1], None, None),
223    ({'radius_effective_mode': 2}, 0.1, None, None, r_effs[2], None, None),
224    ({'radius_effective_mode': 3}, 0.1, None, None, r_effs[3], None, None),
225    ({'radius_effective_mode': 4}, 0.1, None, None, r_effs[4], None, None),
226    ({'radius_effective_mode': 5}, 0.1, None, None, r_effs[5], None, None),
227    ({'radius_effective_mode': 6}, 0.1, None, None, r_effs[6], None, None),
228    ({'radius_effective_mode': 7}, 0.1, None, None, r_effs[7], None, None),
229])
230del r_effs, cyl_vol
231# pylint: enable=bad-whitespace, line-too-long
232
233# ADDED by:  RKH  ON: 18Mar2016 renamed sld's etc
Note: See TracBrowser for help on using the repository browser.