source: sasmodels/explore/asymint.py @ 2bccb5a

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 2bccb5a was 4f611f1, checked in by Paul Kienzle <pkienzle@…>, 7 years ago

play with cylinder integral

  • Property mode set to 100644
File size: 5.4 KB
Line 
1"""
2Asymmetric shape integration
3"""
4
5from __future__ import print_function, division
6
7import os, sys
8sys.path.insert(0, os.path.dirname(os.path.dirname(__file__)))
9
10import numpy as np
11from numpy import pi, sin, cos, sqrt, exp, expm1, degrees, log10
12from scipy.integrate import dblquad, simps, romb, romberg
13import pylab
14
15from sasmodels.special import square
16from sasmodels.special import Gauss20Wt, Gauss20Z
17from sasmodels.special import Gauss76Wt, Gauss76Z
18from sasmodels.special import Gauss150Wt, Gauss150Z
19from sasmodels.special import sas_2J1x_x, sas_sinx_x, sas_3j1x_x
20
21SLD = 3.0
22SLD_SOLVENT = 6.3
23CONTRAST = SLD - SLD_SOLVENT
24
25def make_cylinder(radius, length):
26    def cylinder(qa, qb, qc):
27        qab = sqrt(qa**2 + qb**2)
28        return sas_2J1x_x(qab*radius) * sas_sinx_x(qc*0.5*length)
29    volume = pi*radius**2*length
30    norm = 1e-4*volume*CONTRAST**2
31    return norm, cylinder
32
33def make_sphere(radius):
34    def sphere(qa, qb, qc):
35        qab = sqrt(qa**2 + qb**2)
36        q = sqrt(qab**2 + qc**2)
37        return sas_3j1x_x(q*radius)
38    volume = 4*pi*radius**3/3
39    norm = 1e-4*volume*CONTRAST**2
40    return norm, sphere
41   
42#NORM, KERNEL = make_cylinder(radius=10., length=100000.)
43NORM, KERNEL = make_cylinder(radius=10., length=30.)
44#NORM, KERNEL = make_sphere(radius=50.)
45
46THETA_LOW, THETA_HIGH = 0, pi
47PHI_LOW, PHI_HIGH = 0, 2*pi
48SCALE = 1
49
50
51def kernel(q, theta, phi):
52    """
53    S(q) kernel for paracrystal forms.
54    """
55    qab = q*sin(theta)
56    qa = qab*cos(phi)
57    qb = qab*sin(phi)
58    qc = q*cos(theta)
59    return NORM*KERNEL(qa, qb, qc)**2
60
61
62def scipy_dblquad(q):
63    """
64    Compute the integral using scipy dblquad.  This gets the correct answer
65    eventually, but it is slow.
66    """
67    evals = [0]
68    def integrand(theta, phi):
69        evals[0] += 1
70        Zq = kernel(q, theta=theta, phi=phi)
71        return Zq*sin(theta)
72    ans = dblquad(integrand, THETA_LOW, THETA_HIGH, lambda x: PHI_LOW, lambda x: PHI_HIGH)[0]*SCALE/(4*pi)
73    print("dblquad evals =", evals[0])
74    return ans
75
76
77def scipy_romberg_2d(q):
78    """
79    Compute the integral using romberg integration.  This function does not
80    complete in a reasonable time.  No idea if it is accurate.
81    """
82    def inner(phi, theta):
83        return kernel(q, theta=theta, phi=phi)
84    def outer(theta):
85        return romberg(inner, PHI_LOW, PHI_HIGH, divmax=100, args=(theta,))*sin(theta)
86    return romberg(outer, THETA_LOW, THETA_HIGH, divmax=100)*SCALE/(4*pi)
87
88
89def semi_romberg(q, n=100):
90    """
91    Use 1D romberg integration in phi and regular simpsons rule in theta.
92    """
93    evals = [0]
94    def inner(phi, theta):
95        evals[0] += 1
96        return kernel(q, theta=theta, phi=phi)
97    theta = np.linspace(THETA_LOW, THETA_HIGH, n)
98    f_phi = [romberg(inner, PHI_LOW, PHI_HIGH, divmax=100, args=(t,))
99             for t in theta]
100    ans = simps(sin(theta)*np.array(f_phi), dx=theta[1]-theta[0])
101    print("semi romberg evals =", evals[0])
102    return ans*SCALE/(4*pi)
103
104def gauss_quad(q, n=150):
105    """
106    Compute the integral using gaussian quadrature for n = 20, 76 or 150.
107    """
108    if n == 20:
109        z, w = Gauss20Z, Gauss20Wt
110    elif n == 76:
111        z, w = Gauss76Z, Gauss76Wt
112    else:
113        z, w = Gauss150Z, Gauss150Wt
114    theta = (THETA_HIGH-THETA_LOW)*(z + 1)/2 + THETA_LOW
115    phi = (PHI_HIGH-PHI_LOW)*(z + 1)/2 + PHI_LOW
116    Atheta, Aphi = np.meshgrid(theta, phi)
117    Aw = w[None, :] * w[:, None]
118    sin_theta = np.fmax(abs(sin(Atheta)), 1e-6)
119    Zq = kernel(q=q, theta=Atheta, phi=Aphi)
120    print("gauss %d evals ="%n, n**2)
121    return np.sum(Zq*Aw*sin_theta)*SCALE/(4*pi)
122
123
124def gridded_integrals(q, n=300):
125    """
126    Compute the integral on a regular grid using rectangular, trapezoidal,
127    simpsons, and romberg integration.  Romberg integration requires that
128    the grid be of size n = 2**k + 1.
129    """
130    theta = np.linspace(THETA_LOW, THETA_HIGH, n)
131    phi = np.linspace(PHI_LOW, PHI_HIGH, n)
132    Atheta, Aphi = np.meshgrid(theta, phi)
133    Zq = kernel(q=q, theta=Atheta, phi=Aphi)
134    Zq *= abs(sin(Atheta))
135    dx, dy = theta[1]-theta[0], phi[1]-phi[0]
136    print("rect", n, np.sum(Zq)*dx*dy*SCALE/(4*pi))
137    print("trapz", n, np.trapz(np.trapz(Zq, dx=dx), dx=dy)*SCALE/(4*pi))
138    print("simpson", n, simps(simps(Zq, dx=dx), dx=dy)*SCALE/(4*pi))
139    print("romb", n, romb(romb(Zq, dx=dx), dx=dy)*SCALE/(4*pi))
140    print("gridded %d evals ="%n, n**2)
141
142def plot(q, n=300):
143    """
144    Plot the 2D surface that needs to be integrated in order to compute
145    the BCC S(q) at a particular q, dnn and d_factor.  *n* is the number
146    of points in the grid.
147    """
148    theta = np.linspace(THETA_LOW, THETA_HIGH, n)
149    phi = np.linspace(PHI_LOW, PHI_HIGH, n)
150    Atheta, Aphi = np.meshgrid(theta, phi)
151    Zq = kernel(q=q, theta=Atheta, phi=Aphi)
152    #Zq *= abs(sin(Atheta))
153    pylab.pcolor(degrees(theta), degrees(phi), log10(np.fmax(Zq, 1.e-6)))
154    pylab.axis('tight')
155    pylab.title("%s Z(q) for q=%g" % (KERNEL.__name__, q))
156    pylab.xlabel("theta (degrees)")
157    pylab.ylabel("phi (degrees)")
158    cbar = pylab.colorbar()
159    cbar.set_label('log10 S(q)')
160    pylab.show()
161
162if __name__ == "__main__":
163    Q = 0.8
164    print("gauss", 20, gauss_quad(Q, n=20))
165    print("gauss", 76, gauss_quad(Q, n=76))
166    print("gauss", 150, gauss_quad(Q, n=150))
167    #print("dblquad", scipy_dblquad(Q))
168    #print("semi romberg", semi_romberg(Q))
169    #gridded_integrals(Q, n=2**8+1)
170    #gridded_integrals(Q, n=2**10+1)
171    #gridded_integrals(Q, n=2**13+1)
172    #print("romberg", scipy_romberg(Q))
173    plot(Q, n=400)
Note: See TracBrowser for help on using the repository browser.