Skip to content

Commit c0eb28e

Browse files
author
Paul Kienzle
committed
adaptive integration for cylinder using qr heuristic
1 parent 87b504a commit c0eb28e

File tree

3 files changed

+11568
-0
lines changed

3 files changed

+11568
-0
lines changed

‎sasmodels/models/cylinderp.c‎

Lines changed: 104 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,104 @@
1+
static double
2+
form_volume(double radius, double length)
3+
{
4+
return M_PI*radius*radius*length;
5+
}
6+
7+
static double
8+
_fq(double qab, double qc, double radius, double length)
9+
{
10+
return sas_2J1x_x(qab*radius) * sas_sinx_x(qc*0.5*length);
11+
}
12+
13+
static double
14+
radius_from_excluded_volume(double radius, double length)
15+
{
16+
return 0.5*cbrt(0.75*radius*(2.0*radius*length
17+
+ (radius + length)*(M_PI*radius + length)));
18+
}
19+
20+
static double
21+
radius_from_volume(double radius, double length)
22+
{
23+
return cbrt(M_PI*radius*radius*length/M_4PI_3);
24+
}
25+
26+
static double
27+
radius_from_diagonal(double radius, double length)
28+
{
29+
return sqrt(radius*radius + 0.25*length*length);
30+
}
31+
32+
static double
33+
radius_effective(int mode, double radius, double length)
34+
{
35+
switch (mode) {
36+
default:
37+
case 1:
38+
return radius_from_excluded_volume(radius, length);
39+
case 2:
40+
return radius_from_volume(radius, length);
41+
case 3:
42+
return radius;
43+
case 4:
44+
return 0.5*length;
45+
case 5:
46+
return (radius < 0.5*length ? radius : 0.5*length);
47+
case 6:
48+
return (radius > 0.5*length ? radius : 0.5*length);
49+
case 7:
50+
return radius_from_diagonal(radius,length);
51+
}
52+
}
53+
54+
55+
static void
56+
Fq(double q,
57+
double *F1,
58+
double *F2,
59+
double sld,
60+
double solvent_sld,
61+
double radius,
62+
double length)
63+
{
64+
// translate a point in [-1,1] to a point in [0, pi/2]
65+
const double zm = M_PI_4;
66+
const double zb = M_PI_4;
67+
68+
double total_F1 = 0.0;
69+
double total_F2 = 0.0;
70+
71+
double qr = q*(radius > length ? radius : length);
72+
double *w, *z;
73+
int n = gauss_weights(qr, &w, &z);
74+
for (int i=0; i<n ;i++) {
75+
const double theta = z[i]*zm + zb;
76+
double sin_theta, cos_theta; // slots to hold sincos function output
77+
// theta (theta,phi) the projection of the cylinder on the detector plane
78+
SINCOS(theta , sin_theta, cos_theta);
79+
const double form = _fq(q*sin_theta, q*cos_theta, radius, length);
80+
total_F1 += w[i] * form * sin_theta;
81+
total_F2 += w[i] * form * form * sin_theta;
82+
}
83+
// translate dx in [-1,1] to dx in [lower,upper]
84+
total_F1 *= zm;
85+
total_F2 *= zm;
86+
const double s = (sld - solvent_sld) * form_volume(radius, length);
87+
*F1 = 1e-2 * s * total_F1;
88+
*F2 = 1e-4 * s * s * total_F2;
89+
}
90+
91+
92+
93+
static double
94+
Iqac(double qab, double qc,
95+
double sld,
96+
double solvent_sld,
97+
double radius,
98+
double length)
99+
{
100+
const double form = _fq(qab, qc, radius, length);
101+
const double s = (sld-solvent_sld) * form_volume(radius, length);
102+
return 1.0e-4 * square(s * form);
103+
}
104+

‎sasmodels/models/cylinderp.py‎

