1 | from __future__ import print_function, division |
---|
2 | |
---|
3 | import numpy as np |
---|
4 | from numpy import pi, sin, cos, exp, expm1, degrees, log10 |
---|
5 | from scipy.integrate import dblquad, simps, romb, romberg |
---|
6 | import pylab |
---|
7 | |
---|
8 | from sasmodels.special import square |
---|
9 | from sasmodels.special import Gauss20Wt, Gauss20Z |
---|
10 | from sasmodels.special import Gauss76Wt, Gauss76Z |
---|
11 | from sasmodels.special import Gauss150Wt, Gauss150Z |
---|
12 | |
---|
13 | def square(x): |
---|
14 | return x**2 |
---|
15 | |
---|
16 | Q = 0.001 |
---|
17 | DNN = 220. |
---|
18 | D_FACTOR = 0.06 |
---|
19 | RADIUS = 40.0 |
---|
20 | SLD = 3.0 |
---|
21 | SLD_SOLVENT = 6.3 |
---|
22 | |
---|
23 | def kernel(q, dnn, d_factor, theta, phi): |
---|
24 | qab = q*sin(theta) |
---|
25 | qa = qab*cos(phi) |
---|
26 | qb = qab*sin(phi) |
---|
27 | qc = q*cos(theta) |
---|
28 | |
---|
29 | if 0: # sc |
---|
30 | a1, a2, a3 = qa, qb, qc |
---|
31 | dcos = dnn |
---|
32 | if 1: # bcc |
---|
33 | a1 = +qa - qc + qb |
---|
34 | a2 = +qa + qc - qb |
---|
35 | a3 = -qa + qc + qb |
---|
36 | dcos = dnn/2 |
---|
37 | |
---|
38 | arg = 0.5*square(dnn*d_factor)*(a1**2 + a2**2 + a3**2) |
---|
39 | exp_arg = exp(-arg) |
---|
40 | den = [((exp_arg - 2*cos(dcos*a))*exp_arg + 1.0) for a in (a1, a2, a3)] |
---|
41 | Sq = -expm1(-2*arg)**3/np.prod(den, axis=0) |
---|
42 | return Sq |
---|
43 | |
---|
44 | |
---|
45 | def scipy_dblquad(q=Q, dnn=DNN, d_factor=D_FACTOR): |
---|
46 | def integrand(theta, phi): |
---|
47 | Sq = kernel(q=q, dnn=dnn, d_factor=d_factor, theta=theta, phi=phi) |
---|
48 | return Sq*sin(theta) |
---|
49 | return dblquad(integrand, 0, pi/2, lambda x: 0, lambda x: pi/2)[0]*8/(4*pi) |
---|
50 | |
---|
51 | |
---|
52 | def scipy_romberg(q=Q, dnn=DNN, d_factor=D_FACTOR): |
---|
53 | def inner(phi, theta): |
---|
54 | return kernel(q=q, dnn=dnn, d_factor=d_factor, theta=theta, phi=phi) |
---|
55 | def outer(theta): |
---|
56 | return romberg(inner, 0, pi/2, divmax=100, args=(theta,))*sin(theta) |
---|
57 | return romberg(outer, 0, pi/2, divmax=100)*8/(4*pi) |
---|
58 | |
---|
59 | |
---|
60 | def gauss_quad(q=Q, dnn=DNN, d_factor=D_FACTOR, n=150): |
---|
61 | if n == 20: |
---|
62 | z, w = Gauss20Z, Gauss20Wt |
---|
63 | elif n == 76: |
---|
64 | z, w = Gauss76Z, Gauss76Wt |
---|
65 | else: |
---|
66 | z, w = Gauss150Z, Gauss150Wt |
---|
67 | theta = pi/4*(z + 1) |
---|
68 | phi = pi/4*(z + 1) |
---|
69 | Atheta, Aphi = np.meshgrid(theta, phi) |
---|
70 | Aw = w[None, :] * w[:, None] |
---|
71 | sin_theta = np.fmax(abs(sin(Atheta)), 1e-6) |
---|
72 | Sq = kernel(q=q, dnn=dnn, d_factor=d_factor, theta=Atheta, phi=Aphi) |
---|
73 | return np.sum(Sq*Aw*sin_theta)*8/(4*pi) |
---|
74 | |
---|
75 | |
---|
76 | def gridded_integrals(q=0.001, dnn=DNN, d_factor=D_FACTOR, n=300): |
---|
77 | theta = np.linspace(0, pi/2, n) |
---|
78 | phi = np.linspace(0, pi/2, n) |
---|
79 | Atheta, Aphi = np.meshgrid(theta, phi) |
---|
80 | Sq = kernel(q=Q, dnn=dnn, d_factor=d_factor, theta=Atheta, phi=Aphi) |
---|
81 | Sq *= abs(sin(Atheta)) |
---|
82 | dx, dy = theta[1]-theta[0], phi[1]-phi[0] |
---|
83 | print("rect", n, np.sum(Sq)*dx*dy*8/(4*pi)) |
---|
84 | print("trapz", n, np.trapz(np.trapz(Sq, dx=dx), dx=dy)*8/(4*pi)) |
---|
85 | print("simpson", n, simps(simps(Sq, dx=dx), dx=dy)*8/(4*pi)) |
---|
86 | print("romb", n, romb(romb(Sq, dx=dx), dx=dy)*8/(4*pi)) |
---|
87 | |
---|
88 | def plot(q=0.001, dnn=DNN, d_factor=D_FACTOR, n=300): |
---|
89 | #theta = np.linspace(0, pi, n) |
---|
90 | #phi = np.linspace(0, 2*pi, n) |
---|
91 | theta = np.linspace(0, pi/2, n) |
---|
92 | phi = np.linspace(0, pi/2, n) |
---|
93 | Atheta, Aphi = np.meshgrid(theta, phi) |
---|
94 | Sq = kernel(q=Q, dnn=dnn, d_factor=d_factor, theta=Atheta, phi=Aphi) |
---|
95 | Sq *= abs(sin(Atheta)) |
---|
96 | pylab.pcolor(degrees(theta), degrees(phi), log10(np.fmax(Sq, 1.e-6))) |
---|
97 | pylab.title("I(q) for q=%g"%q) |
---|
98 | pylab.xlabel("theta (degrees)") |
---|
99 | pylab.ylabel("phi (degrees)") |
---|
100 | cbar = pylab.colorbar() |
---|
101 | cbar.set_label('log10 I(q)') |
---|
102 | pylab.show() |
---|
103 | |
---|
104 | if __name__ == "__main__": |
---|
105 | print("gauss", 20, gauss_quad(n=20)) |
---|
106 | print("gauss", 76, gauss_quad(n=76)) |
---|
107 | print("gauss", 150, gauss_quad(n=150)) |
---|
108 | print("dblquad", scipy_dblquad()) |
---|
109 | gridded_integrals(n=2**8+1) |
---|
110 | gridded_integrals(n=2**10+1) |
---|
111 | gridded_integrals(n=2**13+1) |
---|
112 | #print("romberg", scipy_romberg()) |
---|
113 | #plot(n=300) |
---|