Changeset 20fe0cd in sasmodels for explore/symint.py

Timestamp:

Oct 23, 2017 9:17:40 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:

6e604f8

Parents:

36b3154

Message:

move paracrystal integration tests with the rest of the non-rotationally symmetric tests

File:

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

explore/symint.py

@@ -10 +10 @@
 import numpy as np
 from numpy import pi, sin, cos, sqrt, exp, expm1, degrees, log10
+from numpy.polynomial.legendre import leggauss
 from scipy.integrate import dblquad, simps, romb, romberg
-from scipy.special.orthogonal import p_roots
 import pylab
 
@@ -21 +21 @@
 
 SLD = 3.0
-SLD_SOLVENT = 6.3
+SLD_SOLVENT = 6
 CONTRAST = SLD - SLD_SOLVENT
 
@@ -29 +29 @@
     cylinder.__doc__ = "cylinder radius=%g, length=%g"%(radius, length)
     volume = pi*radius**2*length
-    norm = 1e-4*volume*CONTRAST**2
+    norm = CONTRAST**2*volume/10000
     return norm, cylinder
 
@@ -37 +37 @@
     long_cylinder.__doc__ = "long cylinder radius=%g, length=%g"%(radius, length)
     volume = pi*radius**2*length
-    norm = 1e-4*volume*CONTRAST**2*pi/length
+    norm = CONTRAST**2*volume/10000*pi/length
     return long_cylinder
 
@@ -46 +46 @@
     sphere.__doc__ = "sphere radius=%g"%(radius,)
     volume = 4*pi*radius**3/3
-    norm = 1e-4*volume*CONTRAST**2
+    norm = CONTRAST**2*volume/10000
     return norm, sphere
 
@@ -53 +53 @@
 
 
-def kernel(q, theta):
+def kernel_1d(q, theta):
     """
     S(q) kernel for paracrystal forms.
@@ -61 +61 @@
     return NORM*KERNEL(qab, qc)**2
 
-
-def gauss_quad(q, n=150):
+def gauss_quad_1d(q, n=150):
     """
     Compute the integral using gaussian quadrature for n = 20, 76 or 150.
     """
-    z, w = p_roots(n)
+    z, w = leggauss(n)
     theta = (THETA_HIGH-THETA_LOW)*(z + 1)/2 + THETA_LOW
     sin_theta = abs(sin(theta))
-    Zq = kernel(q=q, theta=theta)
+    Zq = kernel_1d(q=q, theta=theta)
     return np.sum(Zq*w*sin_theta)*(THETA_HIGH-THETA_LOW)/2
 
-
-def gridded_integrals(q, n=300):
+def gridded_1d(q, n=300):
     """
     Compute the integral on a regular grid using rectangular, trapezoidal,
@@ -80 +78 @@
     """
     theta = np.linspace(THETA_LOW, THETA_HIGH, n)
-    Zq = kernel(q=q, theta=theta)
+    Zq = kernel_1d(q=q, theta=theta)
     Zq *= abs(sin(theta))
     dx = theta[1]-theta[0]
-    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)
+    print("rect-%d"%n, np.sum(Zq)*dx*SCALE)
+    print("trapz-%d"%n, np.trapz(Zq, dx=dx)*SCALE)
+    print("simpson-%d"%n, simps(Zq, dx=dx)*SCALE)
+    print("romb-%d"%n, romb(Zq, dx=dx)*SCALE)
 
-def scipy_romberg(q):
+def scipy_romberg_1d(q):
     """
     Compute the integral using romberg integration.  This function does not
@@ -96 +94 @@
     def outer(theta):
         evals[0] += 1
-        return kernel(q, theta=theta)*abs(sin(theta))
+        return kernel_1d(q, theta=theta)*abs(sin(theta))
     result = romberg(outer, THETA_LOW, THETA_HIGH, divmax=100)*SCALE
     print("scipy romberg", evals[0], result)
 
-def plot(q, n=300):
+def plot_1d(q, n=300):
     """
-    Plot the 2D surface that needs to be integrated in order to compute
-    the BCC S(q) at a particular q, dnn and d_factor.  *n* is the number
-    of points in the grid.
+    Plot the function that needs to be integrated in order to compute
+    the I(q) at a particular q.  *n* is the number of points in the grid.
     """
     theta = np.linspace(THETA_LOW, THETA_HIGH, n)
-    Zq = kernel(q=q, theta=theta)
+    Zq = kernel_1d(q=q, theta=theta)
     Zq *= abs(sin(theta))
     pylab.semilogy(degrees(theta), np.fmax(Zq, 1.e-6), label="Q=%g"%q)
@@ -116 +113 @@
 def Iq_trapz(q, n):
     theta = np.linspace(THETA_LOW, THETA_HIGH, n)
-    Zq = kernel(q=q, theta=theta)
+    Zq = kernel_1d(q=q, theta=theta)
     Zq *= abs(sin(theta))
     dx = theta[1]-theta[0]
@@ -142 +139 @@
     Q = 0.386
     for n in (20, 76, 150, 300, 1000): #, 10000, 30000):
-        print("gauss", n, gauss_quad(Q, n=n))
+        print("gauss-%d"%n, gauss_quad(Q, n=n))
     for k in (8, 10, 13, 16, 19):
-        gridded_integrals(Q, n=2**k+1)
+        gridded_1d(Q, n=2**k+1)
     #print("inf cyl", 0, long_cyl(Q))
     #scipy_romberg(Q)
 
-    plot(0.386, n=2000)
-    plot(0.5, n=2000)
-    plot(0.8, n=2000)
+    plot_1d(0.386, n=2000)
+    plot_1d(0.5, n=2000)
+    plot_1d(0.8, n=2000)
     pylab.legend()
     pylab.figure()