-                      r31eea1f
+                      r20fe0cd
+#!/usr/bin/env python
 """
 Asymmetric shape integration
 …
 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
 …
 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
 …
 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
 …
 """
 def mp_quad(q, shape):
+# 2D integration functions
+def mp_quad_2d(q, shape):
     evals = [0]
     def integrand(theta, phi):
 …
     return evals[0], Iq
 def kernel(q, theta, phi):
+def kernel_2d(q, theta, phi):
     """
     S(q) kernel for paracrystal forms.
 …
     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
 …
     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]
 …
     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,))
 …
 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.
 …
     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]
 …
     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
 …
     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
 …
     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,
 …
     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]
 …
     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
 …
     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)")
 …
     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)

Note: See TracChangeset for help on using the changeset viewer.

SasView

Changeset 20fe0cd in sasmodels for explore/asymint.py

Legend:

explore/asymint.py

Download in other formats: