Changeset 20fe0cd in sasmodels for explore

Timestamp:

Oct 23, 2017 11: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

Location:

explore

Files:

• angles.py (modified) (1 diff, 1 prop)
• asymint.py (modified) (18 diffs, 1 prop)
• bccpy.py (deleted)
• symint.py (modified) (11 diffs)
• transform_angles.py (added)

explore/angles.py

@@ -229 +229 @@
 
 if 0:
-    comment("==== asymmetric (old) ====")
+    comment("==== asymmetric (3.x) ====")
     view_inv = Rz(-psi)*Rx(theta)*Ry(-(pi/2 - phi))
     vprint(q_abc, view_inv*q_xy, "reverse view")
+    print("""  existing code
+    cos_alpha = cos_theta*cos_phi*qxhat + sin_theta*qyhat;
+    cos_mu = (-sin_theta*cos_psi*cos_phi - sin_psi*sin_phi)*qxhat + cos_theta*cos_psi*qyhat;
+    cos_nu = (-cos_phi*sin_psi*sin_theta + sin_phi*cos_psi)*qxhat + sin_psi*cos_theta*qyhat;
+    """)

explore/asymint.py

@@ +1 @@
+#!/usr/bin/env python
 """
 Asymmetric shape integration
@@ -6 +7 @@
 
 import os, sys
-sys.path.insert(0, os.path.dirname(os.path.dirname(__file__)))
+sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.realpath(__file__))))
 
 import numpy as np
@@ -15 +16 @@
 import pylab
 
-from sasmodels.special import square
-from sasmodels.special import Gauss20Wt, Gauss20Z
-from sasmodels.special import Gauss76Wt, Gauss76Z
-from sasmodels.special import Gauss150Wt, Gauss150Z
-from sasmodels.special import sas_2J1x_x, sas_sinx_x, sas_3j1x_x
-
-def mp_3j1x_x(x):
-    return 3*(mp.sin(x)/x - mp.cos(x))/(x*x)
-def mp_2J1x_x(x):
-    return 2*mp.j1(x)/x
-def mp_sinx_x(x):
-    return mp.sin(x)/x
+import sasmodels.special as sp
+
+# Need to parse shape early since it determines the kernel function
+# that will be used for the various integrators
+shape = 'parallelepiped' if len(sys.argv) < 2 else sys.argv[1]
+Qstr = '0.005' if len(sys.argv) < 3 else sys.argv[2]
+
+class MPenv:
+    sqrt = mp.sqrt
+    exp = mp.exp
+    expm1 = mp.expm1
+    cos = mp.cos
+    sin = mp.sin
+    tan = mp.tan
+    def sas_3j1x_x(x):
+        return 3*(mp.sin(x)/x - mp.cos(x))/(x*x)
+    def sas_2J1x_x(x):
+        return 2*mp.j1(x)/x
+    def sas_sinx_x(x):
+        return mp.sin(x)/x
+    pi = mp.pi
+    mpf = mp.mpf
+
+class NPenv:
+    sqrt = np.sqrt
+    exp = np.exp
+    expm1 = np.expm1
+    cos = np.cos
+    sin = np.sin
+    tan = np.tan
+    sas_3j1x_x = sp.sas_3j1x_x
+    sas_2J1x_x = sp.sas_2J1x_x
+    sas_sinx_x = sp.sas_sinx_x
+    pi = np.pi
+    mpf = float
 
 SLD = 3
@@ -32 +56 @@
 CONTRAST = SLD - SLD_SOLVENT
 
-def make_parallelepiped(a, b, c):
-    def Fq(qa, qb, qc):
-        siA = sas_sinx_x(0.5*a*qa)
-        siB = sas_sinx_x(0.5*b*qb)
-        siC = sas_sinx_x(0.5*c*qc)
+# Carefully code models so that mpmath will use full precision.  That means:
+#    * wrap inputs in env.mpf
+#    * don't use floating point constants, only integers
+#    * for division, make sure the numerator or denominator is env.mpf
+#    * use env.pi, env.sas_sinx_x, etc. for functions
+def make_parallelepiped(a, b, c, env=NPenv):
+    a, b, c = env.mpf(a), env.mpf(b), env.mpf(c)
+    def Fq(qa, qb, qc):
+        siA = env.sas_sinx_x(0.5*a*qa/2)
+        siB = env.sas_sinx_x(0.5*b*qb/2)
+        siC = env.sas_sinx_x(0.5*c*qc/2)
         return siA * siB * siC
+    Fq.__doc__ = "parallelepiped a=%g, b=%g c=%g"%(a, b, c)
     volume = a*b*c
-    norm = volume*CONTRAST**2/10**4
-    return norm, Fq
-
-def make_parallelepiped_mp(a, b, c):
-    a, b, c = mp.mpf(a), mp.mpf(b), mp.mpf(c)
-    def Fq(qa, qb, qc):
-        siA = mp_sinx_x(a*qa/2)
-        siB = mp_sinx_x(b*qb/2)
-        siC = mp_sinx_x(c*qc/2)
-        return siA * siB * siC
-    volume = a*b*c
-    norm = (volume*CONTRAST**2)/10000 # mpf since volume=a*b*c is mpf
-    return norm, Fq
-
-def make_triellip(a, b, c):
-    def Fq(qa, qb, qc):
-        qr = sqrt((a*qa)**2 + (b*qb)**2 + (c*qc)**2)
-        return sas_3j1x_x(qr)
-    volume = 4*pi*a*b*c/3
-    norm = volume*CONTRAST**2/10**4
-    return norm, Fq
-
-def make_triellip_mp(a, b, c):
-    a, b, c = mp.mpf(a), mp.mpf(b), mp.mpf(c)
-    def Fq(qa, qb, qc):
-        qr = mp.sqrt((a*qa)**2 + (b*qb)**2 + (c*qc)**2)
-        return mp_3j1x_x(qr)
-    volume = (4*mp.pi*a*b*c)/3
-    norm = (volume*CONTRAST**2)/10000  # mpf since mp.pi is mpf
-    return norm, Fq
-
-def make_cylinder(radius, length):
-    def Fq(qa, qb, qc):
-        qab = sqrt(qa**2 + qb**2)
-        return sas_2J1x_x(qab*radius) * sas_sinx_x(qc*0.5*length)
-    volume = pi*radius**2*length
-    norm = volume*CONTRAST**2/10**4
-    return norm, Fq
-
-def make_cylinder_mp(radius, length):
-    radius, length = mp.mpf(radius), mp.mpf(length)
-    def Fq(qa, qb, qc):
-        qab = mp.sqrt(qa**2 + qb**2)
-        return mp_2J1x_x(qab*radius) * mp_sinx_x((qc*length)/2)
-    volume = mp.pi*radius**2*length
-    norm = (volume*CONTRAST**2)/10000  # mpf since mp.pi is mpf
-    return norm, Fq
-
-def make_sphere(radius):
-    def Fq(qa, qb, qc):
-        q = sqrt(qa**2 + qb**2 + qc**2)
-        return sas_3j1x_x(q*radius)
-    volume = 4*pi*radius**3/3
-    norm = volume*CONTRAST**2/10**4
-    return norm, Fq
-
-def make_sphere_mp(radius):
-    radius = mp.mpf(radius)
-    def Fq(qa, qb, qc):
-        q = mp.sqrt(qa**2 + qb**2 + qc**2)
-        return mp_3j1x_x(q*radius)
-    volume = (4*mp.pi*radius**3)/3
-    norm = (volume*CONTRAST**2)/10000  # mpf since mp.pi is mpf
-    return norm, Fq
-
-shape = 'parallelepiped'
-#shape = 'triellip'
-#shape = 'sphere'
-#shape = 'cylinder'
-if shape == 'cylinder':
+    norm = CONTRAST**2*volume/10000
+    return norm, Fq
+
+def make_triaxial_ellipsoid(a, b, c, env=NPenv):
+    a, b, c = env.mpf(a), env.mpf(b), env.mpf(c)
+    def Fq(qa, qb, qc):
+        qr = env.sqrt((a*qa)**2 + (b*qb)**2 + (c*qc)**2)
+        return env.sas_3j1x_x(qr)
+    Fq.__doc__ = "triaxial ellipse minor=%g, major=%g polar=%g"%(a, b, c)
+    volume = 4*env.pi*a*b*c/3
+    norm = CONTRAST**2*volume/10000
+    return norm, Fq
+
+def make_cylinder(radius, length, env=NPenv):
+    radius, length = env.mpf(radius), env.mpf(length)
+    def Fq(qa, qb, qc):
+        qab = env.sqrt(qa**2 + qb**2)
+        return env.sas_2J1x_x(qab*radius) * env.sas_sinx_x((qc*length)/2)
+    Fq.__doc__ = "cylinder radius=%g, length=%g"%(radius, length)
+    volume = env.pi*radius**2*length
+    norm = CONTRAST**2*volume/10000
+    return norm, Fq
+
+def make_sphere(radius, env=NPenv):
+    radius = env.mpf(radius)
+    def Fq(qa, qb, qc):
+        q = env.sqrt(qa**2 + qb**2 + qc**2)
+        return env.sas_3j1x_x(q*radius)
+    Fq.__doc__ = "sphere radius=%g"%(radius, )
+    volume = 4*pi*radius**3
+    norm = CONTRAST**2*volume/10000
+    return norm, Fq
+
+def make_paracrystal(radius, dnn, d_factor, lattice='bcc', env=NPenv):
+    radius, dnn, d_factor = env.mpf(radius), env.mpf(dnn), env.mpf(d_factor)
+    def sc(qa, qb, qc):
+        return qa, qb, qc
+    def bcc(qa, qb, qc):
+        a1 = (+qa + qb + qc)/2
+        a2 = (-qa - qb + qc)/2
+        a3 = (-qa + qb - qc)/2
+        return a1, a2, a3
+    def fcc(qa, qb, qc):
+        a1 = ( 0. + qb + qc)/2
+        a2 = (-qa + 0. + qc)/2
+        a3 = (-qa + qb + 0.)/2
+        return a1, a2, a3
+    lattice_fn = {'sc': sc, 'bcc': bcc, 'fcc': fcc}[lattice]
+    radius, dnn, d_factor = env.mpf(radius), env.mpf(dnn), env.mpf(d_factor)
+    def Fq(qa, qb, qc):
+        a1, a2, a3 = lattice_fn(qa, qb, qc)
+        # Note: paper says that different directions can have different
+        # distoration factors.  Easy enough to add to the code.
+        # This would definitely break 8-fold symmetry.
+        arg = -(dnn*d_factor)**2*(a1**2 + a2**2 + a3**2)/2
+        exp_arg = env.exp(arg)
+        den = [((exp_arg - 2*env.cos(dnn*a))*exp_arg + 1) for a in (a1, a2, a3)]
+        Sq = -env.expm1(2*arg)**3/(den[0]*den[1]*den[2])
+
+        q = env.sqrt(qa**2 + qb**2 + qc**2)
+        Fq = env.sas_3j1x_x(q*radius)
+        # the kernel computes F(q)**2, but we need S(q)*F(q)**2
+        return env.sqrt(Sq)*Fq
+    Fq.__doc__ = "%s paracrystal a=%g da=%g r=%g"%(lattice, dnn, d_factor, radius)
+    def sphere_volume(r): return 4*env.pi*r**3/3
+    Vf = {
+        'sc': sphere_volume(radius/dnn),
+        'bcc': 2*sphere_volume(env.sqrt(3)/2*radius/dnn),
+        'fcc': 4*sphere_volume(radius/dnn/env.sqrt(2)),
+    }[lattice]
+    volume = sphere_volume(radius)
+    norm = CONTRAST**2*volume*Vf/10000
+    return norm, Fq
+
+if shape == 'sphere':
+    RADIUS = 50  # integer for the sake of mpf
+    NORM, KERNEL = make_sphere(radius=RADIUS)
+    NORM_MP, KERNEL_MP = make_sphere(radius=RADIUS, env=MPenv)
+elif shape == 'cylinder':
     #RADIUS, LENGTH = 10, 100000
     RADIUS, LENGTH = 10, 300  # integer for the sake of mpf
     NORM, KERNEL = make_cylinder(radius=RADIUS, length=LENGTH)
-    NORM_MP, KERNEL_MP = make_cylinder_mp(radius=RADIUS, length=LENGTH)
-elif shape == 'triellip':
+    NORM_MP, KERNEL_MP = make_cylinder(radius=RADIUS, length=LENGTH, env=MPenv)
+elif shape == 'triaxial_ellipsoid':
     #A, B, C = 4450, 14000, 47
     A, B, C = 445, 140, 47  # integer for the sake of mpf
     NORM, KERNEL = make_triellip(A, B, C)
-    NORM_MP, KERNEL_MP = make_triellip_mp(A, B, C)
+    NORM_MP, KERNEL_MP = make_triellip(A, B, C, env=MPenv)
 elif shape == 'parallelepiped':
     #A, B, C = 4450, 14000, 47
     A, B, C = 445, 140, 47  # integer for the sake of mpf
     NORM, KERNEL = make_parallelepiped(A, B, C)
-    NORM_MP, KERNEL_MP = make_parallelepiped_mp(A, B, C)
-elif shape == 'sphere':
-    RADIUS = 50  # integer for the sake of mpf
-    NORM, KERNEL = make_sphere(radius=RADIUS)
-    NORM_MP, KERNEL_MP = make_sphere_mp(radius=RADIUS)
+    NORM_MP, KERNEL_MP = make_parallelepiped(A, B, C, env=MPenv)
+elif shape == 'paracrystal':
+    LATTICE = 'bcc'
+    #LATTICE = 'fcc'
+    #LATTICE = 'sc'
+    DNN, D_FACTOR = 220, '0.06'  # mpmath needs to initialize floats from string
+    RADIUS = 40  # integer for the sake of mpf
+    NORM, KERNEL = make_paracrystal(
+        radius=RADIUS, dnn=DNN, d_factor=D_FACTOR, lattice=LATTICE)
+    NORM_MP, KERNEL_MP = make_paracrystal(
+        radius=RADIUS, dnn=DNN, d_factor=D_FACTOR, lattice=LATTICE, env=MPenv)
 else:
     raise ValueError("Unknown shape %r"%shape)
 
+# Note: hardcoded in mp_quad
 THETA_LOW, THETA_HIGH = 0, pi
 PHI_LOW, PHI_HIGH = 0, 2*pi
@@ -145 +192 @@
 """
 
-
-def mp_quad(q, shape):
+# 2D integration functions
+def mp_quad_2d(q, shape):
     evals = [0]
     def integrand(theta, phi):
@@ -160 +207 @@
     return evals[0], Iq
 
-def kernel(q, theta, phi):
+def kernel_2d(q, theta, phi):
     """
     S(q) kernel for paracrystal forms.
@@ -170 +217 @@
     return NORM*KERNEL(qa, qb, qc)**2
 
-def scipy_dblquad(q):
+def scipy_dblquad_2d(q):
     """
     Compute the integral using scipy dblquad.  This gets the correct answer
@@ -178 +225 @@
     def integrand(phi, theta):
         evals[0] += 1
-        Zq = kernel(q, theta=theta, phi=phi)
+        Zq = kernel_2d(q, theta=theta, phi=phi)
         return Zq*sin(theta)
     ans = dblquad(integrand, THETA_LOW, THETA_HIGH, lambda x: PHI_LOW, lambda x: PHI_HIGH)[0]
@@ -192 +239 @@
     def inner(phi, theta):
         evals[0] += 1
-        return kernel(q, theta=theta, phi=phi)
+        return kernel_2d(q, theta=theta, phi=phi)
     def outer(theta):
         Zq = romberg(inner, PHI_LOW, PHI_HIGH, divmax=100, args=(theta,))
