Index: explore/J1.py
===================================================================
--- explore/J1.py	(revision 95ce773c779daaa58a2e556193cd79b2e13d71d8)
+++ 	(revision )
@@ -1,106 +1,0 @@
-r"""
-Show numerical precision of $2 J_1(x)/x$.
-"""
-import sys; sys.path.insert(0, '..')
-
-import numpy as np
-try:
-    from mpmath import mp
-except:
-    # CRUFT: mpmath split out into its own package
-    from sympy.mpmath import mp
-#import matplotlib; matplotlib.use('TkAgg')
-import pylab
-
-
-SHOW_DIFF = True # True if show diff rather than function value
-#SHOW_DIFF = False # True if show diff rather than function value
-LINEAR_X = False  # True if q is linearly spaced instead of log spaced
-#LINEAR_X = True # True if q is linearly spaced instead of log spaced
-FUNCTION = "J1(x)"
-
-def mp_fn(vec, bits=500):
-    """
-    Direct calculation using sympy multiprecision library.
-    """
-    with mp.workprec(bits):
-        return [_mp_fn(mp.mpf(x)) for x in vec]
-
-def _mp_fn(x):
-    """
-    Actual function that gets evaluated.  The caller just vectorizes.
-    """
-    #return mp.mpf(2)*mp.j1(x)/x
-    return mp.j1(x)
-
-def np_fn(x, dtype):
-    """
-    Direct calculation using scipy.
-    """
-    from scipy.special import j1 as J1
-    x = np.asarray(x, dtype)
-    #return np.asarray(2, dtype)*J1(x)/x
-    return J1(x)
-
-def sasmodels_fn(x, dtype, platform='ocl'):
-    """
-    Calculation using pade approximant.
-    """
-    from sasmodels import core, data, direct_model
-    model = core.load_model('bessel', dtype=dtype)
-    calculator = direct_model.DirectModel(data.empty_data1D(x), model)
-
-    return calculator(background=0)
-
-def plotdiff(x, target, actual, label):
-    """
-    Plot the computed value.
-
-    Use relative error if SHOW_DIFF, otherwise just plot the value directly.
-    """
-    if SHOW_DIFF:
-        err = abs((target-actual)/target)
-        #err = np.clip(err, 0, 1)
-        pylab.loglog(x, err, '-', label=label)
-    else:
-        limits = np.min(target), np.max(target)
-        pylab.semilogx(x, np.clip(actual,*limits),  '-', label=label)
-
-def compare(x, precision, target):
-    r"""
-    Compare the different computation methods using the given precision.
-    """
-    #plotdiff(x, target, mp_fn(x, 11), 'mp 11 bits')
-    plotdiff(x, target, np_fn(x, precision), 'numpy '+precision)
-    plotdiff(x, target, sasmodels_fn(x, precision, 0), 'sasmodels '+precision)
-    pylab.xlabel("qr (1/Ang)")
-    if SHOW_DIFF:
-        pylab.ylabel("relative error")
-    else:
-        pylab.ylabel(FUNCTION)
-        pylab.semilogx(x, target,  '-', label="true value")
-    if LINEAR_X:
-        pylab.xscale('linear')
-
-def main():
-    r"""
-    Compare accuracy of different methods for computing $3 j_1(x)/x$.
-    :return:
-    """
-    if LINEAR_X:
-        qr = np.linspace(1,1000,2000)
-    else:
-        qr = np.logspace(-3,5,400)
-    target = np.asarray(mp_fn(qr), 'double')
-    pylab.subplot(121)
-    compare(qr, 'single', target)
-    pylab.legend(loc='best')
-    pylab.subplot(122)
-    compare(qr, 'double', target)
-    pylab.legend(loc='best')
-    pylab.suptitle(FUNCTION)
-
-if __name__ == "__main__":
-    #print "\n".join(str(x) for x in mp_J1c([1e-6,1e-5,1e-4,1e-3]))
-    main()
-    pylab.show()
Index: explore/J1c.py
===================================================================
--- explore/J1c.py	(revision 95ce773c779daaa58a2e556193cd79b2e13d71d8)
+++ 	(revision )
@@ -1,103 +1,0 @@
-r"""
-Show numerical precision of $2 J_1(x)/x$.
-"""
-import sys; sys.path.insert(0, '..')
-
-import numpy as np
-try:
-    from mpmath import mp
-except:
-    # CRUFT: mpmath split out into its own package
-    from sympy.mpmath import mp
-#import matplotlib; matplotlib.use('TkAgg')
-import pylab
-
-
-SHOW_DIFF = True # True if show diff rather than function value
-#SHOW_DIFF = False # True if show diff rather than function value
-LINEAR_X = False  # True if q is linearly spaced instead of log spaced
-#LINEAR_X = True # True if q is linearly spaced instead of log spaced
-FUNCTION = "2*J1(x)/x"
-
-def mp_fn(vec, bits=500):
-    """
-    Direct calculation using sympy multiprecision library.
-    """
-    with mp.workprec(bits):
-        return [_mp_fn(mp.mpf(x)) for x in vec]
-
-def _mp_fn(x):
-    """
-    Actual function that gets evaluated.  The caller just vectorizes.
-    """
-    return mp.mpf(2)*mp.j1(x)/x
-
-def np_fn(x, dtype):
-    """
-    Direct calculation using scipy.
-    """
-    from scipy.special import j1 as J1
-    x = np.asarray(x, dtype)
-    return np.asarray(2, dtype)*J1(x)/x
-
-def sasmodels_fn(x, dtype, platform='ocl'):
-    """
-    Calculation using pade approximant.
-    """
-    from sasmodels import core, data, direct_model
-    model = core.load_model('bessel', dtype=dtype)
-    calculator = direct_model.DirectModel(data.empty_data1D(x), model)
-    return calculator(background=0)
-
-def plotdiff(x, target, actual, label):
-    """
-    Plot the computed value.
-
-    Use relative error if SHOW_DIFF, otherwise just plot the value directly.
-    """
-    if SHOW_DIFF:
-        err = abs((target-actual)/target)
-        #err = np.clip(err, 0, 1)
-        pylab.loglog(x, err, '-', label=label)
-    else:
-        limits = np.min(target), np.max(target)
-        pylab.semilogx(x, np.clip(actual,*limits),  '-', label=label)
-
-def compare(x, precision, target):
-    r"""
-    Compare the different computation methods using the given precision.
-    """
-    #plotdiff(x, target, mp_fn(x, 11), 'mp 11 bits')
-    plotdiff(x, target, np_fn(x, precision), 'numpy '+precision)
-    plotdiff(x, target, sasmodels_fn(x, precision, 0), 'sasmodels '+precision)
-    pylab.xlabel("qr (1/Ang)")
-    if SHOW_DIFF:
-        pylab.ylabel("relative error")
-    else:
-        pylab.ylabel(FUNCTION)
-        pylab.semilogx(x, target,  '-', label="true value")
-    if LINEAR_X:
-        pylab.xscale('linear')
-
-def main():
-    r"""
-    Compare accuracy of different methods for computing $3 j_1(x)/x$.
-    :return:
-    """
-    if LINEAR_X:
-        qr = np.linspace(1,1000,2000)
-    else:
-        qr = np.logspace(-3,5,400)
-    target = np.asarray(mp_fn(qr), 'double')
-    pylab.subplot(121)
-    compare(qr, 'single', target)
-    pylab.legend(loc='best')
-    pylab.subplot(122)
-    compare(qr, 'double', target)
-    pylab.legend(loc='best')
-    pylab.suptitle(FUNCTION)
-
-if __name__ == "__main__":
-    #print "\n".join(str(x) for x in mp_J1c([1e-6,1e-5,1e-4,1e-3]))
-    main()
-    pylab.show()
Index: explore/cyl.py
===================================================================
--- explore/cyl.py	(revision 0a6da3c2133a98c9475f95364de124e647838b54)
+++ 	(revision )
@@ -1,150 +1,0 @@
-r"""
-Show numerical precision of cylinder form.
-
-Using::
-
-    qr = q r sin(t)
-    qh = q h/2 cos(t)
-    F = 2 J_1(qr)/qr sin(qh)/qh
-"""
-
-import numpy as np
-from sympy.mpmath import mp
-#import matplotlib; matplotlib.use('TkAgg')
-import pylab
-
-SHOW_DIFF = True  # True if show diff rather than function value
-LINEAR_X = False  # True if q is linearly spaced instead of log spaced
-
-RADIUS = 20
-LENGTH = 300
-CONTRAST = 5
-THETA = 45
-
-def mp_form(vec, bits=500):
-    """
-    Direct calculation using sympy multiprecision library.
-    """
-    with mp.workprec(bits):
-        return [_mp_f(mp.mpf(x)) for x in vec]
-
-def _mp_f(x):
-    """
-    Helper function for mp_j1c
-    """
-    f = mp.mpf
-    theta = f(THETA)*mp.pi/f(180)
-    qr = x * f(RADIUS)*mp.sin(theta)
-    qh = x * f(LENGTH)/f(2)*mp.cos(theta)
-    return (f(2)*mp.j1(qr)/qr * mp.sin(qh)/qh)**f(2)
-
-def np_form(x, dtype):
-    """
-    Direct calculation using scipy.
-    """
-    from scipy.special import j1
-    f = np.float64 if np.dtype(dtype) == np.float64 else np.float32
-    x = np.asarray(x, dtype)
-    theta = f(THETA)*f(np.pi)/f(180)
-    qr = x * f(RADIUS)*np.sin(theta)
-    qh = x * f(LENGTH)/f(2)*np.cos(theta)
-    return (f(2)*j1(qr)/qr*np.sin(qh)/qh)**f(2)
-
-def sasmodels_form(x, dtype):
-    f = np.float64 if np.dtype(dtype) == np.float64 else np.float32
-    x = np.asarray(x, dtype)
-    theta = f(THETA)*f(np.pi)/f(180)
-    qr = x * f(RADIUS)*np.sin(theta)
-    qh = x * f(LENGTH)/f(2)*np.cos(theta)
-    return (J1c(qr, dtype)*np.sin(qh)/qh)**f(2)
-
-def J1c(x, dtype):
-    x = np.asarray(x, dtype)
-    f = np.float64 if np.dtype(dtype) == np.float64 else np.float32
-    return np.asarray([_J1c(xi, f) for xi in x], dtype)
-
-def _J1c(x, f):
-    ax = abs(x)
-    if ax < f(8.0):
-        y = x*x
-        ans1 = f(2)*(f(72362614232.0)
-                  + y*(f(-7895059235.0)
-                  + y*(f(242396853.1)
-                  + y*(f(-2972611.439)
-                  + y*(f(15704.48260)
-                  + y*(f(-30.16036606)))))))
-        ans2 = (f(144725228442.0)
-                  + y*(f(2300535178.0)
-                  + y*(f(18583304.74)
-                  + y*(f(99447.43394)
-                  + y*(f(376.9991397)
-                  + y)))))
-        return ans1/ans2
-    else:
-        y = f(64.0)/(ax*ax)
-        xx = ax - f(2.356194491)
-        ans1 = (f(1.0)
-                  + y*(f(0.183105e-2)
-                  + y*(f(-0.3516396496e-4)
-                  + y*(f(0.2457520174e-5)
-                  + y*f(-0.240337019e-6)))))
-        ans2 = (f(0.04687499995)
-                  + y*(f(-0.2002690873e-3)
-                  + y*(f(0.8449199096e-5)
-                  + y*(f(-0.88228987e-6)
-                  + y*f(0.105787412e-6)))))
-        sn, cn = np.sin(xx), np.cos(xx)
-        ans = np.sqrt(f(0.636619772)/ax) * (cn*ans1 - (f(8.0)/ax)*sn*ans2) * f(2)/x
-        return -ans if (x < f(0.0)) else ans
-
-def plotdiff(x, target, actual, label):
-    """
-    Plot the computed value.
-
-    Use relative error if SHOW_DIFF, otherwise just plot the value directly.
-    """
-    if SHOW_DIFF:
-        err = np.clip(abs((target-actual)/target), 0, 1)
-        pylab.loglog(x, err, '-', label=label)
-    else:
-        limits = np.min(target), np.max(target)
-        pylab.semilogx(x, np.clip(actual,*limits),  '-', label=label)
-
-def compare(x, precision):
-    r"""
-    Compare the different computation methods using the given precision.
-    """
-    target = np.asarray(mp_form(x), 'double')
-    plotdiff(x, target, mp_form(x, bits=11), '11-bit')
-    plotdiff(x, target, np_form(x, precision), 'direct '+precision)
-    plotdiff(x, target, sasmodels_form(x, precision), 'sasmodels '+precision)
-    pylab.xlabel("qr (1/Ang)")
-    if SHOW_DIFF:
-        pylab.ylabel("relative error")
-    else:
-        pylab.ylabel("2 J1(x)/x")
-        pylab.semilogx(x, target,  '-', label="true value")
-    if LINEAR_X:
-        pylab.xscale('linear')
-
-def main():
-    r"""
-    Compare accuracy of different methods for computing $3 j_1(x)/x$.
-    :return:
-    """
-    if LINEAR_X:
-        qr = np.linspace(1,1000,2000)
-    else:
-        qr = np.logspace(-3,5,400)
-    pylab.subplot(121)
-    compare(qr, 'single')
-    pylab.legend(loc='best')
-    pylab.subplot(122)
-    compare(qr, 'double')
-    pylab.legend(loc='best')
-    pylab.suptitle('2 J1(x)/x')
-
-if __name__ == "__main__":
-    #print "\n".join(str(x) for x in mp_J1c([1e-6,1e-5,1e-4,1e-3]))
-    main()
-    pylab.show()
Index: explore/lgamma.py
===================================================================
--- explore/lgamma.py	(revision ad9af31d75a4f3938e9694c12b5015b4c6efec6b)
+++ 	(revision )
@@ -1,106 +1,0 @@
-r"""
-Show numerical precision of $2 J_1(x)/x$.
-"""
-import sys; sys.path.insert(0, '..')
-
-import numpy as np
-try:
-    from mpmath import mp
-except:
-    # CRUFT: mpmath split out into its own package
-    from sympy.mpmath import mp
-#import matplotlib; matplotlib.use('TkAgg')
-import pylab
-
-
-SHOW_DIFF = True # True if show diff rather than function value
-#SHOW_DIFF = False # True if show diff rather than function value
-LINEAR_X = False  # True if q is linearly spaced instead of log spaced
-#LINEAR_X = True # True if q is linearly spaced instead of log spaced
-FUNCTION = "Log_Gamma(x)"
-
-def mp_fn(vec, bits=500):
-    """
-    Direct calculation using sympy multiprecision library.
-    """
-    with mp.workprec(bits):
-        return [_mp_fn(mp.mpf(x)) for x in vec]
-
-def _mp_fn(x):
-    """
-    Actual function that gets evaluated.  The caller just vectorizes.
-    """
-    #return mp.mpf(2)*mp.j1(x)/x
-    return mp.exp(mp.loggamma(x))
-
-def np_fn(x, dtype):
-    """
-    Direct calculation using scipy.
-    """
-    from scipy.special import gammaln as logGamma
-    x = np.asarray(x, dtype)
-    #return np.asarray(2, dtype)*J1(x)/x
-    return np.exp(logGamma(x))
-
-def sasmodels_fn(x, dtype, platform='ocl'):
-    """
-    Calculation using pade approximant.
-    """
-    from sasmodels import core, data, direct_model
-    model = core.load_model('bessel', dtype=dtype)
-    calculator = direct_model.DirectModel(data.empty_data1D(x), model)
-
-    return calculator(background=0)
-
-def plotdiff(x, target, actual, label):
-    """
-    Plot the computed value.
-
-    Use relative error if SHOW_DIFF, otherwise just plot the value directly.
-    """
-    if SHOW_DIFF:
-        err = abs((target-actual)/target)
-        #err = np.clip(err, 0, 1)
-        pylab.loglog(x, err, '-', label=label)
-    else:
-        limits = np.min(target), np.max(target)
-        pylab.semilogx(x, np.clip(actual,*limits),  '-', label=label)
-
-def compare(x, precision, target):
-    r"""
-    Compare the different computation methods using the given precision.
-    """
-    #plotdiff(x, target, mp_fn(x, 11), 'mp 11 bits')
-    plotdiff(x, target, np_fn(x, precision), 'numpy '+precision)
-    plotdiff(x, target, sasmodels_fn(x, precision, 0), 'sasmodels '+precision)
-    pylab.xlabel("qr (1/Ang)")
-    if SHOW_DIFF:
-        pylab.ylabel("relative error")
-    else:
-        pylab.ylabel(FUNCTION)
-        pylab.semilogx(x, target,  '-', label="true value")
-    if LINEAR_X:
-        pylab.xscale('linear')
-
-def main():
-    r"""
-    Compare accuracy of different methods for computing $3 j_1(x)/x$.
-    :return:
-    """
-    if LINEAR_X:
-        qr = np.linspace(1,1000,2000)
-    else:
-        qr = np.logspace(-3,5,400)
-    target = np.asarray(mp_fn(qr), 'double')
-    pylab.subplot(121)
-    compare(qr, 'single', target)
-    pylab.legend(loc='best')
-    pylab.subplot(122)
-    compare(qr, 'double', target)
-    pylab.legend(loc='best')
-    pylab.suptitle(FUNCTION)
-
-if __name__ == "__main__":
-    #print "\n".join(str(x) for x in mp_J1c([1e-6,1e-5,1e-4,1e-3]))
-    main()
-    pylab.show()
Index: explore/log_gamma.py
===================================================================
--- explore/log_gamma.py	(revision e6f14104e82489f43d4acc1eeda0d50ad5408fb6)
+++ 	(revision )
@@ -1,101 +1,0 @@
-r"""
-Show numerical precision of $sin(x)/x$.
-"""
-
-import numpy as np
-import scipy.special
-from sympy.mpmath import mp
-#import matplotlib; matplotlib.use('TkAgg')
-import pylab
-
-mp.dps = 150 # number of digits to use in estimating true value
-
-SHOW_DIFF = True  # True if show diff rather than function value
-LINEAR_X = False  # True if q is linearly spaced instead of log spaced
-
-def mp_gamma(vec):
-    """
-    Direct calculation using sympy multiprecision library.
-    """
-    return [_mp_fn(mp.mpf(x)) for x in vec]
-
-def _mp_fn(x):
-    """
-    Helper funciton for mp_j1c
-    """
-    #return mp.gamma(x)
-    return mp.loggamma(x)
-
-def np_gamma(x, dtype):
-    """
-    Direct calculation using scipy.
-    """
-    x = np.asarray(x, dtype)
-    return scipy.special.gammaln(x)
-    #return scipy.special.gamma(x)
-
-def lanczos_gamma(x, dtype):
-    coeff = np.asarray((
-         76.18009172947146,     -86.50532032941677,
-         24.01409824083091,     -1.231739572450155,
-          0.1208650973866179e-2,-0.5395239384953e-5), dtype)
-
-    x = np.asarray(x, dtype)
-    tmp  = x + np.asarray(5.5, dtype)
-    tmp -= (x + np.asarray(0.5, dtype))*np.log(tmp)
-    ser  = np.ones_like(x)*np.asarray(1.000000000190015, dtype)
-    for k,c in enumerate(coeff):
-        ser += c/(x + np.asarray(k+1, dtype))
-    return -tmp+np.log(np.asarray(2.5066282746310005, dtype)*ser/x);
-
-def plotdiff(x, target, actual, label):
-    """
-    Plot the computed value.
-
-    Use relative error if SHOW_DIFF, otherwise just plot the value directly.
-    """
-    if SHOW_DIFF:
-        err = np.clip(abs((target-actual)/target), 0, 1)
-        pylab.loglog(x, err, '-', label=label)
-    else:
-        limits = np.min(target), np.max(target)
-        pylab.loglog(x, np.clip(actual,*limits),  '-', label=label)
-
-def compare(x, precision):
-    r"""
-    Compare the different computation methods using the given precision.
-    """
-    target = np.asarray(mp_gamma(x), 'double')
-    direct = np_gamma(x, precision)
-    approx = lanczos_gamma(x, precision)
-    plotdiff(x, target, direct, 'scipy '+precision)
-    plotdiff(x, target, approx, 'sasmodels '+precision)
-    pylab.xlabel("x (arbitrary)")
-    if SHOW_DIFF:
-        pylab.ylabel("relative error")
-    else:
-        pylab.ylabel("ln(gamma(x))")
-        pylab.loglog(x, target,  '-', label="true value")
-    if LINEAR_X:
-        pylab.xscale('linear')
-
-def main():
-    r"""
-    Compare accuracy of different methods for computing $3 j_1(x)/x$.
-    :return:
-    """
-    if LINEAR_X:
-        qr = np.linspace(1,1000,2000)
-    else:
-        qr = np.logspace(-3,5,400)
-    pylab.subplot(121)
-    compare(qr, 'single')
-    pylab.legend(loc='best')
-    pylab.subplot(122)
-    compare(qr, 'double')
-    pylab.legend(loc='best')
-    pylab.suptitle('ln gamma')
-
-if __name__ == "__main__":
-    main()
-    pylab.show()
Index: explore/precision/J1.py
===================================================================
--- explore/precision/J1.py	(revision 452b168ea2996d516e225564dd3b83395fb32911)
+++ explore/precision/J1.py	(revision 452b168ea2996d516e225564dd3b83395fb32911)
@@ -0,0 +1,104 @@
+r"""
+Show numerical precision of $2 J_1(x)/x$.
+"""
+import sys; sys.path.insert(0, '..')
+
+import numpy as np
+try:
+    from mpmath import mp
+except:
+    # CRUFT: mpmath split out into its own package
+    from sympy.mpmath import mp
+#import matplotlib; matplotlib.use('TkAgg')
+import pylab
+
+
+SHOW_DIFF = True # True if show diff rather than function value
+#SHOW_DIFF = False # True if show diff rather than function value
+LINEAR_X = False  # True if q is linearly spaced instead of log spaced
+#LINEAR_X = True # True if q is linearly spaced instead of log spaced
+FUNCTION = "J1(x)"
+
+def mp_fn(vec, bits=500):
+    """
+    Direct calculation using sympy multiprecision library.
+    """
+    with mp.workprec(bits):
+        return [_mp_fn(mp.mpf(x)) for x in vec]
+
+def _mp_fn(x):
+    """
+    Actual function that gets evaluated.  The caller just vectorizes.
+    """
+    return mp.j1(x)
+
+def np_fn(x, dtype):
+    """
+    Direct calculation using scipy.
+    """
+    from scipy.special import j1 as J1
+    x = np.asarray(x, dtype)
+    return J1(x)
+
+def sasmodels_fn(x, dtype, platform='ocl'):
+    """
+    Calculation using pade approximant.
+    """
+    from sasmodels import core, data, direct_model
+    model = core.load_model('bessel', dtype=dtype)
+    calculator = direct_model.DirectModel(data.empty_data1D(x), model)
+
+    return calculator(background=0)
+
+def plotdiff(x, target, actual, label):
+    """
+    Plot the computed value.
+
+    Use relative error if SHOW_DIFF, otherwise just plot the value directly.
+    """
+    if SHOW_DIFF:
+        err = abs((target-actual)/target)
+        #err = np.clip(err, 0, 1)
+        pylab.loglog(x, err, '-', label=label)
+    else:
+        limits = np.min(target), np.max(target)
+        pylab.semilogx(x, np.clip(actual,*limits),  '-', label=label)
+
+def compare(x, precision, target):
+    r"""
+    Compare the different computation methods using the given precision.
+    """
+    #plotdiff(x, target, mp_fn(x, 11), 'mp 11 bits')
+    plotdiff(x, target, np_fn(x, precision), 'numpy '+precision)
+    plotdiff(x, target, sasmodels_fn(x, precision, 0), 'sasmodels '+precision)
+    pylab.xlabel("qr (1/Ang)")
+    if SHOW_DIFF:
+        pylab.ylabel("relative error")
+    else:
+        pylab.ylabel(FUNCTION)
+        pylab.semilogx(x, target,  '-', label="true value")
+    if LINEAR_X:
+        pylab.xscale('linear')
+
+def main():
+    r"""
+    Compare accuracy of different methods for computing $3 j_1(x)/x$.
+    :return:
+    """
+    if LINEAR_X:
+        qr = np.linspace(1,1000,2000)
+    else:
+        qr = np.logspace(-3,5,400)
+    target = np.asarray(mp_fn(qr), 'double')
+    pylab.subplot(121)
+    compare(qr, 'single', target)
+    pylab.legend(loc='best')
+    pylab.subplot(122)
+    compare(qr, 'double', target)
+    pylab.legend(loc='best')
+    pylab.suptitle(FUNCTION)
+
+if __name__ == "__main__":
+    #print "\n".join(str(x) for x in mp_J1c([1e-6,1e-5,1e-4,1e-3]))
+    main()
+    pylab.show()
Index: explore/precision/J1c.py
===================================================================
--- explore/precision/J1c.py	(revision 452b168ea2996d516e225564dd3b83395fb32911)
+++ explore/precision/J1c.py	(revision 452b168ea2996d516e225564dd3b83395fb32911)
@@ -0,0 +1,103 @@
+r"""
+Show numerical precision of $2 J_1(x)/x$.
+"""
+import sys; sys.path.insert(0, '..')
+
+import numpy as np
+try:
+    from mpmath import mp
+except:
+    # CRUFT: mpmath split out into its own package
+    from sympy.mpmath import mp
+#import matplotlib; matplotlib.use('TkAgg')
+import pylab
+
+
+SHOW_DIFF = True # True if show diff rather than function value
+#SHOW_DIFF = False # True if show diff rather than function value
+LINEAR_X = False  # True if q is linearly spaced instead of log spaced
+#LINEAR_X = True # True if q is linearly spaced instead of log spaced
+FUNCTION = "2*J1(x)/x"
+
+def mp_fn(vec, bits=500):
+    """
+    Direct calculation using sympy multiprecision library.
+    """
+    with mp.workprec(bits):
+        return [_mp_fn(mp.mpf(x)) for x in vec]
+
+def _mp_fn(x):
+    """
+    Actual function that gets evaluated.  The caller just vectorizes.
+    """
+    return mp.mpf(2)*mp.j1(x)/x
+
+def np_fn(x, dtype):
+    """
+    Direct calculation using scipy.
+    """
+    from scipy.special import j1 as J1
+    x = np.asarray(x, dtype)
+    return np.asarray(2, dtype)*J1(x)/x
+
+def sasmodels_fn(x, dtype, platform='ocl'):
+    """
+    Calculation using pade approximant.
+    """
+    from sasmodels import core, data, direct_model
+    model = core.load_model('bessel', dtype=dtype)
+    calculator = direct_model.DirectModel(data.empty_data1D(x), model)
+    return calculator(background=0)
+
+def plotdiff(x, target, actual, label):
+    """
+    Plot the computed value.
+
+    Use relative error if SHOW_DIFF, otherwise just plot the value directly.
+    """
+    if SHOW_DIFF:
+        err = abs((target-actual)/target)
+        #err = np.clip(err, 0, 1)
+        pylab.loglog(x, err, '-', label=label)
+    else:
+        limits = np.min(target), np.max(target)
+        pylab.semilogx(x, np.clip(actual,*limits),  '-', label=label)
+
+def compare(x, precision, target):
+    r"""
+    Compare the different computation methods using the given precision.
+    """
+    #plotdiff(x, target, mp_fn(x, 11), 'mp 11 bits')
+    plotdiff(x, target, np_fn(x, precision), 'numpy '+precision)
+    plotdiff(x, target, sasmodels_fn(x, precision, 0), 'sasmodels '+precision)
+    pylab.xlabel("qr (1/Ang)")
+    if SHOW_DIFF:
+        pylab.ylabel("relative error")
+    else:
+        pylab.ylabel(FUNCTION)
+        pylab.semilogx(x, target,  '-', label="true value")
+    if LINEAR_X:
+        pylab.xscale('linear')
+
+def main():
+    r"""
+    Compare accuracy of different methods for computing $3 j_1(x)/x$.
+    :return:
+    """
+    if LINEAR_X:
+        qr = np.linspace(1,1000,2000)
+    else:
+        qr = np.logspace(-3,5,400)
+    target = np.asarray(mp_fn(qr), 'double')
+    pylab.subplot(121)
+    compare(qr, 'single', target)
+    pylab.legend(loc='best')
+    pylab.subplot(122)
+    compare(qr, 'double', target)
+    pylab.legend(loc='best')
+    pylab.suptitle(FUNCTION)
+
+if __name__ == "__main__":
+    #print "\n".join(str(x) for x in mp_J1c([1e-6,1e-5,1e-4,1e-3]))
+    main()
+    pylab.show()
Index: explore/precision/cyl.py
===================================================================
--- explore/precision/cyl.py	(revision 452b168ea2996d516e225564dd3b83395fb32911)
+++ explore/precision/cyl.py	(revision 452b168ea2996d516e225564dd3b83395fb32911)
@@ -0,0 +1,150 @@
+r"""
+Show numerical precision of cylinder form.
+
+Using::
+
+    qr = q r sin(t)
+    qh = q h/2 cos(t)
+    F = 2 J_1(qr)/qr sin(qh)/qh
+"""
+
+import numpy as np
+from sympy.mpmath import mp
+#import matplotlib; matplotlib.use('TkAgg')
+import pylab
+
+SHOW_DIFF = True  # True if show diff rather than function value
+LINEAR_X = False  # True if q is linearly spaced instead of log spaced
+
+RADIUS = 20
+LENGTH = 300
+CONTRAST = 5
+THETA = 45
+
+def mp_form(vec, bits=500):
+    """
+    Direct calculation using sympy multiprecision library.
+    """
+    with mp.workprec(bits):
+        return [_mp_f(mp.mpf(x)) for x in vec]
+
+def _mp_f(x):
+    """
+    Helper function for mp_j1c
+    """
+    f = mp.mpf
+    theta = f(THETA)*mp.pi/f(180)
+    qr = x * f(RADIUS)*mp.sin(theta)
+    qh = x * f(LENGTH)/f(2)*mp.cos(theta)
+    return (f(2)*mp.j1(qr)/qr * mp.sin(qh)/qh)**f(2)
+
+def np_form(x, dtype):
+    """
+    Direct calculation using scipy.
+    """
+    from scipy.special import j1
+    f = np.float64 if np.dtype(dtype) == np.float64 else np.float32
+    x = np.asarray(x, dtype)
+    theta = f(THETA)*f(np.pi)/f(180)
+    qr = x * f(RADIUS)*np.sin(theta)
+    qh = x * f(LENGTH)/f(2)*np.cos(theta)
+    return (f(2)*j1(qr)/qr*np.sin(qh)/qh)**f(2)
+
+def sasmodels_form(x, dtype):
+    f = np.float64 if np.dtype(dtype) == np.float64 else np.float32
+    x = np.asarray(x, dtype)
+    theta = f(THETA)*f(np.pi)/f(180)
+    qr = x * f(RADIUS)*np.sin(theta)
+    qh = x * f(LENGTH)/f(2)*np.cos(theta)
+    return (J1c(qr, dtype)*np.sin(qh)/qh)**f(2)
+
+def J1c(x, dtype):
+    x = np.asarray(x, dtype)
+    f = np.float64 if np.dtype(dtype) == np.float64 else np.float32
+    return np.asarray([_J1c(xi, f) for xi in x], dtype)
+
+def _J1c(x, f):
+    ax = abs(x)
+    if ax < f(8.0):
+        y = x*x
+        ans1 = f(2)*(f(72362614232.0)
+                  + y*(f(-7895059235.0)
+                  + y*(f(242396853.1)
+                  + y*(f(-2972611.439)
+                  + y*(f(15704.48260)
+                  + y*(f(-30.16036606)))))))
+        ans2 = (f(144725228442.0)
+                  + y*(f(2300535178.0)
+                  + y*(f(18583304.74)
+                  + y*(f(99447.43394)
+                  + y*(f(376.9991397)
+                  + y)))))
+        return ans1/ans2
+    else:
+        y = f(64.0)/(ax*ax)
+        xx = ax - f(2.356194491)
+        ans1 = (f(1.0)
+                  + y*(f(0.183105e-2)
+                  + y*(f(-0.3516396496e-4)
+                  + y*(f(0.2457520174e-5)
+                  + y*f(-0.240337019e-6)))))
+        ans2 = (f(0.04687499995)
+                  + y*(f(-0.2002690873e-3)
+                  + y*(f(0.8449199096e-5)
+                  + y*(f(-0.88228987e-6)
+                  + y*f(0.105787412e-6)))))
+        sn, cn = np.sin(xx), np.cos(xx)
+        ans = np.sqrt(f(0.636619772)/ax) * (cn*ans1 - (f(8.0)/ax)*sn*ans2) * f(2)/x
+        return -ans if (x < f(0.0)) else ans
+
+def plotdiff(x, target, actual, label):
+    """
+    Plot the computed value.
+
+    Use relative error if SHOW_DIFF, otherwise just plot the value directly.
+    """
+    if SHOW_DIFF:
+        err = np.clip(abs((target-actual)/target), 0, 1)
+        pylab.loglog(x, err, '-', label=label)
+    else:
+        limits = np.min(target), np.max(target)
+        pylab.semilogx(x, np.clip(actual,*limits),  '-', label=label)
+
+def compare(x, precision):
+    r"""
+    Compare the different computation methods using the given precision.
+    """
+    target = np.asarray(mp_form(x), 'double')
+    plotdiff(x, target, mp_form(x, bits=11), '11-bit')
+    plotdiff(x, target, np_form(x, precision), 'direct '+precision)
+    plotdiff(x, target, sasmodels_form(x, precision), 'sasmodels '+precision)
+    pylab.xlabel("qr (1/Ang)")
+    if SHOW_DIFF:
+        pylab.ylabel("relative error")
+    else:
+        pylab.ylabel("2 J1(x)/x")
+        pylab.semilogx(x, target,  '-', label="true value")
+    if LINEAR_X:
+        pylab.xscale('linear')
+
+def main():
+    r"""
+    Compare accuracy of different methods for computing $3 j_1(x)/x$.
+    :return:
+    """
+    if LINEAR_X:
+        qr = np.linspace(1,1000,2000)
+    else:
+        qr = np.logspace(-3,5,400)
+    pylab.subplot(121)
+    compare(qr, 'single')
+    pylab.legend(loc='best')
+    pylab.subplot(122)
+    compare(qr, 'double')
+    pylab.legend(loc='best')
+    pylab.suptitle('2 J1(x)/x')
+
+if __name__ == "__main__":
+    #print "\n".join(str(x) for x in mp_J1c([1e-6,1e-5,1e-4,1e-3]))
+    main()
+    pylab.show()
Index: explore/precision/lgamma.py
===================================================================
--- explore/precision/lgamma.py	(revision 452b168ea2996d516e225564dd3b83395fb32911)
+++ explore/precision/lgamma.py	(revision 452b168ea2996d516e225564dd3b83395fb32911)
@@ -0,0 +1,106 @@
+r"""
+Show numerical precision of $\ln \Gamma(x)$.
+"""
+import sys; sys.path.insert(0, '..')
+
+import numpy as np
+try:
+    from mpmath import mp
+except:
+    # CRUFT: mpmath split out into its own package
+    from sympy.mpmath import mp
+#import matplotlib; matplotlib.use('TkAgg')
+import pylab
+
+
+SHOW_DIFF = True # True if show diff rather than function value
+#SHOW_DIFF = False # True if show diff rather than function value
+LINEAR_X = False  # True if q is linearly spaced instead of log spaced
+#LINEAR_X = True # True if q is linearly spaced instead of log spaced
+FUNCTION = "Log_Gamma(x)"
+
+def mp_fn(vec, bits=500):
+    """
+    Direct calculation using sympy multiprecision library.
+    """
+    with mp.workprec(bits):
+        return [_mp_fn(mp.mpf(x)) for x in vec]
+
+def _mp_fn(x):
+    """
+    Actual function that gets evaluated.  The caller just vectorizes.
+    """
+    #return mp.mpf(2)*mp.j1(x)/x
+    return mp.exp(mp.loggamma(x))
+
+def np_fn(x, dtype):
+    """
+    Direct calculation using scipy.
+    """
+    from scipy.special import gammaln as logGamma
+    x = np.asarray(x, dtype)
+    #return np.asarray(2, dtype)*J1(x)/x
+    return np.exp(logGamma(x))
+
+def sasmodels_fn(x, dtype, platform='ocl'):
+    """
+    Calculation using pade approximant.
+    """
+    from sasmodels import core, data, direct_model
+    model = core.load_model('bessel', dtype=dtype)
+    calculator = direct_model.DirectModel(data.empty_data1D(x), model)
+
+    return calculator(background=0)
+
+def plotdiff(x, target, actual, label):
+    """
+    Plot the computed value.
+
+    Use relative error if SHOW_DIFF, otherwise just plot the value directly.
+    """
+    if SHOW_DIFF:
+        err = abs((target-actual)/target)
+        #err = np.clip(err, 0, 1)
+        pylab.loglog(x, err, '-', label=label)
+    else:
+        limits = np.min(target), np.max(target)
+        pylab.semilogx(x, np.clip(actual,*limits),  '-', label=label)
+
+def compare(x, precision, target):
+    r"""
+    Compare the different computation methods using the given precision.
+    """
+    #plotdiff(x, target, mp_fn(x, 11), 'mp 11 bits')
+    plotdiff(x, target, np_fn(x, precision), 'numpy '+precision)
+    plotdiff(x, target, sasmodels_fn(x, precision, 0), 'sasmodels '+precision)
+    pylab.xlabel("qr (1/Ang)")
+    if SHOW_DIFF:
+        pylab.ylabel("relative error")
+    else:
+        pylab.ylabel(FUNCTION)
+        pylab.semilogx(x, target,  '-', label="true value")
+    if LINEAR_X:
+        pylab.xscale('linear')
+
+def main():
+    r"""
+    Compare accuracy of different methods for computing $3 j_1(x)/x$.
+    :return:
+    """
+    if LINEAR_X:
+        qr = np.linspace(1,1000,2000)
+    else:
+        qr = np.logspace(-3,5,400)
+    target = np.asarray(mp_fn(qr), 'double')
+    pylab.subplot(121)
+    compare(qr, 'single', target)
+    pylab.legend(loc='best')
+    pylab.subplot(122)
+    compare(qr, 'double', target)
+    pylab.legend(loc='best')
+    pylab.suptitle(FUNCTION)
+
+if __name__ == "__main__":
+    #print "\n".join(str(x) for x in mp_J1c([1e-6,1e-5,1e-4,1e-3]))
+    main()
+    pylab.show()
Index: explore/precision/log_gamma.py
===================================================================
--- explore/precision/log_gamma.py	(revision 452b168ea2996d516e225564dd3b83395fb32911)
+++ explore/precision/log_gamma.py	(revision 452b168ea2996d516e225564dd3b83395fb32911)
@@ -0,0 +1,101 @@
+r"""
+Show numerical precision of $\ln \Gamma(x)$.
+"""
+
+import numpy as np
+import scipy.special
+from sympy.mpmath import mp
+#import matplotlib; matplotlib.use('TkAgg')
+import pylab
+
+mp.dps = 150 # number of digits to use in estimating true value
+
+SHOW_DIFF = True  # True if show diff rather than function value
+LINEAR_X = False  # True if q is linearly spaced instead of log spaced
+
+def mp_gamma(vec):
+    """
+    Direct calculation using sympy multiprecision library.
+    """
+    return [_mp_fn(mp.mpf(x)) for x in vec]
+
+def _mp_fn(x):
+    """
+    Helper funciton for mp_j1c
+    """
+    #return mp.gamma(x)
+    return mp.loggamma(x)
+
+def np_gamma(x, dtype):
+    """
+    Direct calculation using scipy.
+    """
+    x = np.asarray(x, dtype)
+    return scipy.special.gammaln(x)
+    #return scipy.special.gamma(x)
+
+def lanczos_gamma(x, dtype):
+    coeff = np.asarray((
+         76.18009172947146,     -86.50532032941677,
+         24.01409824083091,     -1.231739572450155,
+          0.1208650973866179e-2,-0.5395239384953e-5), dtype)
+
+    x = np.asarray(x, dtype)
+    tmp  = x + np.asarray(5.5, dtype)
+    tmp -= (x + np.asarray(0.5, dtype))*np.log(tmp)
+    ser  = np.ones_like(x)*np.asarray(1.000000000190015, dtype)
+    for k,c in enumerate(coeff):
+        ser += c/(x + np.asarray(k+1, dtype))
+    return -tmp+np.log(np.asarray(2.5066282746310005, dtype)*ser/x);
+
+def plotdiff(x, target, actual, label):
+    """
+    Plot the computed value.
+
+    Use relative error if SHOW_DIFF, otherwise just plot the value directly.
+    """
+    if SHOW_DIFF:
+        err = np.clip(abs((target-actual)/target), 0, 1)
+        pylab.loglog(x, err, '-', label=label)
+    else:
+        limits = np.min(target), np.max(target)
+        pylab.loglog(x, np.clip(actual,*limits),  '-', label=label)
+
+def compare(x, precision):
+    r"""
+    Compare the different computation methods using the given precision.
+    """
+    target = np.asarray(mp_gamma(x), 'double')
+    direct = np_gamma(x, precision)
+    approx = lanczos_gamma(x, precision)
+    plotdiff(x, target, direct, 'scipy '+precision)
+    plotdiff(x, target, approx, 'sasmodels '+precision)
+    pylab.xlabel("x (arbitrary)")
+    if SHOW_DIFF:
+        pylab.ylabel("relative error")
+    else:
+        pylab.ylabel("ln(gamma(x))")
+        pylab.loglog(x, target,  '-', label="true value")
+    if LINEAR_X:
+        pylab.xscale('linear')
+
+def main():
+    r"""
+    Compare accuracy of different methods for computing $3 j_1(x)/x$.
+    :return:
+    """
+    if LINEAR_X:
+        qr = np.linspace(1,1000,2000)
+    else:
+        qr = np.logspace(-3,5,400)
+    pylab.subplot(121)
+    compare(qr, 'single')
+    pylab.legend(loc='best')
+    pylab.subplot(122)
+    compare(qr, 'double')
+    pylab.legend(loc='best')
+    pylab.suptitle('ln gamma')
+
+if __name__ == "__main__":
+    main()
+    pylab.show()
Index: explore/precision/sinc.py
===================================================================
--- explore/precision/sinc.py	(revision 452b168ea2996d516e225564dd3b83395fb32911)
+++ explore/precision/sinc.py	(revision 452b168ea2996d516e225564dd3b83395fb32911)
@@ -0,0 +1,84 @@
+r"""
+Show numerical precision of $sin(x)/x$.
+"""
+
+import numpy as np
+from sympy.mpmath import mp
+#import matplotlib; matplotlib.use('TkAgg')
+import pylab
+
+mp.dps = 150 # number of digits to use in estimating true value
+
+SHOW_DIFF = True  # True if show diff rather than function value
+LINEAR_X = False  # True if q is linearly spaced instead of log spaced
+
+def mp_sinc(vec):
+    """
+    Direct calculation using sympy multiprecision library.
+    """
+    return [_mp_sinc(mp.mpf(x)) for x in vec]
+
+def _mp_sinc(x):
+    """
+    Helper funciton for mp_j1c
+    """
+    return mp.sin(x)/x
+
+def np_sinc(x, dtype):
+    """
+    Direct calculation using scipy.
+    """
+    x = np.asarray(x, dtype)
+    return np.sin(x)/x
+    #return np.asarray(np.sin(np.double(x))/np.double(x),dtype)
+
+def plotdiff(x, target, actual, label):
+    """
+    Plot the computed value.
+
+    Use relative error if SHOW_DIFF, otherwise just plot the value directly.
+    """
+    if SHOW_DIFF:
+        err = np.clip(abs((target-actual)/target), 0, 1)
+        pylab.loglog(x, err, '-', label=label)
+    else:
+        limits = np.min(target), np.max(target)
+        pylab.semilogx(x, np.clip(actual,*limits),  '-', label=label)
+
+def compare(x, precision):
+    r"""
+    Compare the different computation methods using the given precision.
+    """
+    target = np.asarray(mp_sinc(x), 'double')
+    direct = np_sinc(x, precision)
+    plotdiff(x, target, direct, 'direct '+precision)
+    pylab.xlabel("qr (1/Ang)")
+    if SHOW_DIFF:
+        pylab.ylabel("relative error")
+    else:
+        pylab.ylabel("sin(x)/x")
+        pylab.semilogx(x, target,  '-', label="true value")
+    if LINEAR_X:
+        pylab.xscale('linear')
+
+def main():
+    r"""
+    Compare accuracy of different methods for computing $3 j_1(x)/x$.
+    :return:
+    """
+    if LINEAR_X:
+        qr = np.linspace(1,1000,2000)
+    else:
+        qr = np.logspace(-3,5,400)
+    pylab.subplot(121)
+    compare(qr, 'single')
+    pylab.legend(loc='best')
+    pylab.subplot(122)
+    compare(qr, 'double')
+    pylab.legend(loc='best')
+    pylab.suptitle('sin(x)/x')
+
+if __name__ == "__main__":
+    print "\n".join(str(x) for x in mp_sinc([1e-6,1e-5,1e-4,1e-3]))
+    main()
+    pylab.show()
Index: explore/precision/sph_j1c.py
===================================================================
--- explore/precision/sph_j1c.py	(revision 452b168ea2996d516e225564dd3b83395fb32911)
+++ explore/precision/sph_j1c.py	(revision 452b168ea2996d516e225564dd3b83395fb32911)
@@ -0,0 +1,145 @@
+r"""
+Show numerical precision of $3 j_1(x)/x$ vs. its Taylor series expansion.
+
+The choice of the number of terms in the series and the cutoff value for
+switching between series and direct calculation depends on the numeric
+precision.
+
+Point where direct calculation reaches machine precision::
+
+    single machine precision eps 3e-8 at qr=1.1 (see below)
+    double machine precision eps 4e-16 at qr=1.1
+
+Point where Taylor series reaches machine precision (eps), where taylor
+series matches direct calculation (cross) and the error at that point::
+
+    prec   n eps  cross  error
+    single 3 0.28  0.4   6.2e-7
+    single 4 0.68  0.7   2.3e-7
+    single 5 1.18  1.2   7.5e-8
+    double 3 0.01  0.03  2.3e-13
+    double 4 0.06  0.1   3.1e-14
+    double 5 0.16  0.2   5.0e-15
+
+Note: relative error on single precision starts increase on the direct
+method at qr=1.1, rising from 3e-8 to 5e-5 by qr=1e3.  This should be
+safe for the sans range, with objects of 100 nm supported to a q of 0.1
+while maintaining 5 digits of precision.  For usans/sesans, the objects
+are larger but the q is smaller, so again it should be fine.
+"""
+
+import numpy as np
+from sympy.mpmath import mp
+#import matplotlib; matplotlib.use('TkAgg')
+import pylab
+
+mp.dps = 150 # number of digits to use in estimating true value
+
+SHOW_DIFF = True  # True if show diff rather than function value
+LINEAR_X = False  # True if q is linearly spaced instead of log spaced
+
+def mp_j1c(vec):
+    """
+    Direct calculation using sympy multiprecision library.
+    """
+    return [_mp_j1c(mp.mpf(x)) for x in vec]
+
+def _mp_j1c(x):
+    """
+    Helper funciton for mp_j1c
+    """
+    return mp.mpf(3)*(mp.sin(x)/x - mp.cos(x))/(x*x)
+
+def np_j1c(x, dtype):
+    """
+    Direct calculation using numpy.
+    """
+    x = np.asarray(x, dtype)
+    return np.asarray(3, dtype)*(np.sin(x) - x*np.cos(x))/(x*x*x)
+
+def lin_j1c(x, dtype):
+    r"""
+    Direct calculation using linear combination of sin/cos.
+
+    Use the following trig identity:
+
+    .. math::
+
+        a \sin(x) + b \cos(x) = c \sin(x + \phi)
+
+    where $c = \surd(a^2+b^2)$ and $\phi = \tan^{-1}(b/a) to calculate the
+    numerator $\sin(x) - x\cos(x)$.
+    """
+    x = np.asarray(x, dtype)
+    c = np.sqrt(np.asarray(1,dtype) + x*x)
+    phi = np.arctan2(-np.asarray(x,dtype),np.asarray(1,dtype))
+    return np.asarray(3, dtype)*(c*np.sin(x+phi))/(x*x*x)
+
+def taylor_j1c(x, dtype, n):
+    """
+    Calculation using taylor series.
+    """
+    # Generate coefficients using the precision of the target value.
+    cinv = [3991680, -45360, 840, -30, 3]
+    x = np.asarray(x, dtype)
+    p = np.asarray(3, dtype)/np.array(cinv, dtype)
+    return np.polyval(p[-n:], x*x)
+
+def plotdiff(x, target, actual, label):
+    """
+    Plot the computed value.
+
+    Use relative error if SHOW_DIFF, otherwise just plot the value directly.
+    """
+    if SHOW_DIFF:
+        err = np.clip(abs((target-actual)/target), 0, 1)
+        pylab.loglog(x, err, '-', label=label)
+    else:
+        limits = np.min(target), np.max(target)
+        pylab.semilogx(x, np.clip(actual,*limits),  '-', label=label)
+
+def compare(x, precision):
+    r"""
+    Compare the different computation methods using the given precision.
+    """
+    target = np.asarray(mp_j1c(x), 'double')
+    direct = np_j1c(x, precision)
+    comb = lin_j1c(x, precision)
+    taylor3 = taylor_j1c(x, precision, 3)
+    taylor4 = taylor_j1c(x, precision, 4)
+    taylor5 = taylor_j1c(x, precision, 5)
+    plotdiff(x, target, direct, 'direct '+precision)
+    #plotdiff(x, target, comb, 'c sin(x+phi) '+precision)
+    plotdiff(x, target, taylor3, 'taylor 3 '+precision)
+    plotdiff(x, target, taylor4, 'taylor 4 '+precision)
+    plotdiff(x, target, taylor5, 'taylor 5 '+precision)
+    pylab.xlabel("qr (1/Ang)")
+    if SHOW_DIFF:
+        pylab.ylabel("relative error")
+    else:
+        pylab.ylabel("3 j1(x)/x")
+        pylab.semilogx(x, target,  '-', label="true value")
+    if LINEAR_X:
+        pylab.xscale('linear')
+
+def main():
+    r"""
+    Compare accuracy of different methods for computing $3 j_1(x)/x$.
+    :return:
+    """
+    if LINEAR_X:
+        qr = np.linspace(1,1000,2000)
+    else:
+        qr = np.logspace(-3,5,400)
+    pylab.subplot(121)
+    compare(qr, 'single')
+    pylab.legend(loc='best')
+    pylab.subplot(122)
+    compare(qr, 'double')
+    pylab.legend(loc='best')
+    pylab.suptitle('3 j1(x)/x')
+
+if __name__ == "__main__":
+    #print "\n".join(str(x) for x in mp_j1c([1e-6,1e-5,1e-4,1e-3]))
+    main()
+    pylab.show()
Index: explore/sinc.py
===================================================================
--- explore/sinc.py	(revision e6f14104e82489f43d4acc1eeda0d50ad5408fb6)
+++ 	(revision )
@@ -1,84 +1,0 @@
-r"""
-Show numerical precision of $sin(x)/x$.
-"""
-
-import numpy as np
-from sympy.mpmath import mp
-#import matplotlib; matplotlib.use('TkAgg')
-import pylab
-
-mp.dps = 150 # number of digits to use in estimating true value
-
-SHOW_DIFF = True  # True if show diff rather than function value
-LINEAR_X = False  # True if q is linearly spaced instead of log spaced
-
-def mp_sinc(vec):
-    """
-    Direct calculation using sympy multiprecision library.
-    """
-    return [_mp_sinc(mp.mpf(x)) for x in vec]
-
-def _mp_sinc(x):
-    """
-    Helper funciton for mp_j1c
-    """
-    return mp.sin(x)/x
-
-def np_sinc(x, dtype):
-    """
-    Direct calculation using scipy.
-    """
-    x = np.asarray(x, dtype)
-    return np.sin(x)/x
-    #return np.asarray(np.sin(np.double(x))/np.double(x),dtype)
-
-def plotdiff(x, target, actual, label):
-    """
-    Plot the computed value.
-
-    Use relative error if SHOW_DIFF, otherwise just plot the value directly.
-    """
-    if SHOW_DIFF:
-        err = np.clip(abs((target-actual)/target), 0, 1)
-        pylab.loglog(x, err, '-', label=label)
-    else:
-        limits = np.min(target), np.max(target)
-        pylab.semilogx(x, np.clip(actual,*limits),  '-', label=label)
-
-def compare(x, precision):
-    r"""
-    Compare the different computation methods using the given precision.
-    """
-    target = np.asarray(mp_sinc(x), 'double')
-    direct = np_sinc(x, precision)
-    plotdiff(x, target, direct, 'direct '+precision)
-    pylab.xlabel("qr (1/Ang)")
-    if SHOW_DIFF:
-        pylab.ylabel("relative error")
-    else:
-        pylab.ylabel("sin(x)/x")
-        pylab.semilogx(x, target,  '-', label="true value")
-    if LINEAR_X:
-        pylab.xscale('linear')
-
-def main():
-    r"""
-    Compare accuracy of different methods for computing $3 j_1(x)/x$.
-    :return:
-    """
-    if LINEAR_X:
-        qr = np.linspace(1,1000,2000)
-    else:
-        qr = np.logspace(-3,5,400)
-    pylab.subplot(121)
-    compare(qr, 'single')
-    pylab.legend(loc='best')
-    pylab.subplot(122)
-    compare(qr, 'double')
-    pylab.legend(loc='best')
-    pylab.suptitle('sin(x)/x')
-
-if __name__ == "__main__":
-    print "\n".join(str(x) for x in mp_sinc([1e-6,1e-5,1e-4,1e-3]))
-    main()
-    pylab.show()
Index: explore/sph_j1c.py
===================================================================
--- explore/sph_j1c.py	(revision e6f14104e82489f43d4acc1eeda0d50ad5408fb6)
+++ 	(revision )
@@ -1,145 +1,0 @@
-r"""
-Show numerical precision of $3 j_1(x)/x$ vs. its Taylor series expansion.
-
-The choice of the number of terms in the series and the cutoff value for
-switching between series and direct calculation depends on the numeric
-precision.
-
-Point where direct calculation reaches machine precision::
-
-    single machine precision eps 3e-8 at qr=1.1 (see below)
-    double machine precision eps 4e-16 at qr=1.1
-
-Point where Taylor series reaches machine precision (eps), where taylor
-series matches direct calculation (cross) and the error at that point::
-
-    prec   n eps  cross  error
-    single 3 0.28  0.4   6.2e-7
-    single 4 0.68  0.7   2.3e-7
-    single 5 1.18  1.2   7.5e-8
-    double 3 0.01  0.03  2.3e-13
-    double 4 0.06  0.1   3.1e-14
-    double 5 0.16  0.2   5.0e-15
-
-Note: relative error on single precision starts increase on the direct
-method at qr=1.1, rising from 3e-8 to 5e-5 by qr=1e3.  This should be
-safe for the sans range, with objects of 100 nm supported to a q of 0.1
-while maintaining 5 digits of precision.  For usans/sesans, the objects
-are larger but the q is smaller, so again it should be fine.
-"""
-
-import numpy as np
-from sympy.mpmath import mp
-#import matplotlib; matplotlib.use('TkAgg')
-import pylab
-
-mp.dps = 150 # number of digits to use in estimating true value
-
-SHOW_DIFF = True  # True if show diff rather than function value
-LINEAR_X = False  # True if q is linearly spaced instead of log spaced
-
-def mp_j1c(vec):
-    """
-    Direct calculation using sympy multiprecision library.
-    """
-    return [_mp_j1c(mp.mpf(x)) for x in vec]
-
-def _mp_j1c(x):
-    """
-    Helper funciton for mp_j1c
-    """
-    return mp.mpf(3)*(mp.sin(x)/x - mp.cos(x))/(x*x)
-
-def np_j1c(x, dtype):
-    """
-    Direct calculation using numpy.
-    """
-    x = np.asarray(x, dtype)
-    return np.asarray(3, dtype)*(np.sin(x) - x*np.cos(x))/(x*x*x)
-
-def lin_j1c(x, dtype):
-    r"""
-    Direct calculation using linear combination of sin/cos.
-
-    Use the following trig identity:
-
-    .. math::
-
-        a \sin(x) + b \cos(x) = c \sin(x + \phi)
-
-    where $c = \surd(a^2+b^2)$ and $\phi = \tan^{-1}(b/a) to calculate the
-    numerator $\sin(x) - x\cos(x)$.
-    """
-    x = np.asarray(x, dtype)
-    c = np.sqrt(np.asarray(1,dtype) + x*x)
-    phi = np.arctan2(-np.asarray(x,dtype),np.asarray(1,dtype))
-    return np.asarray(3, dtype)*(c*np.sin(x+phi))/(x*x*x)
-
-def taylor_j1c(x, dtype, n):
-    """
-    Calculation using taylor series.
-    """
-    # Generate coefficients using the precision of the target value.
-    cinv = [3991680, -45360, 840, -30, 3]
-    x = np.asarray(x, dtype)
-    p = np.asarray(3, dtype)/np.array(cinv, dtype)
-    return np.polyval(p[-n:], x*x)
-
-def plotdiff(x, target, actual, label):
-    """
-    Plot the computed value.
-
-    Use relative error if SHOW_DIFF, otherwise just plot the value directly.
-    """
-    if SHOW_DIFF:
-        err = np.clip(abs((target-actual)/target), 0, 1)
-        pylab.loglog(x, err, '-', label=label)
-    else:
-        limits = np.min(target), np.max(target)
-        pylab.semilogx(x, np.clip(actual,*limits),  '-', label=label)
-
-def compare(x, precision):
-    r"""
-    Compare the different computation methods using the given precision.
-    """
-    target = np.asarray(mp_j1c(x), 'double')
-    direct = np_j1c(x, precision)
-    comb = lin_j1c(x, precision)
-    taylor3 = taylor_j1c(x, precision, 3)
-    taylor4 = taylor_j1c(x, precision, 4)
-    taylor5 = taylor_j1c(x, precision, 5)
-    plotdiff(x, target, direct, 'direct '+precision)
-    #plotdiff(x, target, comb, 'c sin(x+phi) '+precision)
-    plotdiff(x, target, taylor3, 'taylor 3 '+precision)
-    plotdiff(x, target, taylor4, 'taylor 4 '+precision)
-    plotdiff(x, target, taylor5, 'taylor 5 '+precision)
-    pylab.xlabel("qr (1/Ang)")
-    if SHOW_DIFF:
-        pylab.ylabel("relative error")
-    else:
-        pylab.ylabel("3 j1(x)/x")
-        pylab.semilogx(x, target,  '-', label="true value")
-    if LINEAR_X:
-        pylab.xscale('linear')
-
-def main():
-    r"""
-    Compare accuracy of different methods for computing $3 j_1(x)/x$.
-    :return:
-    """
-    if LINEAR_X:
-        qr = np.linspace(1,1000,2000)
-    else:
-        qr = np.logspace(-3,5,400)
-    pylab.subplot(121)
-    compare(qr, 'single')
-    pylab.legend(loc='best')
-    pylab.subplot(122)
-    compare(qr, 'double')
-    pylab.legend(loc='best')
-    pylab.suptitle('3 j1(x)/x')
-
-if __name__ == "__main__":
-    #print "\n".join(str(x) for x in mp_j1c([1e-6,1e-5,1e-4,1e-3]))
-    main()
-    pylab.show()
