[55e82f0] | 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 | |
---|
[35d2300] | 13 | w(\theta) = e^{a{\cos^2 \theta}} |
---|
[55e82f0] | 14 | |
---|
| 15 | This provides a close match to the gaussian distribution for |
---|
[35d2300] | 16 | low angles, but the tails are limited to $\pm 90^\circ$. For $a \ll 1$ |
---|
[55e82f0] | 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 |
---|
[35d2300] | 22 | $w(\theta) = e^{-sin^2(\theta)/(2\a^2)}. |
---|
| 23 | |
---|
[b50e28b] | 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 | |
---|
[35d2300] | 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 | |
---|
[cd23f31] | 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 | |
---|
[35d2300] | 69 | References |
---|
| 70 | ---------- |
---|
| 71 | |
---|
| 72 | [1] Hardouin, et al., 1995. SANS study of a semiflexible main chain |
---|
[b50e28b] | 73 | liquid crystalline polyether. *Macromolecules* 28, 5427-5433. |
---|
[55e82f0] | 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): |
---|
[35d2300] | 80 | # use the width parameter as the value for Maier-Saupe "a" |
---|
[55e82f0] | 81 | # and find the equivalent width sigma |
---|
[35d2300] | 82 | a = sigma |
---|
| 83 | sigma = 1./sqrt(2.*a) |
---|
[55e82f0] | 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) |
---|
[35d2300] | 96 | return degrees(x), exp(-a*sin(x)**2) |
---|
| 97 | |
---|
| 98 | |
---|
| 99 | |
---|
| 100 | def 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 | |
---|
[cd23f31] | 113 | def P_2_inv(S): |
---|
| 114 | from scipy.optimize import fsolve |
---|
| 115 | return fsolve(lambda x: P_2(x) - S, 1.0)[0] |
---|
| 116 | |
---|
[35d2300] | 117 | def 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 | |
---|
| 126 | if __name__ == "__main__": |
---|
| 127 | import sys |
---|
| 128 | a = float(sys.argv[1]) |
---|
[a5516b1] | 129 | sigma = 1/(2*radians(a)**2) |
---|
[35d2300] | 130 | #print("P_2", P_2(a), "difference from integral", P_2(a) - P_2_numerical(a)) |
---|
[a5516b1] | 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)) |
---|