Lines changed: 251 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,251 @@
1+
# cylinder model
2+
# Note: model title and parameter table are inserted automatically
3+
r"""
4+
5+
For information about polarised and magnetic scattering, see
6+
the :ref:`magnetism` documentation.
7+
8+
Definition
9+
----------
10+
11+
The output of the 2D scattering intensity function for oriented cylinders is
12+
given by (Guinier, 1955)
13+
14+
.. math::
15+
16+
I(q,\alpha) = \frac{\text{scale}}{V} F^2(q,\alpha) + \text{background}
17+
18+
where
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+
27+
and $\alpha$ is the angle between the axis of the cylinder and $\vec q$,
28+
$V =\pi R^2L$ is the volume of the cylinder, $L$ is the length of the cylinder,
29+
$R$ is the radius of the cylinder, and $\Delta\rho$ (contrast) is the
30+
scattering length density difference between the scatterer and the solvent.
31+
$J_1$ is the first order Bessel function.
32+
33+
For randomly oriented particles:
34+
35+
.. math::
36+
37+
P(q)=F^2(q)=\int_{0}^{\pi/2}{F^2(q,\alpha)\sin(\alpha)d\alpha}
38+
39+
40+
The output of the 1D scattering intensity function for randomly oriented
41+
cylinders is thus given by
42+
43+
.. math::
44+
45+
I(q) = \frac{\text{scale}}{V}
46+
\int_0^{\pi/2} F^2(q,\alpha) \sin \alpha\ d\alpha + \text{background}
47+
48+
49+
NB: The 2nd virial coefficient of the cylinder is calculated based on the
50+
radius and length values, and used as the effective radius for $S(q)$
51+
when $P(q) \cdot S(q)$ is applied.
52+
53+
For 2d scattering from oriented cylinders, we define the direction of the
54+
axis of the cylinder using two angles $\theta$ (note this is not the same as
55+
the scattering angle used in q) and $\phi$. Those angles are defined in
56+
:numref:`cylinder-angle-definition` , for further details see
57+
:ref:`orientation`.
58+
59+
.. _cylinder-angle-definition:
60+
61+
.. figure:: img/cylinder_angle_definition.png
62+
63+
Angles $\theta$ and $\phi$ orient the cylinder relative to the beam line
64+
coordinates, where the beam is along the $z$ axis. Rotation $\theta$,
65+
initially in the $xz$ plane, is carried out first, then rotation $\phi$
66+
about the $z$ axis. Orientation distributions are described as rotations
67+
about two perpendicular axes $\delta_1$ and $\delta_2$ in the frame of
68+
the cylinder itself, which when $\theta = \phi = 0$ are parallel to the
69+
$Y$ and $X$ axes.
70+
71+
.. figure:: img/cylinder_angle_projection.png
72+
73+
Examples for oriented cylinders.
74+
75+
The $\theta$ and $\phi$ parameters to orient the cylinder only appear in the
76+
model when fitting 2d data.
77+
78+
Validation
79+
----------
80+
81+
Validation of the code was done by comparing the output of the 1D model
82+
to the output of the software provided by the NIST (Kline, 2006).
83+
The implementation of the intensity for fully oriented cylinders was done
84+
by averaging over a uniform distribution of orientations using
85+
86+
.. math::
87+
88+
P(q) = \int_0^{\pi/2} d\phi
89+
\int_0^\pi p(\theta) P_0(q,\theta) \sin \theta\ d\theta
90+
91+
92+
where $p(\theta,\phi) = 1$ is the probability distribution for the orientation
93+
and $P_0(q,\theta)$ is the scattering intensity for the fully oriented
94+
system, and then comparing to the 1D result.
95+
96+
References
97+
----------
98+
99+
#. J. Pedersen, *Adv. Colloid Interface Sci.*, 70 (1997) 171-210
100+
#. G. Fournet, *Bull. Soc. Fr. Mineral. Cristallogr.*, 74 (1951) 39-113
101+
#. L. Onsager, *Ann. New York Acad. Sci.*, 51 (1949) 627-659
102+
103+
Authorship and Verification
104+
----------------------------
105+
106+
* **Author:**
107+
* **Last Modified by:** Paul Butler (docs only) November 10, 2022
108+
* **Last Reviewed by:**
109+
"""
110+
111+
import numpy as np # type: ignore
112+
from numpy import pi, inf # type: ignore
113+
114+
name = "cylinder"
115+
title = "Right circular cylinder with uniform scattering length density."
116+
description = """
117+
f(q,alpha) = 2*(sld - sld_solvent)*V*sin(qLcos(alpha)/2))
118+
/[qLcos(alpha)/2]*J1(qRsin(alpha))/[qRsin(alpha)]
119+
120+
P(q,alpha)= scale/V*f(q,alpha)^(2)+background
121+
V: Volume of the cylinder
122+
R: Radius of the cylinder
123+
L: Length of the cylinder
124+
J1: The Bessel function
125+
alpha: angle between the axis of the
126+
cylinder and the q-vector for 1D
127+
:the ouput is P(q)=scale/V*integral
128+
from pi/2 to zero of...
129+
f(q,alpha)^(2)*sin(alpha)*dalpha + background
130+
"""
131+
category = "shape:cylinder"
132+
133+
# pylint: disable=bad-whitespace, line-too-long
134+
# [ "name", "units", default, [lower, upper], "type", "description"],
135+
parameters = [
136+
["sld", "1e-6/Ang^2", 4, [-inf, inf], "sld", "Cylinder scattering length density"],
137+
["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "sld", "Solvent scattering length density"],
138+
["radius", "Ang", 20, [0, inf], "volume", "Cylinder radius"],
139+
["length", "Ang", 400, [0, inf], "volume", "Cylinder length"],
140+
["theta", "degrees", 60, [-360, 360], "orientation", "cylinder axis to beam angle"],
141+
["phi", "degrees", 60, [-360, 360], "orientation", "rotation about beam"],
142+
]
143+
# pylint: enable=bad-whitespace, line-too-long
144+
145+
source = ["lib/polevl.c", "lib/sas_J1.c", "lib/adaptive.c", "cylinderp.c"]
146+
valid = "radius >= 0.0 && length >= 0.0"
147+
have_Fq = True
148+
radius_effective_modes = [
149+
"excluded volume", "equivalent volume sphere", "radius",
150+
"half length", "half min dimension", "half max dimension", "half diagonal",
151+
]
152+
153+
def random():
154+
"""Return a random parameter set for the model."""
155+
volume = 10**np.random.uniform(5, 12)
156+
length = 10**np.random.uniform(-2, 2)*volume**0.333
157+
radius = np.sqrt(volume/length/np.pi)
158+
pars = dict(
159+
#scale=1,
160+
#background=0,
161+
length=length,
162+
radius=radius,
163+
)
164+
return pars
165+
166+
167+
# Test 1-D and 2-D models
168+
qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5)
169+
theta, phi = 80.1534480601659, 10.1510817110481 # (10, 10) in sasview 3.x
170+
tests = [
171+
[{}, 0.2, 0.042761386790780453],
172+
[{}, [0.2], [0.042761386790780453]],
173+
[{"scale": 1., "background": 0.}, [0.01, 0.05, 0.2],
174+
# these numerical results NOT independently verified
175+
[3.01823887e+02, 5.36661653e+01, 4.17613868e-02]],
176+
[{"scale": 1., "background": 0.},
177+
# the longer list here checks F1, F2=I(Q)*V, R_eff, volume, volume_ratio
178+
0.05, 2214.9614083, 26975556.88749548, 73.34013315261608,
179+
502654.8245743669, 1.0],
180+
#2345678901234567890123456789012345678901234567890123456789012345678901234567890
181+
[{"@S": "hardsphere", # MONODISPERSE
182+
"scale": 5., "background": 0., "volfraction": 0.2,
183+
"structure_factor_mode": 0, # normal decoupling approx
184+
"radius_effective_mode": 1, # Reff "excluded volume"
185+
}, [0.01, 0.05, 0.2], [7.35571916e+01, 5.78147797e+01, 4.15623248e-2]
186+
],
187+
[{"@S": "hardsphere",
188+
"scale": 5., "background": 0., "volfraction": 0.2,
189+
"structure_factor_mode": 1, # beta approx
190+
"radius_effective_mode": 1, # Reff "excluded volume"
191+
}, [0.01, 0.05, 0.2], [8.29729770e+01, 5.44206752e+01, 4.17598382e-2]
192+
],
193+
[{"@S": "hardsphere", # POLYDISPERSE
194+
"scale": 5., "background": 0., "volfraction": 0.2,
195+
"radius_pd": 0.2, "radius_pd_n": 15, "length_pd": 0.1,
196+
"structure_factor_mode": 0, # normal decoupling approx
197+
"radius_effective_mode": 1, # Reff "excluded volume"
198+
}, [0.01, 0.05, 0.2], [87.50350037, 63.95202427, 0.27889988]
199+
],
200+
[{"@S": "hardsphere",
201+
"scale": 5., "background": 0., "volfraction": 0.2,
202+
"radius_pd": 0.2, "radius_pd_n": 15, "length_pd": 0.1,
203+
"structure_factor_mode": 1, # beta approx
204+
"radius_effective_mode": 1, # Reff "excluded volume"
205+
}, [0.01, 0.05, 0.2], [132.2101165, 59.89948174, 0.28048784]
206+
],
207+
#
208+
[{'theta': theta, 'phi': phi}, (qx, qy), 0.03514647218513852],
209+
[{'theta': theta, 'phi': phi}, [(qx, qy)], [0.03514647218513852]],
210+
]
211+
del qx, qy, theta, phi # not necessary to delete, but cleaner
212+
213+
def _extend_with_reff_tests(radius, length):
214+
"""Test R_eff and form volume calculations"""
215+
# V and Vr are the same for each R_eff mode
216+
V = pi*radius**2*length # shell volume = form volume for solid objects
217+
Vr = 1.0 # form:shell volume ratio
218+
# Use test value for I(0.2) from above to check Fsq value. Need to
219+
# remove scale and background before testing.
220+
q = 0.2
221+
scale, background = V, 0.001
222+
Fsq = (0.042761386790780453 - background)*scale
223+
F = None # Need target value for <F>
224+
# Various values for R_eff, depending on mode
225+
r_effs = [
226+
0.,
227+
0.5*(0.75*radius*(2.0*radius*length
228+
+ (radius + length)*(pi*radius + length)))**(1./3.),
229+
(0.75*radius**2*length)**(1./3.),
230+
radius,
231+
length/2.,
232+
min(radius, length/2.),
233+
max(radius, length/2.),
234+
np.sqrt(4*radius**2 + length**2)/2.,
235+
]
236+
tests.extend([
237+
({'radius_effective_mode': 0}, q, F, Fsq, r_effs[0], V, Vr),
238+
({'radius_effective_mode': 1}, q, F, Fsq, r_effs[1], V, Vr),
239+
({'radius_effective_mode': 2}, q, F, Fsq, r_effs[2], V, Vr),
240+
({'radius_effective_mode': 3}, q, F, Fsq, r_effs[3], V, Vr),
241+
({'radius_effective_mode': 4}, q, F, Fsq, r_effs[4], V, Vr),
242+
({'radius_effective_mode': 5}, q, F, Fsq, r_effs[5], V, Vr),
243+
({'radius_effective_mode': 6}, q, F, Fsq, r_effs[6], V, Vr),
244+
({'radius_effective_mode': 7}, q, F, Fsq, r_effs[7], V, Vr),
245+
])
246+
247+
# Test Reff and volume with default model parameters
248+
_extend_with_reff_tests(parameters[2][2], parameters[3][2])
249+
del _extend_with_reff_tests
250+
251+
# ADDED by: RKH ON: 18Mar2016 renamed sld's etc

0 commit comments

Comments
 (0)