source: sasmodels/explore/symint.py @ ed5b109

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since ed5b109 was a7db3c05, checked in by dirk, 7 years ago

Fix broken explore/symint.py (calling of gauss_quad_1d)

  • Property mode set to 100644
File size: 5.1 KB
Line 
1"""
2Explore integration of rotationally symmetric shapes
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 numpy.polynomial.legendre import leggauss
13from scipy.integrate import dblquad, simps, romb, romberg
14import pylab
15
16from sasmodels.special import square
17from sasmodels.special import Gauss20Wt, Gauss20Z
18from sasmodels.special import Gauss76Wt, Gauss76Z
19from sasmodels.special import Gauss150Wt, Gauss150Z
20from sasmodels.special import sas_2J1x_x, sas_sinx_x, sas_3j1x_x
21
22SLD = 3.0
23SLD_SOLVENT = 6
24CONTRAST = SLD - SLD_SOLVENT
25
26def make_cylinder(radius, length):
27    def cylinder(qab, qc):
28        return sas_2J1x_x(qab*radius) * sas_sinx_x(qc*0.5*length)
29    cylinder.__doc__ = "cylinder radius=%g, length=%g"%(radius, length)
30    volume = pi*radius**2*length
31    norm = CONTRAST**2*volume/10000
32    return norm, cylinder
33
34def make_long_cylinder(radius, length):
35    def long_cylinder(q):
36        return norm/q * sas_2J1x_x(q*radius)**2
37    long_cylinder.__doc__ = "long cylinder radius=%g, length=%g"%(radius, length)
38    volume = pi*radius**2*length
39    norm = CONTRAST**2*volume/10000*pi/length
40    return long_cylinder
41
42def make_sphere(radius):
43    def sphere(qab, qc):
44        q = sqrt(qab**2 + qc**2)
45        return sas_3j1x_x(q*radius)
46    sphere.__doc__ = "sphere radius=%g"%(radius,)
47    volume = 4*pi*radius**3/3
48    norm = CONTRAST**2*volume/10000
49    return norm, sphere
50
51THETA_LOW, THETA_HIGH = 0, pi/2
52SCALE = 1
53
54
55def kernel_1d(q, theta):
56    """
57    S(q) kernel for paracrystal forms.
58    """
59    qab = q*sin(theta)
60    qc = q*cos(theta)
61    return NORM*KERNEL(qab, qc)**2
62
63def gauss_quad_1d(q, n=150):
64    """
65    Compute the integral using gaussian quadrature for n = 20, 76 or 150.
66    """
67    z, w = leggauss(n)
68    theta = (THETA_HIGH-THETA_LOW)*(z + 1)/2 + THETA_LOW
69    sin_theta = abs(sin(theta))
70    Zq = kernel_1d(q=q, theta=theta)
71    return np.sum(Zq*w*sin_theta)*(THETA_HIGH-THETA_LOW)/2
72
73def gridded_1d(q, n=300):
74    """
75    Compute the integral on a regular grid using rectangular, trapezoidal,
76    simpsons, and romberg integration.  Romberg integration requires that
77    the grid be of size n = 2**k + 1.
78    """
79    theta = np.linspace(THETA_LOW, THETA_HIGH, n)
80    Zq = kernel_1d(q=q, theta=theta)
81    Zq *= abs(sin(theta))
82    dx = theta[1]-theta[0]
83    print("rect-%d"%n, np.sum(Zq)*dx*SCALE)
84    print("trapz-%d"%n, np.trapz(Zq, dx=dx)*SCALE)
85    print("simpson-%d"%n, simps(Zq, dx=dx)*SCALE)
86    print("romb-%d"%n, romb(Zq, dx=dx)*SCALE)
87
88def scipy_romberg_1d(q):
89    """
90    Compute the integral using romberg integration.  This function does not
91    complete in a reasonable time.  No idea if it is accurate.
92    """
93    evals = [0]
94    def outer(theta):
95        evals[0] += 1
96        return kernel_1d(q, theta=theta)*abs(sin(theta))
97    result = romberg(outer, THETA_LOW, THETA_HIGH, divmax=100)*SCALE
98    print("scipy romberg", evals[0], result)
99
100def plot_1d(q, n=300):
101    """
102    Plot the function that needs to be integrated in order to compute
103    the I(q) at a particular q.  *n* is the number of points in the grid.
104    """
105    theta = np.linspace(THETA_LOW, THETA_HIGH, n)
106    Zq = kernel_1d(q=q, theta=theta)
107    Zq *= abs(sin(theta))
108    pylab.semilogy(degrees(theta), np.fmax(Zq, 1.e-6), label="Q=%g"%q)
109    pylab.title("%s I(q, theta) sin(theta)" % (KERNEL.__doc__,))
110    pylab.xlabel("theta (degrees)")
111    pylab.ylabel("Iq 1/cm")
112
113def Iq_trapz(q, n):
114    theta = np.linspace(THETA_LOW, THETA_HIGH, n)
115    Zq = kernel_1d(q=q, theta=theta)
116    Zq *= abs(sin(theta))
117    dx = theta[1]-theta[0]
118    return np.trapz(Zq, dx=dx)*SCALE
119
120def plot_Iq(q, n, form="trapz"):
121    if form == "trapz":
122        Iq = np.array([Iq_trapz(qk, n) for qk in q])
123    elif form == "gauss":
124        Iq = np.array([gauss_quad_1d(qk, n) for qk in q])
125    pylab.loglog(q, Iq, label="%s, n=%d"%(form, n))
126    pylab.xlabel("q (1/A)")
127    pylab.ylabel("Iq (1/cm)")
128    pylab.title(KERNEL.__doc__ + " I(q) circular average")
129    return Iq
130
131radius = 10.
132length = 1e5
133NORM, KERNEL = make_cylinder(radius=radius, length=length)
134long_cyl = make_long_cylinder(radius=radius, length=length)
135#NORM, KERNEL = make_sphere(radius=50.)
136
137
138if __name__ == "__main__":
139    Q = 0.386
140    for n in (20, 76, 150, 300, 1000): #, 10000, 30000):
141        print("gauss-%d"%n, gauss_quad_1d(Q, n=n))
142    for k in (8, 10, 13, 16, 19):
143        gridded_1d(Q, n=2**k+1)
144    #print("inf cyl", 0, long_cyl(Q))
145    #scipy_romberg(Q)
146
147    plot_1d(0.386, n=2000)
148    plot_1d(0.5, n=2000)
149    plot_1d(0.8, n=2000)
150    pylab.legend()
151    pylab.figure()
152
153    q = np.logspace(-3, 0, 400)
154    I1 = long_cyl(q)
155    I2 = plot_Iq(q, n=2**19+1, form="trapz")
156    #plot_Iq(q, n=2**16+1, form="trapz")
157    #plot_Iq(q, n=2**10+1, form="trapz")
158    plot_Iq(q, n=1024, form="gauss")
159    #plot_Iq(q, n=300, form="gauss")
160    #plot_Iq(q, n=150, form="gauss")
161    #plot_Iq(q, n=76, form="gauss")
162    pylab.loglog(q, long_cyl(q), label="limit")
163    pylab.legend()
164
165    pylab.figure()
166    pylab.semilogx(q, (I2 - I1)/I1)
167
168    pylab.show()
Note: See TracBrowser for help on using the repository browser.