source: sasmodels/example/weights/maier_saupe_eq.py

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

Relate order parameter P_2 to 'a' value in Maier-Saupe dist

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