@@ -200 +247 @@
 
 
-def semi_romberg(q, n=100):
+def semi_romberg_2d(q, n=100):
     """
     Use 1D romberg integration in phi and regular simpsons rule in theta.
@@ -207 +254 @@
     def inner(phi, theta):
         evals[0] += 1
-        return kernel(q, theta=theta, phi=phi)
+        return kernel_2d(q, theta=theta, phi=phi)
     theta = np.linspace(THETA_LOW, THETA_HIGH, n)
     Zq = [romberg(inner, PHI_LOW, PHI_HIGH, divmax=100, args=(t,)) for t in theta]
@@ -213 +260 @@
     return evals[0], ans*SCALE/(4*pi)
 
-def gauss_quad(q, n=150):
+def gauss_quad_2d(q, n=150):
     """
     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
-    elif n == 150:
-        z, w = Gauss150Z, Gauss150Wt
-    else:
-        z, w = leggauss(n)
+    z, w = leggauss(n)
     theta = (THETA_HIGH-THETA_LOW)*(z + 1)/2 + THETA_LOW
     phi = (PHI_HIGH-PHI_LOW)*(z + 1)/2 + PHI_LOW
@@ -230 +270 @@
     Aw = w[None, :] * w[:, None]
     sin_theta = abs(sin(Atheta))
