source: sasmodels/explore/symint.py @ acc6c8f

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

explore: code cleanup in long cylinder accuracy assessment

  • 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 scipy.integrate import dblquad, simps, romb, romberg
13from scipy.special.orthogonal import p_roots
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.3
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 = 1e-4*volume*CONTRAST**2
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 = 1e-4*volume*CONTRAST**2*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 = 1e-4*volume*CONTRAST**2
49    return norm, sphere
50
51THETA_LOW, THETA_HIGH = 0, pi/2
52SCALE = 1
53
54
55def kernel(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
63
64def gauss_quad(q, n=150):
65    """
66    Compute the integral using gaussian quadrature for n = 20, 76 or 150.
67    """
68    z, w = p_roots(n)
69    theta = (THETA_HIGH-THETA_LOW)*(z + 1)/2 + THETA_LOW
70    sin_theta = abs(sin(theta))
71    Zq = kernel(q=q, theta=theta)
72    return np.sum(Zq*w*sin_theta)*(THETA_HIGH-THETA_LOW)/2
73
74
75def gridded_integrals(q, n=300):
76    """
77    Compute the integral on a regular grid using rectangular, trapezoidal,
78    simpsons, and romberg integration.  Romberg integration requires that
79    the grid be of size n = 2**k + 1.
80    """
81    theta = np.linspace(THETA_LOW, THETA_HIGH, n)
82    Zq = kernel(q=q, theta=theta)
83    Zq *= abs(sin(theta))
84    dx = theta[1]-theta[0]
85    print("rect", n, np.sum(Zq)*dx*SCALE)
86    print("trapz", n, np.trapz(Zq, dx=dx)*SCALE)
87    print("simpson", n, simps(Zq, dx=dx)*SCALE)
88    print("romb", n, romb(Zq, dx=dx)*SCALE)
89
90def scipy_romberg(q):
91    """
92    Compute the integral using romberg integration.  This function does not
93    complete in a reasonable time.  No idea if it is accurate.
94    """
95    evals = [0]
96    def outer(theta):
97        evals[0] += 1
98        return kernel(q, theta=theta)*abs(sin(theta))
99    result = romberg(outer, THETA_LOW, THETA_HIGH, divmax=100)*SCALE
100    print("scipy romberg", evals[0], result)
101
102def plot(q, n=300):
103    """
104    Plot the 2D surface that needs to be integrated in order to compute
105    the BCC S(q) at a particular q, dnn and d_factor.  *n* is the number
106    of points in the grid.
107    """
108    theta = np.linspace(THETA_LOW, THETA_HIGH, n)
109    Zq = kernel(q=q, theta=theta)
110    Zq *= abs(sin(theta))
111    pylab.semilogy(degrees(theta), np.fmax(Zq, 1.e-6), label="Q=%g"%q)
112    pylab.title("%s I(q, theta) sin(theta)" % (KERNEL.__doc__,))
113    pylab.xlabel("theta (degrees)")
114    pylab.ylabel("Iq 1/cm")
115
116def Iq_trapz(q, n):
117    theta = np.linspace(THETA_LOW, THETA_HIGH, n)
118    Zq = kernel(q=q, theta=theta)
119    Zq *= abs(sin(theta))
120    dx = theta[1]-theta[0]
121    return np.trapz(Zq, dx=dx)*SCALE
122
123def plot_Iq(q, n, form="trapz"):
124    if form == "trapz":
125        Iq = np.array([Iq_trapz(qk, n) for qk in q])
126    elif form == "gauss":
127        Iq = np.array([gauss_quad(qk, n) for qk in q])
128    pylab.loglog(q, Iq, label="%s, n=%d"%(form, n))
129    pylab.xlabel("q (1/A)")
130    pylab.ylabel("Iq (1/cm)")
131    pylab.title(KERNEL.__doc__ + " I(q) circular average")
132    return Iq
133
134radius = 10.
135length = 1e5
136NORM, KERNEL = make_cylinder(radius=radius, length=length)
137long_cyl = make_long_cylinder(radius=radius, length=length)
138#NORM, KERNEL = make_sphere(radius=50.)
139
140
141if __name__ == "__main__":
142    Q = 0.386
143    for n in (20, 76, 150, 300, 1000): #, 10000, 30000):
144        print("gauss", n, gauss_quad(Q, n=n))
145    for k in (8, 10, 13, 16, 19):
146        gridded_integrals(Q, n=2**k+1)
147    #print("inf cyl", 0, long_cyl(Q))
148    #scipy_romberg(Q)
149
150    plot(0.386, n=2000)
151    plot(0.5, n=2000)
152    plot(0.8, n=2000)
153    pylab.legend()
154    pylab.figure()
155
156    q = np.logspace(-3, 0, 400)
157    I1 = long_cyl(q)
158    I2 = plot_Iq(q, n=2**19+1, form="trapz")
159    #plot_Iq(q, n=2**16+1, form="trapz")
160    #plot_Iq(q, n=2**10+1, form="trapz")
161    plot_Iq(q, n=1024, form="gauss")
162    #plot_Iq(q, n=300, form="gauss")
163    #plot_Iq(q, n=150, form="gauss")
164    #plot_Iq(q, n=76, form="gauss")
165    pylab.loglog(q, long_cyl(q), label="limit")
166    pylab.legend()
167
168    pylab.figure()
169    pylab.semilogx(q, (I2 - I1)/I1)
170
171    pylab.show()
Note: See TracBrowser for help on using the repository browser.