1 | from __future__ import print_function |
---|
2 | import numpy as np |
---|
3 | from numpy import exp, sin, degrees, radians, pi, sqrt |
---|
4 | |
---|
5 | from sasmodels.weights import Dispersion as BaseDispersion |
---|
6 | |
---|
7 | class 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 | The order parameter $P_2$ is defined as |
---|
25 | |
---|
26 | .. math: |
---|
27 | |
---|
28 | P(a, \beta) = \frac{e^{a \cos^2 \beta}}{ |
---|
29 | 4\pi \int_0^{\pi/2} e^{a \cos^2 \beta} \sin \beta\,d\beta} |
---|
30 | |
---|
31 | P_2 = 4\pi \int_0^{\pi/2} \frac{3 \cos^2 \beta - 1)}{2} |
---|
32 | P(a, \beta) \sin \beta\,d\beta |
---|
33 | |
---|
34 | where $a$ is the distribution width $\sigma$ handed to the weights function. |
---|
35 | There is a somewhat complicated closed form solution |
---|
36 | |
---|
37 | .. math: |
---|
38 | |
---|
39 | P_2 = \frac{3e^a}{2\sqrt{\pi a} E_a} - \frac{3}{4a} - \frac{1}{2} |
---|
40 | |
---|
41 | where $E_a = \mathrm{erfi}(\sqrt{a})$ is the imaginary error function, |
---|
42 | which can be coded in python as:: |
---|
43 | |
---|
44 | from numpy import pi, sqrt, exp |
---|
45 | from scipy.special import erfi |
---|
46 | |
---|
47 | def P_2(a): |
---|
48 | # protect against overflow with a Taylor series at infinity |
---|
49 | if a <= 700: |
---|
50 | r = exp(a)/sqrt(pi*a)/erfi(sqrt(a)) |
---|
51 | else: |
---|
52 | r = 1/((((6.5525/a + 1.875)/a + 0.75)/a + 0.5)/a + 1) |
---|
53 | return 1.5*r - 0.75/a - 0.5 |
---|
54 | |
---|
55 | References |
---|
56 | ---------- |
---|
57 | |
---|
58 | [1] Hardouin, et al., 1995. SANS study of a semiflexible main chain |
---|
59 | liquid crystalline polyether. *Macromolecules* 28, 5427â5433. |
---|
60 | """ |
---|
61 | type = "maier_saupe" |
---|
62 | default = dict(npts=35, width=1, nsigmas=None) |
---|
63 | |
---|
64 | # Note: center is always zero for orientation distributions |
---|
65 | def _weights(self, center, sigma, lb, ub): |
---|
66 | # use the width parameter as the value for Maier-Saupe "a" |
---|
67 | # and find the equivalent width sigma |
---|
68 | a = sigma |
---|
69 | sigma = 1./sqrt(2.*a) |
---|
70 | |
---|
71 | # Limit width to +/- 90 degrees |
---|
72 | width = min(self.nsigmas*sigma, pi/2) |
---|
73 | x = np.linspace(-width, width, self.npts) |
---|
74 | |
---|
75 | # Truncate the distribution in case the parameter value is limited |
---|
76 | x[(x >= radians(lb)) & (x <= radians(ub))] |
---|
77 | |
---|
78 | # Return orientation in degrees with Maier-Saupe weights |
---|
79 | # Note: weights are normalized to sum to 1, so we can scale |
---|
80 | # by an arbitrary scale factor c = exp(m) to get: |
---|
81 | # w = exp(m*cos(x)**2)/c = exp(-m*sin(x)**2) |
---|
82 | return degrees(x), exp(-a*sin(x)**2) |
---|
83 | |
---|
84 | |
---|
85 | |
---|
86 | def P_2(a): |
---|
87 | from numpy import pi, sqrt, exp |
---|
88 | from scipy.special import erfi |
---|
89 | |
---|
90 | # Double precision e^x overflows so do a Taylor expansion around infinity |
---|
91 | # for large values of a. With five terms it is accurate to 12 digits |
---|
92 | # at a = 700, and 7 digits at a = 75. |
---|
93 | if a <= 700: |
---|
94 | r = exp(a)/sqrt(pi*a)/erfi(sqrt(a)) |
---|
95 | else: |
---|
96 | r = 1/((((6.5525/a + 1.875)/a + 0.75)/a + 0.5)/a + 1) |
---|
97 | return 1.5*r - 0.75/a - 0.5 |
---|
98 | |
---|
99 | def P_2_numerical(a): |
---|
100 | from scipy.integrate import romberg |
---|
101 | from numpy import cos, sin, pi, exp |
---|
102 | def num(beta): |
---|
103 | return (1.5 * cos(beta)**2 - 0.5) * exp(a * cos(beta)**2) * sin(beta) |
---|
104 | def denom(beta): |
---|
105 | return exp(a * cos(beta)**2) * sin(beta) |
---|
106 | return romberg(num, 0, pi/2) / romberg(denom, 0, pi/2) |
---|
107 | |
---|
108 | if __name__ == "__main__": |
---|
109 | import sys |
---|
110 | a = float(sys.argv[1]) |
---|
111 | #print("P_2", P_2(a), "difference from integral", P_2(a) - P_2_numerical(a)) |
---|
112 | print("P_2", P_2(a)) |
---|