-    Zq = kernel(q=q, theta=Atheta, phi=Aphi)
+    Zq = kernel_2d(q=q, theta=Atheta, phi=Aphi)
     # change from [-1,1] x [-1,1] range to [0, pi] x [0, 2 pi] range
     dxdy_stretch = (THETA_HIGH-THETA_LOW)/2 * (PHI_HIGH-PHI_LOW)/2
@@ -236 +276 @@
     return n**2, Iq
 
-
-def gridded_integrals(q, n=300):
+def gridded_2d(q, n=300):
     """
     Compute the integral on a regular grid using rectangular, trapezoidal,
@@ -246 +285 @@
     phi = np.linspace(PHI_LOW, PHI_HIGH, n)
     Atheta, Aphi = np.meshgrid(theta, phi)
-    Zq = kernel(q=q, theta=Atheta, phi=Aphi)
+    Zq = kernel_2d(q=q, theta=Atheta, phi=Aphi)
     Zq *= abs(sin(Atheta))
     dx, dy = theta[1]-theta[0], phi[1]-phi[0]
@@ -254 +293 @@
     print("romb-%d"%n, n**2, romb(romb(Zq, dx=dx), dx=dy)*SCALE/(4*pi))
 
-def plot(q, n=300):
+def plot_2d(q, n=300):
     """
     Plot the 2D surface that needs to be integrated in order to compute
