# Changeset 35d2300 in sasmodels

Ignore:
Timestamp:
Sep 13, 2018 8:08:42 PM (12 months ago)
Branches:
master
Children:
b50e28b
Parents:
a5cb9bc
Message:

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

Location:
example/weights
Files:
3 edited

Unmodified
Removed
• ## example/weights/cyclic_gaussian.py

 r55e82f0 This is eqivalent to a Maier-Saupe distribution with order parameter $P_2 = 1/(2 \sigma^2)$, with $\sigma$ in radians. parameter $a = 1/(2 \sigma^2)$, with $\sigma$ in radians. """ type = "cyclic_gaussian"
• ## example/weights/maier_saupe.py

 r55e82f0 .. 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$ 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" # 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 # 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))
• ## example/weights/maier_saupe_eq.py

 r55e82f0 .. 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$ # 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 P2 = sigma sigma = 1./sqrt(2.*P2) # use the width parameter as the value for Maier-Saupe "a" a = sigma sigma = 1./sqrt(2.*a) # Create a lookup table for finding n points equally spaced # we can scale by an arbitrary scale factor c = exp(m) to get: #     w = exp(m*cos(x)**2)/c = exp(-m*sin(x)**2) yp = np.cumsum(exp(-P2*sin(xp)**2)) yp = np.cumsum(exp(-a*sin(xp)**2)) yp /= yp[-1]
Note: See TracChangeset for help on using the changeset viewer.