1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
|
import numpy as np
from numpy import exp, sin, degrees, radians, pi, sqrt
from sasmodels.weights import Dispersion as BaseDispersion
class Dispersion(BaseDispersion):
r"""
Maier-Saupe dispersion on orientation (equal weights).
.. math:
w(\theta) = e^{a{\cos^2 \theta}}
This provides a close match to the gaussian distribution for
low angles, but the tails are limited to $\pm 90^\circ$. For $a \ll 1$
the distribution is approximately uniform. The usual polar coordinate
projection applies, with $\theta$ weights scaled by $\cos \theta$
and $\phi$ weights unscaled.
This is equivalent to a cyclic gaussian distribution
$w(\theta) = e^{-sin^2(\theta)/(2\sigma^2)}.
The $\theta$ points are spaced such that each interval has an
equal contribution to the distribution.
This works surprisingly poorly. Try::
$ sascomp cylinder -2d theta=45 phi=20 phi_pd_type=maier_saupe_eq \
phi_pd_n=100,1000 radius=50 length=2*radius -midq phi_pd=5
Leaving it here for others to improve.
"""
type = "maier_saupe_eq"
default = dict(npts=35, width=1, nsigmas=None)
# Note: center is always zero for orientation distributions
def _weights(self, center, sigma, lb, ub):
# use the width parameter as the value for Maier-Saupe "a"
a = sigma
sigma = 1./sqrt(2.*a)
# Create a lookup table for finding n points equally spaced
# in the cumulative density function.
# Limit width to +/-90 degrees.
width = min(self.nsigmas*sigma, pi/2)
xp = np.linspace(-width, width, max(self.npts*10, 100))
# Compute CDF. Since we normalized the sum of the weights to 1,
# we can scale by an arbitrary scale factor c = exp(m) to get:
# w = exp(m*cos(x)**2)/c = exp(-m*sin(x)**2)
yp = np.cumsum(exp(-a*sin(xp)**2))
yp /= yp[-1]
# Find the mid-points of the equal-weighted intervals in the CDF
y = np.linspace(0, 1, self.npts+2)[1:-1]
x = np.interp(y, yp, xp)
wx = np.ones(self.npts)
# Truncate the distribution in case the parameter value is limited
index = (x >= lb) & (x <= ub)
x, wx = x[index], wx[index]
return degrees(x), wx
|
