Changeset d327066 in sasmodels

Timestamp:

Jun 29, 2017 12:42:23 PM (7 years ago)

Author:

Paul Kienzle <pkienzle@…>

Branches:

master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

fac23b3

Parents:

a2fcbd8

Message:

play with cylinder integral

File:

• explore/symint.py (modified) (11 diffs)

explore/symint.py

@@ -11 +11 @@
 from numpy import pi, sin, cos, sqrt, exp, expm1, degrees, log10
 from scipy.integrate import dblquad, simps, romb, romberg
+from scipy.special.orthogonal import p_roots
 import pylab
 
@@ -26 +27 @@
     def cylinder(qab, qc):
         return sas_2J1x_x(qab*radius) * sas_sinx_x(qc*0.5*length)
+    cylinder.__doc__ = "cylinder radius=%g, length=%g"%(radius, length)
     volume = pi*radius**2*length
     norm = 1e-4*volume*CONTRAST**2
@@ -34 +36 @@
         q = sqrt(qab**2 + qc**2)
         return sas_3j1x_x(q*radius)
+    sphere.__doc__ = "sphere radius=%g"%(radius,)
     volume = 4*pi*radius**3/3
     norm = 1e-4*volume*CONTRAST**2
-    return norm, sphere 
-    
-THETA_LOW, THETA_HIGH = 0, pi
+    return norm, sphere
+
+THETA_LOW, THETA_HIGH = 0, pi/2
 SCALE = 1
 
@@ -55 +58 @@
     Compute the integral using gaussian quadrature for n = 20, 76 or 150.
     """
-    if n == 20:
-        z, w = Gauss20Z, Gauss20Wt
-    elif n == 76:
-        z, w = Gauss76Z, Gauss76Wt
-    else:
-        z, w = Gauss150Z, Gauss150Wt
+    z, w = p_roots(n)
     theta = (THETA_HIGH-THETA_LOW)*(z + 1)/2 + THETA_LOW
     sin_theta = abs(sin(theta))
     Zq = kernel(q=q, theta=theta)
-    return np.sum(Zq*w*sin_theta)*SCALE/2
+    return np.sum(Zq*w*sin_theta)*(THETA_HIGH-THETA_LOW)/2
 
 
@@ -77 +75 @@
     Zq *= abs(sin(theta))
     dx = theta[1]-theta[0]
-    print("rect", n, np.sum(Zq)*dx*SCALE/pi)
-    print("trapz", n, np.trapz(Zq, dx=dx)*SCALE/pi)
-    print("simpson", n, simps(Zq, dx=dx)*SCALE/pi)
-    print("romb", n, romb(Zq, dx=dx)*SCALE/pi)
+    print("rect", n, np.sum(Zq)*dx*SCALE)
+    print("trapz", n, np.trapz(Zq, dx=dx)*SCALE)
+    print("simpson", n, simps(Zq, dx=dx)*SCALE)
+    print("romb", n, romb(Zq, dx=dx)*SCALE)
 
 def scipy_romberg(q):
@@ -91 +89 @@
         evals[0] += 1
         return kernel(q, theta=theta)*abs(sin(theta))
-    result = romberg(outer, THETA_LOW, THETA_HIGH, divmax=100)*SCALE/pi
+    result = romberg(outer, THETA_LOW, THETA_HIGH, divmax=100)*SCALE
     print("scipy romberg", evals[0], result)
 
@@ -104 +102 @@
     Zq *= abs(sin(theta))
     pylab.semilogy(degrees(theta), np.fmax(Zq, 1.e-6), label="Q=%g"%q)
-    pylab.title("%s I(q, theta) sin(theta)" % (KERNEL.__name__,))
+    pylab.title("%s I(q, theta) sin(theta)" % (KERNEL.__doc__,))
     pylab.xlabel("theta (degrees)")
     pylab.ylabel("Iq 1/cm")
@@ -113 +111 @@
     Zq *= abs(sin(theta))
     dx = theta[1]-theta[0]
-    return np.trapz(Zq, dx=dx)*SCALE/pi
+    return np.trapz(Zq, dx=dx)*SCALE
 
 def plot_Iq(q, n, form="trapz"):
@@ -123 +121 @@
     pylab.xlabel("q (1/A)")
     pylab.ylabel("Iq (1/cm)")
+    pylab.title(KERNEL.__doc__ + " I(q) circular average")
+    return q, I
 
 NORM, KERNEL = make_cylinder(radius=10., length=100000.)
+#NORM, KERNEL = make_cylinder(radius=10., length=50000.)
+#NORM, KERNEL = make_cylinder(radius=10., length=20000.)
 #NORM, KERNEL = make_cylinder(radius=10., length=10000.)
+#NORM, KERNEL = make_cylinder(radius=10., length=5000.)
+#NORM, KERNEL = make_cylinder(radius=10., length=1000.)
+#NORM, KERNEL = make_cylinder(radius=10., length=500.)
+#NORM, KERNEL = make_cylinder(radius=10., length=100.)
 #NORM, KERNEL = make_cylinder(radius=10., length=30.)
 #NORM, KERNEL = make_sphere(radius=50.)
@@ -131 +137 @@
 if __name__ == "__main__":
     Q = 0.8
-    print("gauss", 20, gauss_quad(Q, n=20))
-    print("gauss", 76, gauss_quad(Q, n=76))
-    print("gauss", 150, gauss_quad(Q, n=150))
-    gridded_integrals(Q, n=2**8+1)
-    gridded_integrals(Q, n=2**10+1)
-    gridded_integrals(Q, n=2**13+1)
-    gridded_integrals(Q, n=2**16+1)
-    gridded_integrals(Q, n=2**19+1)
+    for n in (20, 76, 150, 300, 1000): #, 10000, 30000):
+        print("gauss", n, gauss_quad(Q, n=n))
+    for k in (8, 10, 13, 16, 19):
+        gridded_integrals(Q, n=2**k+1)
     #scipy_romberg(Q)
     plot(0.5, n=2000)
@@ -145 +147 @@
     pylab.legend()
     pylab.figure()
-    plot_Iq(np.logspace(-3,0,200), n=2**16+1, form="trapz")
-    plot_Iq(np.logspace(-3,0,200), n=2**10+1, form="trapz")
-    plot_Iq(np.logspace(-3,0,200), n=150, form="gauss")
+    #plot_Iq(np.logspace(-3, 0, 400), n=2**19+1, form="trapz")
+    q1, I1 = plot_Iq(np.logspace(-3, 0, 400), n=2**16+1, form="trapz")
+    #plot_Iq(np.logspace(-3, 0, 400), n=2**10+1, form="trapz")
+    q2, I2 = plot_Iq(np.logspace(-3, 0, 400), n=1000, form="gauss")
+    #plot_Iq(np.logspace(-3, 0, 400), n=300, form="gauss")
+    plot_Iq(np.logspace(-3, 0, 400), n=150, form="gauss")
+    plot_Iq(np.logspace(-3, 0, 400), n=76, form="gauss")
     pylab.legend()
+    #pylab.figure()
+    #pylab.semilogx(q1, (I2 - I1)/I1)
     pylab.show()