Changeset 35d2300 in sasmodels for example/weights/maier_saupe.py

Timestamp:

Sep 13, 2018 6:08:42 PM (6 years ago)

Author:

Paul Kienzle <pkienzle@…>

Branches:

master

Children:

b50e28b

Parents:

a5cb9bc

Message:

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

File:

• example/weights/maier_saupe.py (modified) (4 diffs)

example/weights/maier_saupe.py

@@ -11 +11 @@
     .. math:
 
-        w(\theta) = e^{P_2{\cos^2 \theta}}
+        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 $P_2 \ll 1$
+    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$
@@ -20 +20 @@
 
     This is equivalent to a cyclic gaussian distribution
-    $w(\theta) = e^{-sin^2(\theta)/(2\P_2^2)}.
+    $w(\theta) = e^{-sin^2(\theta)/(2\a^2)}.
+
+    The order parameter $P_2$ is defined as
+
+    .. math:
+
+        P(a, \beta) = \frac{e^{a \cos^2 \beta}}{
+                       4\pi \int_0^{\pi/2} e^{a \cos^2 \beta} \sin \beta\,d\beta}
+
+        P_2 = 4\pi \int_0^{\pi/2} \frac{3 \cos^2 \beta - 1)}{2}
+                        P(a, \beta) \sin \beta\,d\beta
+
+    where $a$ is the distribution width $\sigma$ handed to the weights function.
+    There is a somewhat complicated closed form solution
+
+    .. math:
+
+        P_2 = \frac{3e^a}{2\sqrt{\pi a} E_a} - \frac{3}{4a} - \frac{1}{2}
+
+    where $E_a = \mathrm{erfi}(\sqrt{a})$ is the imaginary error function,
+    which can be coded in python as::
+
+        from numpy import pi, sqrt, exp
+        from scipy.special import erfi
+
+        def P_2(a):
+            # protect against overflow with a Taylor series at infinity
+            if a <= 700:
+                r = exp(a)/sqrt(pi*a)/erfi(sqrt(a))
+            else:
+                r = 1/((((6.5525/a + 1.875)/a + 0.75)/a + 0.5)/a + 1)
+            return 1.5*r - 0.75/a - 0.5
+
+    References
+    ----------
+
+    [1] Hardouin, et al., 1995. SANS study of a semiflexible main chain
+    liquid crystalline polyether. *Macromolecules* 28, 5427–5433.
     """
     type = "maier_saupe"
@@ -27 +64 @@
     # 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 P_2
+        # use the width parameter as the value for Maier-Saupe "a"
         # and find the equivalent width sigma
-        P2 = sigma
-        sigma = 1./sqrt(2.*P2)
+        a = sigma
+        sigma = 1./sqrt(2.*a)
 
         # Limit width to +/- 90 degrees
@@ -43 +80 @@
         # by an arbitrary scale factor c = exp(m) to get:
         #     w = exp(m*cos(x)**2)/c = exp(-m*sin(x)**2)
-        return degrees(x), exp(-P2*sin(x)**2)
+        return degrees(x), exp(-a*sin(x)**2)
+
+
+
+def P_2(a):
+    from numpy import pi, sqrt, exp
+    from scipy.special import erfi
+
+    # Double precision e^x overflows so do a Taylor expansion around infinity
+    # for large values of a.  With five terms it is accurate to 12 digits
+    # at a = 700, and 7 digits at a = 75.
+    if a <= 700:
+        r = exp(a)/sqrt(pi*a)/erfi(sqrt(a))
+    else:
+        r = 1/((((6.5525/a + 1.875)/a + 0.75)/a + 0.5)/a + 1)
+    return 1.5*r - 0.75/a - 0.5
+
+def P_2_numerical(a):
+    from scipy.integrate import romberg
+    from numpy import cos, sin, pi, exp
+    def num(beta):
+        return (1.5 * cos(beta)**2 - 0.5) * exp(a * cos(beta)**2) * sin(beta)
+    def denom(beta):
+        return exp(a * cos(beta)**2) * sin(beta)
+    return romberg(num, 0, pi/2) / romberg(denom, 0, pi/2)
+
+if __name__ == "__main__":
+    import sys
+    a = float(sys.argv[1])
+    #print("P_2", P_2(a), "difference from integral", P_2(a) - P_2_numerical(a))
+    print("P_2", P_2(a))