source: sasmodels/sasmodels/models/_cylpy.py @ 59ee4db

Last change on this file since 59ee4db was 59ee4db, checked in by Paul Kienzle <pkienzle@…>, 6 years ago

replace sincos with sin,cos in python cylinder example

  • Property mode set to 100644
File size: 7.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 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
98J. S. Pedersen, Adv. Colloid Interface Sci. 70, 171-210 (1997).
99G. Fournet, Bull. Soc. Fr. Mineral. Cristallogr. 74, 39-113 (1951).
100"""
101
102import numpy as np  # type: ignore
103from sasmodels.special import sas_2J1x_x, sas_sinx_x, pi, inf, sincos, square, gauss76 as gauss
104from sasmodels.special import sin, cos
105
106name = "cylinder"
107title = "Right circular cylinder with uniform scattering length density."
108description = """
109     f(q,alpha) = 2*(sld - sld_solvent)*V*sin(qLcos(alpha)/2))
110                /[qLcos(alpha)/2]*J1(qRsin(alpha))/[qRsin(alpha)]
111
112            P(q,alpha)= scale/V*f(q,alpha)^(2)+background
113            V: Volume of the cylinder
114            R: Radius of the cylinder
115            L: Length of the cylinder
116            J1: The bessel function
117            alpha: angle between the axis of the
118            cylinder and the q-vector for 1D
119            :the ouput is P(q)=scale/V*integral
120            from pi/2 to zero of...
121            f(q,alpha)^(2)*sin(alpha)*dalpha + background
122"""
123category = "shape:cylinder"
124
125#             [ "name", "units", default, [lower, upper], "type", "description"],
126parameters = [["sld", "1e-6/Ang^2", 4, [-inf, inf], "sld",
127               "Cylinder scattering length density"],
128              ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "sld",
129               "Solvent scattering length density"],
130              ["radius", "Ang", 20, [0, inf], "volume",
131               "Cylinder radius"],
132              ["length", "Ang", 400, [0, inf], "volume",
133               "Cylinder length"],
134              ["theta", "degrees", 60, [-360, 360], "orientation",
135               "cylinder axis to beam angle"],
136              ["phi", "degrees", 60, [-360, 360], "orientation",
137               "rotation about beam"],
138             ]
139
140def invalid(v):
141    return v.radius < 0 or v.length < 0
142
143def form_volume(radius, length):
144    return pi*radius**2*length
145
146def fq(qab, qc, radius, length):
147    return sas_2J1x_x(qab*radius) * sas_sinx_x(qc*0.5*length)
148
149def orient_avg_1D(q, radius, length):
150    # translate a point in [-1,1] to a point in [0, pi/2]
151    zm = pi/4
152    zb = pi/4
153
154    total = 0.0
155    for i in range(gauss.n):
156        theta = gauss.z[i]*zm + zb
157        # TODO: should be using u-substitution of cos(theta) over [0, 1]
158        #sin_theta, cos_theta = sincos(theta)
159        sin_theta = sin(theta)
160        cos_theta = cos(theta)
161        form = fq(q*sin_theta, q*cos_theta, radius, length)
162        total += gauss.w[i] * form * form * sin_theta
163
164    # translate dx in [-1,1] to dx in [lower,upper]
165    return total*zm
166
167def Iq(q, sld, solvent_sld, radius, length):
168    s = (sld - solvent_sld) * form_volume(radius, length)
169    return 1.0e-4 * s * s * orient_avg_1D(q, radius, length)
170
171def Iqxy(qab, qc, sld, solvent_sld, radius, length):
172    s = (sld-solvent_sld) * form_volume(radius, length)
173    form = fq(qab, qc, radius, length)
174    return 1.0e-4 * square(s * form)
175
176def ER(radius, length):
177    """
178    Return equivalent radius (ER)
179    """
180    ddd = 0.75 * radius * (2 * radius * length + (length + radius) * (length + pi * radius))
181    return 0.5 * (ddd) ** (1. / 3.)
182
183def random():
184    volume = 10**np.random.uniform(5, 12)
185    length = 10**np.random.uniform(-2, 2)*volume**0.333
186    radius = np.sqrt(volume/length/np.pi)
187    pars = dict(
188        #scale=1,
189        #background=0,
190        length=length,
191        radius=radius,
192    )
193    return pars
194
195
196# parameters for demo
197demo = dict(scale=1, background=0,
198            sld=6, sld_solvent=1,
199            radius=20, length=300,
200            theta=60, phi=60,
201            radius_pd=.2, radius_pd_n=9,
202            length_pd=.2, length_pd_n=10,
203            theta_pd=10, theta_pd_n=5,
204            phi_pd=10, phi_pd_n=5)
205
206qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5)
207# After redefinition of angles, find new tests values.  Was 10 10 in old coords
208tests = [
209    [{}, 0.2, 0.042761386790780453],
210    [{}, [0.2], [0.042761386790780453]],
211    #  new coords
212    [{'theta':80.1534480601659, 'phi':10.1510817110481}, (qx, qy), 0.03514647218513852],
213    [{'theta':80.1534480601659, 'phi':10.1510817110481}, [(qx, qy)], [0.03514647218513852]],
214    # old coords
215    #[{'theta':10.0, 'phi':10.0}, (qx, qy), 0.03514647218513852],
216    #[{'theta':10.0, 'phi':10.0}, [(qx, qy)], [0.03514647218513852]],
217]
218del qx, qy  # not necessary to delete, but cleaner
219# ADDED by:  RKH  ON: 18Mar2016 renamed sld's etc
Note: See TracBrowser for help on using the repository browser.