@@ -263 +302 @@
     phi = np.linspace(PHI_LOW, PHI_HIGH, n)
     Atheta, Aphi = np.meshgrid(theta, phi)
-    Zq = kernel(q=q, theta=Atheta, phi=Aphi)
+    Zq = kernel_2d(q=q, theta=Atheta, phi=Aphi)
     #Zq *= abs(sin(Atheta))
     pylab.pcolor(degrees(theta), degrees(phi), log10(np.fmax(Zq, 1.e-6)))
     pylab.axis('tight')
-    pylab.title("%s Z(q) for q=%g" % (KERNEL.__name__, q))
+    pylab.title("%s I(q,t) sin(t) for q=%g" % (KERNEL.__doc__, q))
     pylab.xlabel("theta (degrees)")
     pylab.ylabel("phi (degrees)")
@@ -274 +313 @@
     pylab.show()
 
-if __name__ == "__main__":
-    Qstr = '0.005'
-    #Qstr = '0.8'
-    #Qstr = '0.0003'
+def main(Qstr):
     Q = float(Qstr)
     if shape == 'sphere':
         print("exact", NORM*sas_3j1x_x(Q*RADIUS)**2)
-    print("gauss-20", *gauss_quad(Q, n=20))
-    print("gauss-76", *gauss_quad(Q, n=76))
-    print("gauss-150", *gauss_quad(Q, n=150))
-    print("gauss-500", *gauss_quad(Q, n=500))
-    print("dblquad", *scipy_dblquad(Q))
-    print("semi-romberg-100", *semi_romberg(Q, n=100))
-    print("romberg", *scipy_romberg_2d(Q))
-    #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**15+1)
-    with mp.workprec(100):
-        print("mpmath", *mp_quad(mp.mpf(Qstr), shape))
-    #plot(Q, n=200)
+    print("gauss-20", *gauss_quad_2d(Q, n=20))
+    print("gauss-76", *gauss_quad_2d(Q, n=76))
+    print("gauss-150", *gauss_quad_2d(Q, n=150))
+    print("gauss-500", *gauss_quad_2d(Q, n=500))
+    #gridded_2d(Q, n=2**8+1)
+    gridded_2d(Q, n=2**10+1)
+    #gridded_2d(Q, n=2**13+1)
+    #gridded_2d(Q, n=2**15+1)
+    if shape != 'paracrystal':  # adaptive forms are too slow!
+        print("dblquad", *scipy_dblquad_2d(Q))
+        print("semi-romberg-100", *semi_romberg_2d(Q, n=100))
+        print("romberg", *scipy_romberg_2d(Q))
+        with mp.workprec(100):
+            print("mpmath", *mp_quad_2d(mp.mpf(Qstr), shape))
+    plot_2d(Q, n=200)
+
+if __name__ == "__main__":
+    main(Qstr)

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()