Changeset 35d2300 in sasmodels
 Timestamp:
 Sep 13, 2018 8:08:42 PM (4 months ago)
 Branches:
 ticket608userdefinedweights
 Children:
 b50e28b
 Parents:
 a5cb9bc
 Location:
 example/weights
 Files:

 3 edited
Legend:
 Unmodified
 Added
 Removed

example/weights/cyclic_gaussian.py
r55e82f0 r35d2300 19 19 20 20 This is eqivalent to a MaierSaupe distribution with order 21 parameter $ P_2= 1/(2 \sigma^2)$, with $\sigma$ in radians.21 parameter $a = 1/(2 \sigma^2)$, with $\sigma$ in radians. 22 22 """ 23 23 type = "cyclic_gaussian" 
example/weights/maier_saupe.py
r55e82f0 r35d2300 11 11 .. math: 12 12 13 w(\theta) = e^{ P_2{\cos^2 \theta}}13 w(\theta) = e^{a{\cos^2 \theta}} 14 14 15 15 This provides a close match to the gaussian distribution for 16 low angles, but the tails are limited to $\pm 90^\circ$. For $ P_2\ll 1$16 low angles, but the tails are limited to $\pm 90^\circ$. For $a \ll 1$ 17 17 the distribution is approximately uniform. The usual polar coordinate 18 18 projection applies, with $\theta$ weights scaled by $\cos \theta$ … … 20 20 21 21 This is equivalent to a cyclic gaussian distribution 22 $w(\theta) = e^{sin^2(\theta)/(2\P_2^2)}. 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. 23 60 """ 24 61 type = "maier_saupe" … … 27 64 # Note: center is always zero for orientation distributions 28 65 def _weights(self, center, sigma, lb, ub): 29 # use the width parameter as the value for MaierSaupe P_266 # use the width parameter as the value for MaierSaupe "a" 30 67 # and find the equivalent width sigma 31 P2= sigma32 sigma = 1./sqrt(2.* P2)68 a = sigma 69 sigma = 1./sqrt(2.*a) 33 70 34 71 # Limit width to +/ 90 degrees … … 43 80 # by an arbitrary scale factor c = exp(m) to get: 44 81 # w = exp(m*cos(x)**2)/c = exp(m*sin(x)**2) 45 return degrees(x), exp(P2*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)) 
example/weights/maier_saupe_eq.py
r55e82f0 r35d2300 10 10 .. math: 11 11 12 w(\theta) = e^{ P_2{\cos^2 \theta}}12 w(\theta) = e^{a{\cos^2 \theta}} 13 13 14 14 This provides a close match to the gaussian distribution for 15 low angles, but the tails are limited to $\pm 90^\circ$. For $ P_2\ll 1$15 low angles, but the tails are limited to $\pm 90^\circ$. For $a \ll 1$ 16 16 the distribution is approximately uniform. The usual polar coordinate 17 17 projection applies, with $\theta$ weights scaled by $\cos \theta$ … … 36 36 # Note: center is always zero for orientation distributions 37 37 def _weights(self, center, sigma, lb, ub): 38 # use the width parameter as the value for MaierSaupe P_239 P2= sigma40 sigma = 1./sqrt(2.* P2)38 # use the width parameter as the value for MaierSaupe "a" 39 a = sigma 40 sigma = 1./sqrt(2.*a) 41 41 42 42 # Create a lookup table for finding n points equally spaced … … 49 49 # we can scale by an arbitrary scale factor c = exp(m) to get: 50 50 # w = exp(m*cos(x)**2)/c = exp(m*sin(x)**2) 51 yp = np.cumsum(exp( P2*sin(xp)**2))51 yp = np.cumsum(exp(a*sin(xp)**2)) 52 52 yp /= yp[1] 53 53
Note: See TracChangeset
for help on using the changeset viewer.