source: sasmodels/sasmodels/models/cylinder.py @ 3330bb4

core_shell_microgelscostrafo411magnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 3330bb4 was 3330bb4, checked in by ajj, 7 years ago

Normalise line endings

  • Property mode set to 100644
File size: 5.5 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.jpg
64
65    Definition of the angles for oriented cylinders.
66
67The $\theta$ and $\phi$ parameters only appear in the model when fitting 2d data.
68
69Validation
70----------
71
72Validation of the code was done by comparing the output of the 1D model
73to the output of the software provided by the NIST (Kline, 2006).
74The implementation of the intensity for fully oriented cylinders was done
75by averaging over a uniform distribution of orientations using
76
77.. math::
78
79    P(q) = \int_0^{\pi/2} d\phi
80        \int_0^\pi p(\theta) P_0(q,\theta) \sin \theta\ d\theta
81
82
83where $p(\theta,\phi) = 1$ is the probability distribution for the orientation
84and $P_0(q,\theta)$ is the scattering intensity for the fully oriented
85system, and then comparing to the 1D result.
86
87References
88----------
89
90J. S. Pedersen, Adv. Colloid Interface Sci. 70, 171-210 (1997).
91G. Fournet, Bull. Soc. Fr. Mineral. Cristallogr. 74, 39-113 (1951).
92"""
93
94import numpy as np  # type: ignore
95from numpy import pi, inf  # type: ignore
96
97name = "cylinder"
98title = "Right circular cylinder with uniform scattering length density."
99description = """
100     f(q,alpha) = 2*(sld - sld_solvent)*V*sin(qLcos(alpha)/2))
101                /[qLcos(alpha)/2]*J1(qRsin(alpha))/[qRsin(alpha)]
102
103            P(q,alpha)= scale/V*f(q,alpha)^(2)+background
104            V: Volume of the cylinder
105            R: Radius of the cylinder
106            L: Length of the cylinder
107            J1: The bessel function
108            alpha: angle between the axis of the
109            cylinder and the q-vector for 1D
110            :the ouput is P(q)=scale/V*integral
111            from pi/2 to zero of...
112            f(q,alpha)^(2)*sin(alpha)*dalpha + background
113"""
114category = "shape:cylinder"
115
116#             [ "name", "units", default, [lower, upper], "type", "description"],
117parameters = [["sld", "1e-6/Ang^2", 4, [-inf, inf], "sld",
118               "Cylinder scattering length density"],
119              ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "sld",
120               "Solvent scattering length density"],
121              ["radius", "Ang", 20, [0, inf], "volume",
122               "Cylinder radius"],
123              ["length", "Ang", 400, [0, inf], "volume",
124               "Cylinder length"],
125              ["theta", "degrees", 60, [-inf, inf], "orientation",
126               "latitude"],
127              ["phi", "degrees", 60, [-inf, inf], "orientation",
128               "longitude"],
129             ]
130
131source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c",  "cylinder.c"]
132
133def ER(radius, length):
134    """
135        Return equivalent radius (ER)
136    """
137    ddd = 0.75 * radius * (2 * radius * length + (length + radius) * (length + pi * radius))
138    return 0.5 * (ddd) ** (1. / 3.)
139
140# parameters for demo
141demo = dict(scale=1, background=0,
142            sld=6, sld_solvent=1,
143            radius=20, length=300,
144            theta=60, phi=60,
145            radius_pd=.2, radius_pd_n=9,
146            length_pd=.2, length_pd_n=10,
147            theta_pd=10, theta_pd_n=5,
148            phi_pd=10, phi_pd_n=5)
149
150qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5)
151# After redefinition of angles, find new tests values.  Was 10 10 in old coords
152tests = [[{}, 0.2, 0.042761386790780453],
153        [{}, [0.2], [0.042761386790780453]],
154#  new coords   
155        [{'theta':80.1534480601659, 'phi':10.1510817110481}, (qx, qy), 0.03514647218513852],
156        [{'theta':80.1534480601659, 'phi':10.1510817110481}, [(qx, qy)], [0.03514647218513852]],
157# old coords   [{'theta':10.0, 'phi':10.0}, (qx, qy), 0.03514647218513852],
158#              [{'theta':10.0, 'phi':10.0}, [(qx, qy)], [0.03514647218513852]],
159        ]
160del qx, qy  # not necessary to delete, but cleaner
161# ADDED by:  RKH  ON: 18Mar2016 renamed sld's etc
Note: See TracBrowser for help on using the repository browser.