Index: explore/precision.py
===================================================================
--- explore/precision.py	(revision eb2946f252d2a59f270879ef73814a9c1f221ae6)
+++ explore/precision.py	(revision eb2946f252d2a59f270879ef73814a9c1f221ae6)
@@ -0,0 +1,473 @@
+#!/usr/bin/env python
+r"""
+Show numerical precision of $2 J_1(x)/x$.
+"""
+from __future__ import division, print_function
+
+import sys
+import os
+sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), '..')))
+
+import numpy as np
+from numpy import pi, inf
+import scipy.special
+try:
+    from mpmath import mp
+except ImportError:
+    # CRUFT: mpmath split out into its own package
+    from sympy.mpmath import mp
+#import matplotlib; matplotlib.use('TkAgg')
+import pylab
+
+from sasmodels import core, data, direct_model, modelinfo
+
+class Comparator(object):
+    def __init__(self, name, mp_function, np_function, ocl_function, xaxis, limits):
+        self.name = name
+        self.mp_function = mp_function
+        self.np_function = np_function
+        self.ocl_function = ocl_function
+        self.xaxis = xaxis
+        self.limits = limits
+
+    def __repr__(self):
+        return "Comparator(%s)"%self.name
+
+    def call_mpmath(self, vec, bits=500):
+        """
+        Direct calculation using mpmath extended precision library.
+        """
+        with mp.workprec(bits):
+            return [self.mp_function(mp.mpf(x)) for x in vec]
+
+    def call_numpy(self, x, dtype):
+        """
+        Direct calculation using numpy/scipy.
+        """
+        x = np.asarray(x, dtype)
+        return self.np_function(x)
+
+    def call_ocl(self, x, dtype, platform='ocl'):
+        """
+        Calculation using sasmodels ocl libraries.
+        """
+        x = np.asarray(x, dtype)
+        model = core.build_model(self.ocl_function, dtype=dtype)
+        calculator = direct_model.DirectModel(data.empty_data1D(x), model)
+        return calculator(background=0)
+
+    def run(self, xrange="log", diff=True):
+        r"""
+        Compare accuracy of different methods for computing f.
+
+        *xrange* is log=[10^-3,10^5], linear=[1,1000], zoom[1000,1010],
+        or neg=[-100,100].
+
+        *diff* is False if showing function value rather than relative error.
+
+        *x_bits* is the precision with which the x values are specified.  The
+        default 23 should reproduce the equivalent of a single precisio
+        """
+        linear = xrange != "log"
+        if xrange == "zoom":
+            lin_min, lin_max, lin_steps = 1000, 1010, 2000
+        elif xrange == "neg":
+            lin_min, lin_max, lin_steps = -100.1, 100.1, 2000
+        else:
+            lin_min, lin_max, lin_steps = 1, 1000, 2000
+        lin_min = max(lin_min, self.limits[0])
+        lin_max = min(lin_max, self.limits[1])
+        log_min, log_max, log_steps = -3, 5, 400
+        with mp.workprec(500):
+            if linear:
+                qrf = np.linspace(lin_min, lin_max, lin_steps, dtype='single')
+                qr = [mp.mpf(float(v)) for v in qrf]
+                #qr = mp.linspace(lin_min, lin_max, lin_steps)
+            else:
+                qrf = np.logspace(log_min, log_max, log_steps, dtype='single')
+                qr = [mp.mpf(float(v)) for v in qrf]
+                #qr = [10**v for v in mp.linspace(log_min, log_max, log_steps)]
+
+        target = self.call_mpmath(qr, bits=500)
+        pylab.subplot(121)
+        self.compare(qr, 'single', target, linear, diff)
+        pylab.legend(loc='best')
+        pylab.subplot(122)
+        self.compare(qr, 'double', target, linear, diff)
+        pylab.legend(loc='best')
+        pylab.suptitle(self.name + " compared to 500-bit mpmath")
+
+    def compare(self, x, precision, target, linear=False, diff=True):
+        r"""
+        Compare the different computation methods using the given precision.
+        """
+        if precision == 'single':
+            #n=11; plotdiff(x, target, self.call_mpmath(x, n), 'mp %d bits'%n, diff=diff)
+            #n=23; plotdiff(x, target, self.call_mpmath(x, n), 'mp %d bits'%n, diff=diff)
+            pass
+        elif precision == 'double':
+            #n=53; plotdiff(x, target, self.call_mpmath(x, n), 'mp %d bits'%n, diff=diff)
+            #n=83; plotdiff(x, target, self.call_mpmath(x, n), 'mp %d bits'%n, diff=diff)
+            pass
+        plotdiff(x, target, self.call_numpy(x, precision), 'numpy '+precision, diff=diff)
+        plotdiff(x, target, self.call_ocl(x, precision, 0), 'OpenCL '+precision, diff=diff)
+        pylab.xlabel(self.xaxis)
+        if diff:
+            pylab.ylabel("relative error")
+        else:
+            pylab.ylabel(self.name)
+            pylab.semilogx(x, target, '-', label="true value")
+        if linear:
+            pylab.xscale('linear')
+
+def plotdiff(x, target, actual, label, diff=True):
+    """
+    Plot the computed value.
+
+    Use relative error if SHOW_DIFF, otherwise just plot the value directly.
+    """
+    if diff:
+        err = np.array([abs((t-a)/t) for t, a in zip(target, actual)], 'd')
+        #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 make_ocl(function, name, source=[]):
+    class Kernel(object):
+        pass
+    Kernel.__file__ = name+".py"
+    Kernel.name = name
+    Kernel.parameters = []
+    Kernel.source = source
+    Kernel.Iq = function
+    model_info = modelinfo.make_model_info(Kernel)
+    return model_info
+
+
+# =============== FUNCTION DEFINITIONS ================
+
+FUNCTIONS = {}
+def add_function(name, mp_function, np_function, ocl_function,
+                 shortname=None, xaxis="x", limits=(-inf, inf)):
+    if shortname is None:
+        shortname = name.replace('(x)', '').replace(' ', '')
+    FUNCTIONS[shortname] = Comparator(name, mp_function, np_function, ocl_function, xaxis, limits)
+
+add_function(
+    name="J0(x)",
+    mp_function=mp.j0,
+    np_function=scipy.special.j0,
+    ocl_function=make_ocl("return sas_J0(q);", "sas_J0", ["lib/polevl.c", "lib/sas_J0.c"]),
+)
+add_function(
+    name="J1(x)",
+    mp_function=mp.j1,
+    np_function=scipy.special.j1,
+    ocl_function=make_ocl("return sas_J1(q);", "sas_J1", ["lib/polevl.c", "lib/sas_J1.c"]),
+)
+add_function(
+    name="JN(-3, x)",
+    mp_function=lambda x: mp.besselj(-3, x),
+    np_function=lambda x: scipy.special.jn(-3, x),
+    ocl_function=make_ocl("return sas_JN(-3, q);", "sas_JN",
+                          ["lib/polevl.c", "lib/sas_J0.c", "lib/sas_J1.c", "lib/sas_JN.c"]),
+    shortname="J-3",
+)
+add_function(
+    name="JN(3, x)",
+    mp_function=lambda x: mp.besselj(3, x),
+    np_function=lambda x: scipy.special.jn(3, x),
+    ocl_function=make_ocl("return sas_JN(3, q);", "sas_JN",
+                          ["lib/polevl.c", "lib/sas_J0.c", "lib/sas_J1.c", "lib/sas_JN.c"]),
+    shortname="J3",
+)
+add_function(
+    name="JN(2, x)",
+    mp_function=lambda x: mp.besselj(2, x),
+    np_function=lambda x: scipy.special.jn(2, x),
+    ocl_function=make_ocl("return sas_JN(2, q);", "sas_JN",
+                          ["lib/polevl.c", "lib/sas_J0.c", "lib/sas_J1.c", "lib/sas_JN.c"]),
+    shortname="J2",
+)
+add_function(
+    name="2 J1(x)/x",
+    mp_function=lambda x: 2*mp.j1(x)/x,
+    np_function=lambda x: 2*scipy.special.j1(x)/x,
+    ocl_function=make_ocl("return sas_2J1x_x(q);", "sas_2J1x_x", ["lib/polevl.c", "lib/sas_J1.c"]),
+)
+add_function(
+    name="J1(x)",
+    mp_function=mp.j1,
+    np_function=scipy.special.j1,
+    ocl_function=make_ocl("return sas_J1(q);", "sas_J1", ["lib/polevl.c", "lib/sas_J1.c"]),
+)
+add_function(
+    name="Si(x)",
+    mp_function=mp.si,
+    np_function=lambda x: scipy.special.sici(x)[0],
+    ocl_function=make_ocl("return sas_Si(q);", "sas_Si", ["lib/sas_Si.c"]),
+)
+#import fnlib
+#add_function(
+#    name="fnlibJ1",
+#    mp_function=mp.j1,
+#    np_function=fnlib.J1,
+#    ocl_function=make_ocl("return sas_J1(q);", "sas_J1", ["lib/polevl.c", "lib/sas_J1.c"]),
+#)
+add_function(
+    name="sin(x)",
+    mp_function=mp.sin,
+    np_function=np.sin,
+    #ocl_function=make_ocl("double sn, cn; SINCOS(q,sn,cn); return sn;", "sas_sin"),
+    ocl_function=make_ocl("return sin(q);", "sas_sin"),
+)
+add_function(
+    name="sin(x)/x",
+    mp_function=lambda x: mp.sin(x)/x if x != 0 else 1,
+    ## scipy sinc function is inaccurate and has an implied pi*x term
+    #np_function=lambda x: scipy.special.sinc(x/pi),
+    ## numpy sin(x)/x needs to check for x=0
+    np_function=lambda x: np.sin(x)/x,
+    ocl_function=make_ocl("return sas_sinx_x(q);", "sas_sinc"),
+)
+add_function(
+    name="cos(x)",
+    mp_function=mp.cos,
+    np_function=np.cos,
+    #ocl_function=make_ocl("double sn, cn; SINCOS(q,sn,cn); return cn;", "sas_cos"),
+    ocl_function=make_ocl("return cos(q);", "sas_cos"),
+)
+add_function(
+    name="gamma(x)",
+    mp_function=mp.gamma,
+    np_function=scipy.special.gamma,
+    ocl_function=make_ocl("return sas_gamma(q);", "sas_gamma", ["lib/sas_gamma.c"]),
+    limits=(-3.1,10),
+)
+add_function(
+    name="erf(x)",
+    mp_function=mp.erf,
+    np_function=scipy.special.erf,
+    ocl_function=make_ocl("return sas_erf(q);", "sas_erf", ["lib/polevl.c", "lib/sas_erf.c"]),
+    limits=(-5.,5.),
+)
+add_function(
+    name="erfc(x)",
+    mp_function=mp.erfc,
+    np_function=scipy.special.erfc,
+    ocl_function=make_ocl("return sas_erfc(q);", "sas_erfc", ["lib/polevl.c", "lib/sas_erf.c"]),
+    limits=(-5.,5.),
+)
+add_function(
+    name="arctan(x)",
+    mp_function=mp.atan,
+    np_function=np.arctan,
+    ocl_function=make_ocl("return atan(q);", "sas_arctan"),
+)
+add_function(
+    name="3 j1(x)/x",
+    mp_function=lambda x: 3*(mp.sin(x)/x - mp.cos(x))/(x*x),
+    # Note: no taylor expansion near 0
+    np_function=lambda x: 3*(np.sin(x)/x - np.cos(x))/(x*x),
+    ocl_function=make_ocl("return sas_3j1x_x(q);", "sas_j1c", ["lib/sas_3j1x_x.c"]),
+)
+add_function(
+    name="fmod_2pi",
+    mp_function=lambda x: mp.fmod(x, 2*mp.pi),
+    np_function=lambda x: np.fmod(x, 2*np.pi),
+    ocl_function=make_ocl("return fmod(q, 2*M_PI);", "sas_fmod"),
+)
+
+RADIUS=3000
+LENGTH=30
+THETA=45
+def mp_cyl(x):
+    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_cyl(x):
+    f = np.float64 if x.dtype == np.float64 else np.float32
+    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)*scipy.special.j1(qr)/qr*np.sin(qh)/qh)**f(2)
+ocl_cyl = """\
+    double THETA = %(THETA).15e*M_PI_180;
+    double qr = q*%(RADIUS).15e*sin(THETA);
+    double qh = q*0.5*%(LENGTH).15e*cos(THETA);
+    return square(sas_2J1x_x(qr)*sas_sinx_x(qh));
+"""%{"LENGTH":LENGTH, "RADIUS": RADIUS, "THETA": THETA}
+add_function(
+    name="cylinder(r=%g, l=%g, theta=%g)"%(RADIUS, LENGTH, THETA),
+    mp_function=mp_cyl,
+    np_function=np_cyl,
+    ocl_function=make_ocl(ocl_cyl, "ocl_cyl", ["lib/polevl.c", "lib/sas_J1.c"]),
+    shortname="cylinder",
+    xaxis="$q/A^{-1}$",
+)
+
+lanczos_gamma = """\
+    const double coeff[] = {
+            76.18009172947146,     -86.50532032941677,
+            24.01409824083091,     -1.231739572450155,
+            0.1208650973866179e-2,-0.5395239384953e-5
+            };
+    const double x = q;
+    double tmp  = x + 5.5;
+    tmp -= (x + 0.5)*log(tmp);
+    double ser = 1.000000000190015;
+    for (int k=0; k < 6; k++) ser += coeff[k]/(x + k+1);
+    return -tmp + log(2.5066282746310005*ser/x);
+"""
+add_function(
+    name="log gamma(x)",
+    mp_function=mp.loggamma,
+    np_function=scipy.special.gammaln,
+    ocl_function=make_ocl(lanczos_gamma, "lgamma"),
+)
+
+# Alternate versions of 3 j1(x)/x, for posterity
+def taylor_3j1x_x(x):
+    """
+    Calculation using taylor series.
+    """
+    # Generate coefficients using the precision of the target value.
+    n = 5
+    cinv = [3991680, -45360, 840, -30, 3]
+    three = x.dtype.type(3)
+    p = three/np.array(cinv, x.dtype)
+    return np.polyval(p[-n:], x*x)
+add_function(
+    name="3 j1(x)/x: taylor",
+    mp_function=lambda x: 3*(mp.sin(x)/x - mp.cos(x))/(x*x),
+    np_function=taylor_3j1x_x,
+    ocl_function=make_ocl("return sas_3j1x_x(q);", "sas_j1c", ["lib/sas_3j1x_x.c"]),
+)
+def trig_3j1x_x(x):
+    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)$.
+    """
+    one = x.dtype.type(1)
+    three = x.dtype.type(3)
+    c = np.sqrt(one + x*x)
+    phi = np.arctan2(-x, one)
+    return three*(c*np.sin(x+phi))/(x*x*x)
+add_function(
+    name="3 j1(x)/x: trig",
+    mp_function=lambda x: 3*(mp.sin(x)/x - mp.cos(x))/(x*x),
+    np_function=trig_3j1x_x,
+    ocl_function=make_ocl("return sas_3j1x_x(q);", "sas_j1c", ["lib/sas_3j1x_x.c"]),
+)
+def np_2J1x_x(x):
+    """
+    numpy implementation of 2J1(x)/x using single precision algorithm
+    """
+    # pylint: disable=bad-continuation
+    f = x.dtype.type
+    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
+add_function(
+    name="2 J1(x)/x:alt",
+    mp_function=lambda x: 2*mp.j1(x)/x,
+    np_function=lambda x: np.asarray([np_2J1x_x(v) for v in x], x.dtype),
+    ocl_function=make_ocl("return sas_2J1x_x(q);", "sas_2J1x_x", ["lib/polevl.c", "lib/sas_J1.c"]),
+)
+
+ALL_FUNCTIONS = set(FUNCTIONS.keys())
+ALL_FUNCTIONS.discard("loggamma")  # OCL version not ready yet
+ALL_FUNCTIONS.discard("3j1/x:taylor")
+ALL_FUNCTIONS.discard("3j1/x:trig")
+ALL_FUNCTIONS.discard("2J1/x:alt")
+
+# =============== MAIN PROGRAM ================
+
+def usage():
+    names = ", ".join(sorted(ALL_FUNCTIONS))
+    print("""\
+usage: precision.py [-f] [--log|--linear|--zoom|--neg] name...
+where
+    -f indicates that the function value should be plotted rather than error,
+    --log indicates log stepping in [10^-3, 10^5]
+    --linear indicates linear stepping in [1, 1000]
+    --zoom indicates linear stepping in [1000, 1010]
+    --neg indicates linear stepping in [-100.1, 100.1]
+and name is "all [first]" or one of:
+    """+names)
+    sys.exit(1)
+
+def main():
+    import sys
+    diff = True
+    xrange = "log"
+    args = sys.argv[1:]
+    if '-f' in args:
+        args.remove('-f')
+        diff = False
+    for k in "log linear zoom neg".split():
+        if '--'+k in args:
+            args.remove('--'+k)
+            xrange = k
+    if not args:
+        usage()
+    if args[0] == "all":
+        cutoff = args[1] if len(args) > 1 else ""
+        args = list(sorted(ALL_FUNCTIONS))
+        args = [k for k in args if k >= cutoff]
+    if any(k not in FUNCTIONS for k in args):
+        usage()
+    multiple = len(args) > 1
+    pylab.interactive(multiple)
+    for k in args:
+        pylab.clf()
+        comparator = FUNCTIONS[k]
+        comparator.run(xrange=xrange, diff=diff)
+        if multiple:
+            raw_input()
+    if not multiple:
+        pylab.show()
+
+if __name__ == "__main__":
+    main()
Index: plore/precision/J1.py
===================================================================
--- explore/precision/J1.py	(revision 452b168ea2996d516e225564dd3b83395fb32911)
+++ 	(revision )
@@ -1,104 +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.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: plore/precision/J1c.py
===================================================================
--- explore/precision/J1c.py	(revision 452b168ea2996d516e225564dd3b83395fb32911)
+++ 	(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: plore/precision/cyl.py
===================================================================
--- explore/precision/cyl.py	(revision 452b168ea2996d516e225564dd3b83395fb32911)
+++ 	(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: plore/precision/lgamma.py
===================================================================
--- explore/precision/lgamma.py	(revision 452b168ea2996d516e225564dd3b83395fb32911)
+++ 	(revision )
@@ -1,106 +1,0 @@
-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: plore/precision/log_gamma.py
===================================================================
--- explore/precision/log_gamma.py	(revision 452b168ea2996d516e225564dd3b83395fb32911)
+++ 	(revision )
@@ -1,101 +1,0 @@
-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: plore/precision/sinc.py
===================================================================
--- explore/precision/sinc.py	(revision 452b168ea2996d516e225564dd3b83395fb32911)
+++ 	(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: plore/precision/sph_j1c.py
===================================================================
--- explore/precision/sph_j1c.py	(revision 452b168ea2996d516e225564dd3b83395fb32911)
+++ 	(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()
Index: sasmodels/models/lib/sas_3j1x_x.c
===================================================================
--- sasmodels/models/lib/sas_3j1x_x.c	(revision 473a9f131b0c8e4ae5fb862b046d8e30d4eddedb)
+++ sasmodels/models/lib/sas_3j1x_x.c	(revision eb2946f252d2a59f270879ef73814a9c1f221ae6)
@@ -46,5 +46,6 @@
 double sas_3j1x_x(double q)
 {
-    if (q < SPH_J1C_CUTOFF) {
+    // 2017-05-18 PAK - support negative q
+    if (fabs(q) < SPH_J1C_CUTOFF) {
         const double q2 = q*q;
         return (1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360.))));// + q2*(3./3991680.)))));
Index: sasmodels/models/lib/sas_J0.c
===================================================================
--- sasmodels/models/lib/sas_J0.c	(revision c8902ac8c2f740f2cd6d59d72c25f8957b7bdd8e)
+++ sasmodels/models/lib/sas_J0.c	(revision eb2946f252d2a59f270879ef73814a9c1f221ae6)
@@ -236,7 +236,8 @@
         xx = x;
 
-    if( x <= 2.0 ) {
+    // 2017-05-18 PAK - support negative x
+    if( xx <= 2.0 ) {
         z = xx * xx;
-        if( x < 1.0e-3 )
+        if( xx < 1.0e-3 )
             return( 1.0 - 0.25*z );
 
@@ -245,5 +246,5 @@
     }
 
-    q = 1.0/x;
+    q = 1.0/xx;
     w = sqrt(q);
 
Index: sasmodels/models/lib/sas_J1.c
===================================================================
--- sasmodels/models/lib/sas_J1.c	(revision 473a9f131b0c8e4ae5fb862b046d8e30d4eddedb)
+++ sasmodels/models/lib/sas_J1.c	(revision eb2946f252d2a59f270879ef73814a9c1f221ae6)
@@ -109,5 +109,5 @@
 {
 
-    double w, z, p, q, xn;
+    double w, z, p, q, xn, abs_x, sign_x;
 
     const double Z1 = 1.46819706421238932572E1;
@@ -116,22 +116,25 @@
     const double SQ2OPI = 0.79788456080286535588;
 
-    w = x;
-    if( x < 0 )
-        w = -x;
-
-    if( w <= 5.0 ) {
-        z = x * x;
+    // 2017-05-18 PAK - mathematica and mpmath use J1(-x) = -J1(x)
+    if (x < 0) {
+        abs_x = -x;
+        sign_x = -1.0;
+    } else {
+        abs_x = x;
+        sign_x = 1.0;
+    }
+
+    if( abs_x <= 5.0 ) {
+        z = abs_x * abs_x;
         w = polevl( z, RPJ1, 3 ) / p1evl( z, RQJ1, 8 );
-        w = w * x * (z - Z1) * (z - Z2);
-        return( w );
+        w = w * abs_x * (z - Z1) * (z - Z2);
+        return( sign_x * w );
     }
 
-    w = 5.0/x;
+    w = 5.0/abs_x;
     z = w * w;
-
     p = polevl( z, PPJ1, 6)/polevl( z, PQJ1, 6 );
     q = polevl( z, QPJ1, 7)/p1evl( z, QQJ1, 7 );
-
-    xn = x - THPIO4;
+    xn = abs_x - THPIO4;
 
     double sn, cn;
@@ -139,5 +142,5 @@
     p = p * cn - w * q * sn;
 
-    return( p * SQ2OPI / sqrt(x) );
+    return( sign_x * p * SQ2OPI / sqrt(abs_x) );
 }
 
@@ -179,8 +182,8 @@
     };
 
-float cephes_j1f(float x)
+float cephes_j1f(float xx)
 {
 
-    float xx, w, z, p, q, xn;
+    float x, w, z, p, q, xn;
 
     const float Z1 = 1.46819706421238932572E1;
@@ -188,12 +191,13 @@
 
 
-    xx = x;
-    if( xx < 0 )
-        xx = -x;
-
-    if( xx <= 2.0 ) {
-        z = xx * xx;
-        p = (z-Z1) * xx * polevl( z, JPJ1, 4 );
-        return( p );
+    // 2017-05-18 PAK - mathematica and mpmath use J1(-x) = -J1(x)
+    x = xx;
+    if( x < 0 )
+        x = -xx;
+
+    if( x <= 2.0 ) {
+        z = x * x;
+        p = (z-Z1) * x * polevl( z, JPJ1, 4 );
+        return( xx < 0. ? -p : p );
     }
 
@@ -204,7 +208,7 @@
     w = q*q;
     xn = q * polevl( w, PH1J1, 7) - THPIO4F;
-    p = p * cos(xn + xx);
-
-    return(p);
+    p = p * cos(xn + x);
+
+    return( xx < 0. ? -p : p );
 }
 #endif
Index: sasmodels/models/lib/sas_erf.c
===================================================================
--- sasmodels/models/lib/sas_erf.c	(revision b3796fa6bfb84621372c42111846c80bd2fa772c)
+++ sasmodels/models/lib/sas_erf.c	(revision eb2946f252d2a59f270879ef73814a9c1f221ae6)
@@ -165,8 +165,9 @@
         // the erf function instead and z < 1.0.
         //return (1.0 - cephes_erf(a));
-        z = x * x;
-        y = x * polevl(z, TD, 4) / p1evl(z, UD, 5);
-
-        return y;
+        // 2017-05-18 PAK - use erf(a) rather than erf(|a|)
+        z = a * a;
+        y = a * polevl(z, TD, 4) / p1evl(z, UD, 5);
+
+        return 1.0 - y;
     }
 
@@ -279,8 +280,9 @@
         //is explicit here for the case < 1.0
         //return (1.0 - sas_erf(a));
-        z = x * x;
-        y = x * polevl( z, TF, 6 );
-
-        return y;
+        // 2017-05-18 PAK - use erf(a) rather than erf(|a|)
+        z = a * a;
+        y = a * polevl( z, TF, 6 );
+
+        return 1.0 - y;
     }
 
