source: sasmodels/example/weights/maier_saupe.py @ a5516b1

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

print sigma, P_2 and P_2 inverse for a given a parameter

  • Property mode set to 100644
File size: 4.5 KB
Line 
1from __future__ import print_function
2import numpy as np
3from numpy import exp, sin, degrees, radians, pi, sqrt
4
5from sasmodels.weights import Dispersion as BaseDispersion
6
7class Dispersion(BaseDispersion):
8    r"""
9    Maier-Saupe dispersion on orientation.
10
11    .. math:
12
13        w(\theta) = e^{a{\cos^2 \theta}}
14
15    This provides a close match to the gaussian distribution for
16    low angles, but the tails are limited to $\pm 90^\circ$.  For $a \ll 1$
17    the distribution is approximately uniform.  The usual polar coordinate
18    projection applies, with $\theta$ weights scaled by $\cos \theta$
19    and $\phi$ weights unscaled.
20
21    This is equivalent to a cyclic gaussian distribution
22    $w(\theta) = e^{-sin^2(\theta)/(2\a^2)}.
23
24    Note that the incorrect distribution is used for $a=0$.  With the
25    distribution parameter labelled as *width*, the value *width=0* is
26    is assumed to be completely oriented and no polydispersity distribution
27    is generated.  However a value of $a=0$ should be completely unoriented.
28    Any small value (e.g., 0.01 or lower) should suffice.
29
30    The order parameter $P_2$ is defined as
31
32    .. math:
33
34        P(a, \beta) = \frac{e^{a \cos^2 \beta}}{
35                       4\pi \int_0^{\pi/2} e^{a \cos^2 \beta} \sin \beta\,d\beta}
36
37        P_2 = 4\pi \int_0^{\pi/2} \frac{3 \cos^2 \beta - 1)}{2}
38                        P(a, \beta) \sin \beta\,d\beta
39
40    where $a$ is the distribution width $\sigma$ handed to the weights function.
41    There is a somewhat complicated closed form solution
42
43    .. math:
44
45        P_2 = \frac{3e^a}{2\sqrt{\pi a} E_a} - \frac{3}{4a} - \frac{1}{2}
46
47    where $E_a = \mathrm{erfi}(\sqrt{a})$ is the imaginary error function,
48    which can be coded in python as::
49
50        from numpy import pi, sqrt, exp
51        from scipy.special import erfi
52
53        def P_2(a):
54            # protect against overflow with a Taylor series at infinity
55            if a <= 700:
56                r = exp(a)/sqrt(pi*a)/erfi(sqrt(a))
57            else:
58                r = 1/((((6.5525/a + 1.875)/a + 0.75)/a + 0.5)/a + 1)
59            return 1.5*r - 0.75/a - 0.5
60
61    Given an order parameter $S = P_2(a)$, one can also solve for the
62    equivalent $a$:
63
64        from scipy.optimize import fsolve
65
66        def P_2_inv(S):
67            return fsolve(lambda x: P_2(x) - S, 1.0)[0]
68
69    References
70    ----------
71
72    [1] Hardouin, et al., 1995. SANS study of a semiflexible main chain
73    liquid crystalline polyether. *Macromolecules* 28, 5427-5433.
74    """
75    type = "maier_saupe"
76    default = dict(npts=35, width=1, nsigmas=None)
77
78    # Note: center is always zero for orientation distributions
79    def _weights(self, center, sigma, lb, ub):
80        # use the width parameter as the value for Maier-Saupe "a"
81        # and find the equivalent width sigma
82        a = sigma
83        sigma = 1./sqrt(2.*a)
84
85        # Limit width to +/- 90 degrees
86        width = min(self.nsigmas*sigma, pi/2)
87        x = np.linspace(-width, width, self.npts)
88
89        # Truncate the distribution in case the parameter value is limited
90        x[(x >= radians(lb)) & (x <= radians(ub))]
91
92        # Return orientation in degrees with Maier-Saupe weights
93        # Note: weights are normalized to sum to 1, so we can scale
94        # by an arbitrary scale factor c = exp(m) to get:
95        #     w = exp(m*cos(x)**2)/c = exp(-m*sin(x)**2)
96        return degrees(x), exp(-a*sin(x)**2)
97
98
99
100def P_2(a):
101    from numpy import pi, sqrt, exp
102    from scipy.special import erfi
103
104    # Double precision e^x overflows so do a Taylor expansion around infinity
105    # for large values of a.  With five terms it is accurate to 12 digits
106    # at a = 700, and 7 digits at a = 75.
107    if a <= 700:
108        r = exp(a)/sqrt(pi*a)/erfi(sqrt(a))
109    else:
110        r = 1/((((6.5525/a + 1.875)/a + 0.75)/a + 0.5)/a + 1)
111    return 1.5*r - 0.75/a - 0.5
112
113def P_2_inv(S):
114    from scipy.optimize import fsolve
115    return fsolve(lambda x: P_2(x) - S, 1.0)[0]
116
117def P_2_numerical(a):
118    from scipy.integrate import romberg
119    from numpy import cos, sin, pi, exp
120    def num(beta):
121        return (1.5 * cos(beta)**2 - 0.5) * exp(a * cos(beta)**2) * sin(beta)
122    def denom(beta):
123        return exp(a * cos(beta)**2) * sin(beta)
124    return romberg(num, 0, pi/2) / romberg(denom, 0, pi/2)
125
126if __name__ == "__main__":
127    import sys
128    a = float(sys.argv[1])
129    sigma = 1/(2*radians(a)**2)
130    #print("P_2", P_2(a), "difference from integral", P_2(a) - P_2_numerical(a))
131    print("a=%g, sigma=%g, P_2=%g, P_2_inv(P_2(a))-a=%g"
132          % (a, sigma, P_2(a), P_2_inv(P_2(a))-a))
Note: See TracBrowser for help on using the repository browser.