1 | """ |
---|
2 | Asymmetric shape integration |
---|
3 | """ |
---|
4 | |
---|
5 | from __future__ import print_function, division |
---|
6 | |
---|
7 | import os, sys |
---|
8 | sys.path.insert(0, os.path.dirname(os.path.dirname(__file__))) |
---|
9 | |
---|
10 | import numpy as np |
---|
11 | import mpmath as mp |
---|
12 | from numpy import pi, sin, cos, sqrt, exp, expm1, degrees, log10 |
---|
13 | from numpy.polynomial.legendre import leggauss |
---|
14 | from scipy.integrate import dblquad, simps, romb, romberg |
---|
15 | import pylab |
---|
16 | |
---|
17 | from sasmodels.special import square |
---|
18 | from sasmodels.special import Gauss20Wt, Gauss20Z |
---|
19 | from sasmodels.special import Gauss76Wt, Gauss76Z |
---|
20 | from sasmodels.special import Gauss150Wt, Gauss150Z |
---|
21 | from sasmodels.special import sas_2J1x_x, sas_sinx_x, sas_3j1x_x |
---|
22 | |
---|
23 | def mp_3j1x_x(x): |
---|
24 | return 3*(mp.sin(x)/x - mp.cos(x))/(x*x) |
---|
25 | def mp_2J1x_x(x): |
---|
26 | return 2*mp.j1(x)/x |
---|
27 | def mp_sinx_x(x): |
---|
28 | return mp.sin(x)/x |
---|
29 | |
---|
30 | SLD = 3 |
---|
31 | SLD_SOLVENT = 6 |
---|
32 | CONTRAST = SLD - SLD_SOLVENT |
---|
33 | |
---|
34 | def make_parallelepiped(a, b, c): |
---|
35 | def Fq(qa, qb, qc): |
---|
36 | siA = sas_sinx_x(0.5*a*qa) |
---|
37 | siB = sas_sinx_x(0.5*b*qb) |
---|
38 | siC = sas_sinx_x(0.5*c*qc) |
---|
39 | return siA * siB * siC |
---|
40 | volume = a*b*c |
---|
41 | norm = volume*CONTRAST**2/10**4 |
---|
42 | return norm, Fq |
---|
43 | |
---|
44 | def make_parallelepiped_mp(a, b, c): |
---|
45 | a, b, c = mp.mpf(a), mp.mpf(b), mp.mpf(c) |
---|
46 | def Fq(qa, qb, qc): |
---|
47 | siA = mp_sinx_x(a*qa/2) |
---|
48 | siB = mp_sinx_x(b*qb/2) |
---|
49 | siC = mp_sinx_x(c*qc/2) |
---|
50 | return siA * siB * siC |
---|
51 | volume = a*b*c |
---|
52 | norm = (volume*CONTRAST**2)/10000 # mpf since volume=a*b*c is mpf |
---|
53 | return norm, Fq |
---|
54 | |
---|
55 | def make_triellip(a, b, c): |
---|
56 | def Fq(qa, qb, qc): |
---|
57 | qr = sqrt((a*qa)**2 + (b*qb)**2 + (c*qc)**2) |
---|
58 | return sas_3j1x_x(qr) |
---|
59 | volume = 4*pi*a*b*c/3 |
---|
60 | norm = volume*CONTRAST**2/10**4 |
---|
61 | return norm, Fq |
---|
62 | |
---|
63 | def make_triellip_mp(a, b, c): |
---|
64 | a, b, c = mp.mpf(a), mp.mpf(b), mp.mpf(c) |
---|
65 | def Fq(qa, qb, qc): |
---|
66 | qr = mp.sqrt((a*qa)**2 + (b*qb)**2 + (c*qc)**2) |
---|
67 | return mp_3j1x_x(qr) |
---|
68 | volume = (4*mp.pi*a*b*c)/3 |
---|
69 | norm = (volume*CONTRAST**2)/10000 # mpf since mp.pi is mpf |
---|
70 | return norm, Fq |
---|
71 | |
---|
72 | def make_cylinder(radius, length): |
---|
73 | def Fq(qa, qb, qc): |
---|
74 | qab = sqrt(qa**2 + qb**2) |
---|
75 | return sas_2J1x_x(qab*radius) * sas_sinx_x(qc*0.5*length) |
---|
76 | volume = pi*radius**2*length |
---|
77 | norm = volume*CONTRAST**2/10**4 |
---|
78 | return norm, Fq |
---|
79 | |
---|
80 | def make_cylinder_mp(radius, length): |
---|
81 | radius, length = mp.mpf(radius), mp.mpf(length) |
---|
82 | def Fq(qa, qb, qc): |
---|
83 | qab = mp.sqrt(qa**2 + qb**2) |
---|
84 | return mp_2J1x_x(qab*radius) * mp_sinx_x((qc*length)/2) |
---|
85 | volume = mp.pi*radius**2*length |
---|
86 | norm = (volume*CONTRAST**2)/10000 # mpf since mp.pi is mpf |
---|
87 | return norm, Fq |
---|
88 | |
---|
89 | def make_sphere(radius): |
---|
90 | def Fq(qa, qb, qc): |
---|
91 | q = sqrt(qa**2 + qb**2 + qc**2) |
---|
92 | return sas_3j1x_x(q*radius) |
---|
93 | volume = 4*pi*radius**3/3 |
---|
94 | norm = volume*CONTRAST**2/10**4 |
---|
95 | return norm, Fq |
---|
96 | |
---|
97 | def make_sphere_mp(radius): |
---|
98 | radius = mp.mpf(radius) |
---|
99 | def Fq(qa, qb, qc): |
---|
100 | q = mp.sqrt(qa**2 + qb**2 + qc**2) |
---|
101 | return mp_3j1x_x(q*radius) |
---|
102 | volume = (4*mp.pi*radius**3)/3 |
---|
103 | norm = (volume*CONTRAST**2)/10000 # mpf since mp.pi is mpf |
---|
104 | return norm, Fq |
---|
105 | |
---|
106 | shape = 'parallelepiped' |
---|
107 | #shape = 'triellip' |
---|
108 | #shape = 'sphere' |
---|
109 | #shape = 'cylinder' |
---|
110 | if shape == 'cylinder': |
---|
111 | #RADIUS, LENGTH = 10, 100000 |
---|
112 | RADIUS, LENGTH = 10, 300 # integer for the sake of mpf |
---|
113 | NORM, KERNEL = make_cylinder(radius=RADIUS, length=LENGTH) |
---|
114 | NORM_MP, KERNEL_MP = make_cylinder_mp(radius=RADIUS, length=LENGTH) |
---|
115 | elif shape == 'triellip': |
---|
116 | #A, B, C = 4450, 14000, 47 |
---|
117 | A, B, C = 445, 140, 47 # integer for the sake of mpf |
---|
118 | NORM, KERNEL = make_triellip(A, B, C) |
---|
119 | NORM_MP, KERNEL_MP = make_triellip_mp(A, B, C) |
---|
120 | elif shape == 'parallelepiped': |
---|
121 | #A, B, C = 4450, 14000, 47 |
---|
122 | A, B, C = 445, 140, 47 # integer for the sake of mpf |
---|
123 | NORM, KERNEL = make_parallelepiped(A, B, C) |
---|
124 | NORM_MP, KERNEL_MP = make_parallelepiped_mp(A, B, C) |
---|
125 | elif shape == 'sphere': |
---|
126 | RADIUS = 50 # integer for the sake of mpf |
---|
127 | NORM, KERNEL = make_sphere(radius=RADIUS) |
---|
128 | NORM_MP, KERNEL_MP = make_sphere_mp(radius=RADIUS) |
---|
129 | else: |
---|
130 | raise ValueError("Unknown shape %r"%shape) |
---|
131 | |
---|
132 | THETA_LOW, THETA_HIGH = 0, pi |
---|
133 | PHI_LOW, PHI_HIGH = 0, 2*pi |
---|
134 | SCALE = 1 |
---|
135 | |
---|
136 | # mathematica code for triaxial_ellipsoid (untested) |
---|
137 | _ = """ |
---|
138 | R[theta_, phi_, a_, b_, c_] := Sqrt[(a Sin[theta]Cos[phi])^2 + (b Sin[theta]Sin[phi])^2 + (c Cos[theta])^2] |
---|
139 | Sphere[q_, r_] := 3 SphericalBesselJ[q r]/(q r) |
---|
140 | V[a_, b_, c_] := 4/3 pi a b c |
---|
141 | Norm[sld_, solvent_, a_, b_, c_] := V[a, b, c] (solvent - sld)^2 |
---|
142 | F[q_, theta_, phi_, a_, b_, c_] := Sphere[q, R[theta, phi, a, b, c]] |
---|
143 | I[q_, sld_, solvent_, a_, b_, c_] := Norm[sld, solvent, a, b, c]/(4 pi) Integrate[F[q, theta, phi, a, b, c]^2 Sin[theta], {phi, 0, 2 pi}, {theta, 0, pi}] |
---|
144 | I[6/10^3, 63/10, 3, 445, 140, 47] |
---|
145 | """ |
---|
146 | |
---|
147 | |
---|
148 | def mp_quad(q, shape): |
---|
149 | evals = [0] |
---|
150 | def integrand(theta, phi): |
---|
151 | evals[0] += 1 |
---|
152 | qab = q*mp.sin(theta) |
---|
153 | qa = qab*mp.cos(phi) |
---|
154 | qb = qab*mp.sin(phi) |
---|
155 | qc = q*mp.cos(theta) |
---|
156 | Zq = KERNEL_MP(qa, qb, qc)**2 |
---|
157 | return Zq*mp.sin(theta) |
---|
158 | ans = mp.quad(integrand, (0, mp.pi), (0, 2*mp.pi)) |
---|
159 | Iq = NORM_MP*ans/(4*mp.pi) |
---|
160 | return evals[0], Iq |
---|
161 | |
---|
162 | def kernel(q, theta, phi): |
---|
163 | """ |
---|
164 | S(q) kernel for paracrystal forms. |
---|
165 | """ |
---|
166 | qab = q*sin(theta) |
---|
167 | qa = qab*cos(phi) |
---|
168 | qb = qab*sin(phi) |
---|
169 | qc = q*cos(theta) |
---|
170 | return NORM*KERNEL(qa, qb, qc)**2 |
---|
171 | |
---|
172 | def scipy_dblquad(q): |
---|
173 | """ |
---|
174 | Compute the integral using scipy dblquad. This gets the correct answer |
---|
175 | eventually, but it is slow. |
---|
176 | """ |
---|
177 | evals = [0] |
---|
178 | def integrand(phi, theta): |
---|
179 | evals[0] += 1 |
---|
180 | Zq = kernel(q, theta=theta, phi=phi) |
---|
181 | return Zq*sin(theta) |
---|
182 | ans = dblquad(integrand, THETA_LOW, THETA_HIGH, lambda x: PHI_LOW, lambda x: PHI_HIGH)[0] |
---|
183 | return evals[0], ans*SCALE/(4*pi) |
---|
184 | |
---|
185 | |
---|
186 | def scipy_romberg_2d(q): |
---|
187 | """ |
---|
188 | Compute the integral using romberg integration. This function does not |
---|
189 | complete in a reasonable time. No idea if it is accurate. |
---|
190 | """ |
---|
191 | evals = [0] |
---|
192 | def inner(phi, theta): |
---|
193 | evals[0] += 1 |
---|
194 | return kernel(q, theta=theta, phi=phi) |
---|
195 | def outer(theta): |
---|
196 | Zq = romberg(inner, PHI_LOW, PHI_HIGH, divmax=100, args=(theta,)) |
---|
197 | return Zq*sin(theta) |
---|
198 | ans = romberg(outer, THETA_LOW, THETA_HIGH, divmax=100) |
---|
199 | return evals[0], ans*SCALE/(4*pi) |
---|
200 | |
---|
201 | |
---|
202 | def semi_romberg(q, n=100): |
---|
203 | """ |
---|
204 | Use 1D romberg integration in phi and regular simpsons rule in theta. |
---|
205 | """ |
---|
206 | evals = [0] |
---|
207 | def inner(phi, theta): |
---|
208 | evals[0] += 1 |
---|
209 | return kernel(q, theta=theta, phi=phi) |
---|
210 | theta = np.linspace(THETA_LOW, THETA_HIGH, n) |
---|
211 | Zq = [romberg(inner, PHI_LOW, PHI_HIGH, divmax=100, args=(t,)) for t in theta] |
---|
212 | ans = simps(np.array(Zq)*sin(theta), dx=theta[1]-theta[0]) |
---|
213 | return evals[0], ans*SCALE/(4*pi) |
---|
214 | |
---|
215 | def gauss_quad(q, n=150): |
---|
216 | """ |
---|
217 | Compute the integral using gaussian quadrature for n = 20, 76 or 150. |
---|
218 | """ |
---|
219 | if n == 20: |
---|
220 | z, w = Gauss20Z, Gauss20Wt |
---|
221 | elif n == 76: |
---|
222 | z, w = Gauss76Z, Gauss76Wt |
---|
223 | elif n == 150: |
---|
224 | z, w = Gauss150Z, Gauss150Wt |
---|
225 | else: |
---|
226 | z, w = leggauss(n) |
---|
227 | theta = (THETA_HIGH-THETA_LOW)*(z + 1)/2 + THETA_LOW |
---|
228 | phi = (PHI_HIGH-PHI_LOW)*(z + 1)/2 + PHI_LOW |
---|
229 | Atheta, Aphi = np.meshgrid(theta, phi) |
---|
230 | Aw = w[None, :] * w[:, None] |
---|
231 | sin_theta = abs(sin(Atheta)) |
---|
232 | Zq = kernel(q=q, theta=Atheta, phi=Aphi) |
---|
233 | # change from [-1,1] x [-1,1] range to [0, pi] x [0, 2 pi] range |
---|
234 | dxdy_stretch = (THETA_HIGH-THETA_LOW)/2 * (PHI_HIGH-PHI_LOW)/2 |
---|
235 | Iq = np.sum(Zq*Aw*sin_theta)*SCALE/(4*pi) * dxdy_stretch |
---|
236 | return n**2, Iq |
---|
237 | |
---|
238 | |
---|
239 | def gridded_integrals(q, n=300): |
---|
240 | """ |
---|
241 | Compute the integral on a regular grid using rectangular, trapezoidal, |
---|
242 | simpsons, and romberg integration. Romberg integration requires that |
---|
243 | the grid be of size n = 2**k + 1. |
---|
244 | """ |
---|
245 | theta = np.linspace(THETA_LOW, THETA_HIGH, n) |
---|
246 | phi = np.linspace(PHI_LOW, PHI_HIGH, n) |
---|
247 | Atheta, Aphi = np.meshgrid(theta, phi) |
---|
248 | Zq = kernel(q=q, theta=Atheta, phi=Aphi) |
---|
249 | Zq *= abs(sin(Atheta)) |
---|
250 | dx, dy = theta[1]-theta[0], phi[1]-phi[0] |
---|
251 | print("rect-%d"%n, n**2, np.sum(Zq)*dx*dy*SCALE/(4*pi)) |
---|
252 | print("trapz-%d"%n, n**2, np.trapz(np.trapz(Zq, dx=dx), dx=dy)*SCALE/(4*pi)) |
---|
253 | print("simpson-%d"%n, n**2, simps(simps(Zq, dx=dx), dx=dy)*SCALE/(4*pi)) |
---|
254 | print("romb-%d"%n, n**2, romb(romb(Zq, dx=dx), dx=dy)*SCALE/(4*pi)) |
---|
255 | |
---|
256 | def plot(q, n=300): |
---|
257 | """ |
---|
258 | Plot the 2D surface that needs to be integrated in order to compute |
---|
259 | the BCC S(q) at a particular q, dnn and d_factor. *n* is the number |
---|
260 | of points in the grid. |
---|
261 | """ |
---|
262 | theta = np.linspace(THETA_LOW, THETA_HIGH, n) |
---|
263 | phi = np.linspace(PHI_LOW, PHI_HIGH, n) |
---|
264 | Atheta, Aphi = np.meshgrid(theta, phi) |
---|
265 | Zq = kernel(q=q, theta=Atheta, phi=Aphi) |
---|
266 | #Zq *= abs(sin(Atheta)) |
---|
267 | pylab.pcolor(degrees(theta), degrees(phi), log10(np.fmax(Zq, 1.e-6))) |
---|
268 | pylab.axis('tight') |
---|
269 | pylab.title("%s Z(q) for q=%g" % (KERNEL.__name__, q)) |
---|
270 | pylab.xlabel("theta (degrees)") |
---|
271 | pylab.ylabel("phi (degrees)") |
---|
272 | cbar = pylab.colorbar() |
---|
273 | cbar.set_label('log10 S(q)') |
---|
274 | pylab.show() |
---|
275 | |
---|
276 | if __name__ == "__main__": |
---|
277 | Qstr = '0.005' |
---|
278 | #Qstr = '0.8' |
---|
279 | #Qstr = '0.0003' |
---|
280 | Q = float(Qstr) |
---|
281 | if shape == 'sphere': |
---|
282 | print("exact", NORM*sas_3j1x_x(Q*RADIUS)**2) |
---|
283 | print("gauss-20", *gauss_quad(Q, n=20)) |
---|
284 | print("gauss-76", *gauss_quad(Q, n=76)) |
---|
285 | print("gauss-150", *gauss_quad(Q, n=150)) |
---|
286 | print("gauss-500", *gauss_quad(Q, n=500)) |
---|
287 | print("dblquad", *scipy_dblquad(Q)) |
---|
288 | print("semi-romberg-100", *semi_romberg(Q, n=100)) |
---|
289 | print("romberg", *scipy_romberg_2d(Q)) |
---|
290 | #gridded_integrals(Q, n=2**8+1) |
---|
291 | gridded_integrals(Q, n=2**10+1) |
---|
292 | #gridded_integrals(Q, n=2**13+1) |
---|
293 | #gridded_integrals(Q, n=2**15+1) |
---|
294 | with mp.workprec(100): |
---|
295 | print("mpmath", *mp_quad(mp.mpf(Qstr), shape)) |
---|
296 | #plot(Q, n=200) |
---|