Index: explore/J1c.py
===================================================================
--- explore/J1c.py	(revision e6f14104e82489f43d4acc1eeda0d50ad5408fb6)
+++ explore/J1c.py	(revision e6f14104e82489f43d4acc1eeda0d50ad5408fb6)
@@ -0,0 +1,121 @@
+r"""
+Show numerical precision of $2 J_1(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_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(2)*mp.j1(x)/x
+
+def np_j1c(x, dtype):
+    """
+    Direct calculation using scipy.
+    """
+    from scipy.special import j1
+    x = np.asarray(x, dtype)
+    return np.asarray(2, dtype)*j1(x)/x
+
+def cephes_j1c(x, dtype, n):
+    """
+    Calculation using pade approximant.
+    """
+    x = np.asarray(x, dtype)
+    ans = np.empty_like(x)
+    ax = abs(x)
+
+    # Branch a
+    a_idx = ax < 8.0
+    a_xsq = x[a_idx]**2
+    a_coeff1 = list(reversed((72362614232.0, -7895059235.0, 242396853.1, -2972611.439, 15704.48260, -30.16036606)))
+    a_coeff2 = list(reversed((144725228442.0, 2300535178.0, 18583304.74, 99447.43394, 376.9991397, 1.0)))
+    a_ans1 = np.polyval(a_coeff1[n:], a_xsq)
+    a_ans2 = np.polyval(a_coeff2[n:], a_xsq)
+    ans[a_idx] = 2*a_ans1/a_ans2
+
+    # Branch b
+    b_idx = ~a_idx
+    b_ax = ax[b_idx]
+    b_x = x[b_idx]
+
+    b_y = 64.0/(b_ax**2)
+    b_xx = b_ax - 2.356194491
+    b_coeff1 = list(reversed((1.0, 0.183105e-2, -0.3516396496e-4, 0.2457520174e-5, -0.240337019e-6)))
+    b_coeff2 = list(reversed((0.04687499995, -0.2002690873e-3, 0.8449199096e-5, -0.88228987e-6, 0.105787412e-6)))
+    b_ans1 = np.polyval(b_coeff1[n:], b_y)
+    b_ans2 = np.polyval(b_coeff2[n:], b_y)
+    b_sn, b_cn = np.sin(b_xx), np.cos(b_xx)
+    ans[b_idx] = np.sign(b_x)*np.sqrt(0.636619772/b_ax) * (b_cn*b_ans1 - (8.0/b_ax)*b_sn*b_ans2)*2.0/b_x
+
+    return 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_J1c(x), 'double')
+    direct = np_j1c(x, precision)
+    approx0 = cephes_j1c(x, precision, 0)
+    approx1 = cephes_j1c(x, precision, 1)
+    plotdiff(x, target, direct, 'direct '+precision)
+    plotdiff(x, target, approx0, 'cephes '+precision)
+    #plotdiff(x, target, approx1, 'reduced '+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/log_gamma.py
===================================================================
--- explore/log_gamma.py	(revision e6f14104e82489f43d4acc1eeda0d50ad5408fb6)
+++ explore/log_gamma.py	(revision e6f14104e82489f43d4acc1eeda0d50ad5408fb6)
@@ -0,0 +1,101 @@
+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/rpa.py
===================================================================
--- explore/rpa.py	(revision e6f14104e82489f43d4acc1eeda0d50ad5408fb6)
+++ explore/rpa.py	(revision e6f14104e82489f43d4acc1eeda0d50ad5408fb6)
@@ -0,0 +1,202 @@
+from collections import namedtuple
+import numpy as np
+from numpy import sqrt, exp, expm1
+
+AVOGADRO = 6.022e23
+
+Polymer = namedtuple("Polymer", "n phi v a b".split())
+def E(q, poly):
+    qrsq = (q*Rg(poly))**2
+    retval = exp(-qrsq)
+    return retval
+
+def F(q, poly):
+    qrsq = (q*Rg(poly))**2
+    retval = -expm1(-qrsq)/qrsq
+    return retval
+
+def P_ii(q, poly):
+    qrsq = (q*Rg(poly))**2
+    retval = 2 * (expm1(-qrsq) + qrsq)/qrsq**2
+    return retval
+
+def P_ij(q, poly_list):
+    i, j = poly_list[0], poly_list[-1]
+    retval = F(q, i) * np.prod([E(q,p) for p in poly_list[1:-1]]) * F(q, j)
+    return retval
+
+def Rg(poly):
+    return sqrt(poly.n/6.)*poly.a
+
+def S0_ii(q, poly):
+    retval = poly.n*poly.phi*poly.v*P_ii(q, poly)
+    return retval
+
+def S0_ij(q, poly_list):
+    i,j = poly_list[0], poly_list[-1]
+    retval = sqrt(i.n*i.phi*i.v*j.n*j.phi*j.v) * P_ij(q, poly_list)
+    return retval
+
+def drho(poly, base):
+    return (poly.b/poly.v -  base.b/base.v)*1e-13*sqrt(AVOGADRO)
+
+def binary_blend(q, C, D, Kcd):
+    """
+    de Gennes, 1979
+    """
+    S0cc = S0_ii(q, C)
+    S0dd = S0_ii(q, D)
+    vcc = 1/S0dd - 2*Kcd #/v0
+    Scc = S0cc/(1 + vcc*S0cc)
+    rhocd = drho(C,D)
+    Iq = rhocd**2 * Scc
+    return Iq
+
+def ternary_blend(q, B, C, D, Kbc, Kbd, Kcd):
+    S0bb = S0_ii(q, B)
+    S0cc = S0_ii(q, C)
+    S0dd = S0_ii(q, D)
+    vbb = 1/S0dd - 2*Kbd
+    vcc = 1/S0dd - 2*Kcd
+    vbc = 1/S0dd + Kbd - Kbc - Kcd
+    rhobd = drho(B,D)
+    rhocd = drho(C,D)
+    det = (1+vbb*S0bb)*(1+vcc*S0cc) - vbc**2*S0bb*S0cc
+    Sbb = S0bb*(1+vcc*S0cc)/det
+    Scc = S0cc*(1+vbb*S0bb)/det
+    Sbc = -S0bb*vbc*S0cc/det
+    Iq = rhobd**2*Sbb + rhocd**2*Scc + 2*rhobd*rhocd*Sbc
+    return Iq
+
+def diblock_copolymer(q, C, D, Kcd):
+    """
+    Leibler, 1980
+    """
+    S0cc = S0_ii(q, C)
+    S0dd = S0_ii(q, D)
+    S0cd = S0_ij(q, [C, D])
+    Scc = (S0cc*S0dd - S0cd**2)/((S0cc+S0dd + 2*S0cd)-2*Kcd*(S0cc+S0dd-2*S0cd))
+    rhocd = drho(C,D)
+    Iq = rhocd**2 * Scc
+    return Iq
+
+def matrix_form(q, case_num, polys,
+                Kab=0., Kac=0., Kad=0., Kbc=0., Kbd=0., Kcd=0.):
+    if case_num < 2:
+        C, D = polys
+        A = B = D
+    elif case_num < 5:
+        B, C, D = polys
+        A = D
+    else:
+        A, B, C, D = polys
+
+    rho = np.matrix([[drho(p, D)] for p in (A,B,C)])
+    S0aa = S0_ii(q, A)
+    S0bb = S0_ii(q, B)
+    S0cc = S0_ii(q, C)
+    S0ab = S0_ij(q, [A, B])
+    S0ac = S0_ij(q, [A, B, C])
+    S0bc = S0_ij(q, [B, C])
+    if case_num == 4: # No a-c interaction in triblock copolymer
+        S0ac *= 0.0
+    elif case_num == 9: # No a-c or a-d interaction in tetrablock copolymer
+        S0ac *= 0.0
+
+    S0 = np.array([[S0aa, S0ab, S0ac], [S0ab, S0bb, S0bc], [S0ac, S0bc, S0cc]])
+    S0dd = S0_ii(q, D)
+    vaa = 1./S0dd - 2*Kad
+    vbb = 1./S0dd - 2*Kbd
+    vcc = 1./S0dd - 2*Kcd
+    vab = 1./S0dd + Kab - Kad - Kbd
+    vac = 1./S0dd + Kac - Kad - Kcd
+    vbc = 1./S0dd + Kbc - Kbd - Kcd
+    v = np.array([[vaa, vab, vac], [vab, vbb, vbc], [vac, vbc, vcc]])
+
+    Iq = np.empty_like(q)
+    for k, qk in enumerate(q):
+        S0_k = S0[:,:,k].M
+        v_k = v[:,:,k].M
+        S_k = np.linalg.inv(1 + S0_k*v_k)*S0_k
+        Iq[k] = rho.T * S_k * rho
+
+def build_pars(case_num, polys, **interactions):
+    def set_poly(x, poly):
+        pars["N"+x] = poly.n
+        pars["Phi"+x] = poly.phi
+        pars["v"+x] = poly.v
+        pars["b"+x] = poly.a
+        pars["L"+x] = poly.b
+    pars = interactions.copy()
+    pars["case_num"] = case_num
+    polys = list(reversed(polys))
+    if len(polys) >= 4: set_poly("a",polys[3])
+    if len(polys) >= 3: set_poly("b",polys[2])
+    if len(polys) >= 2: set_poly("c",polys[1])
+    if len(polys) >= 1: set_poly("d",polys[0])
+    return pars
+
+def sasmodels_rpa(q, pars):
+    from sasmodels.models import rpa
+    from sasmodels.core import load_model
+    from sasmodels.direct_model import DirectModel
+    from sasmodels.data import empty_data1D
+    data = empty_data1D(q, resolution=0.0)
+    model = load_model(rpa, dtype="double", platform="dll")
+    #model = load_model(rpa, dtype="single", platform="ocl")
+    M = DirectModel(data, model)
+    return M(**pars)
+
+def sasview_rpa(q, pars):
+    from sasmodels.models import rpa
+    from sasmodels.compare import eval_sasview
+    from sasmodels.data import empty_data1D
+    data = empty_data1D(q, resolution=0.0)
+    M = eval_sasview(rpa, data)
+    return M(**pars)
+
+def demo():
+    import sys
+    case_num = 0 if len(sys.argv) < 2 else int(sys.argv[1])
+
+    B = Polymer(n=525,phi=0.57,v=97.5,b=-4.99,a=8)
+    C = Polymer(n=525,phi=0.57,v=97.5,b=-4.99,a=8)
+    D = Polymer(n=1105,phi=0.43,v=81.9,b=53.1,a=2)
+    q = np.logspace(-4,-1,400)
+    #q = np.array([0.1])
+    Kab=Kac=Kad=0.0
+    Kcd = 0.0106*0.0035 - 1.84e-5
+    Kbd = Kcd + 2e-4
+    Kbc = (Kcd + 1e-4)*0.5
+    K = dict(Kab=Kab,Kac=Kac,Kad=Kad,Kbc=Kbc,Kbd=Kbd,Kcd=Kcd)
+    if case_num == 0:
+        Iq = binary_blend(q, C, D, Kcd)
+    elif case_num == 1:
+        Iq = diblock_copolymer(q, C, D, Kcd)
+    elif case_num == 2:
+        Iq = ternary_blend(q, B, C, D, Kbc, Kbd, Kcd)
+    else:
+        raise ValueError("Case %d not implmented"%case_num)
+
+    pars = build_pars(case_num, [B, C, D], **K)
+    print "eval sasmodels"
+    Iq_sasmodels = sasmodels_rpa(q, pars)
+    print "eval sasview"
+    Iq_sasview = sasview_rpa(q, pars)
+    print 1./Iq[0], 1./Iq_sasmodels[0], 1./Iq_sasview[0]
+
+    #return
+    import pylab
+    pylab.subplot(121)
+    pylab.loglog(q, Iq, label='direct')
+    pylab.loglog(q, Iq_sasmodels, label='sasmodels')
+    pylab.loglog(q, Iq_sasview, label='sasview')
+    pylab.legend()
+    pylab.subplot(122)
+    pylab.loglog(q, abs(Iq_sasview - Iq_sasmodels)/Iq_sasmodels, label='sasview-sasmodels')
+    pylab.loglog(q, abs(Iq_sasmodels - Iq)/Iq, label='sasmodels-direct')
+    pylab.legend()
+    pylab.show()
+
+if __name__ == "__main__":
+    demo()
Index: explore/sinc.py
===================================================================
--- explore/sinc.py	(revision e6f14104e82489f43d4acc1eeda0d50ad5408fb6)
+++ explore/sinc.py	(revision e6f14104e82489f43d4acc1eeda0d50ad5408fb6)
@@ -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/sph_j1c.py
===================================================================
--- explore/sph_j1c.py	(revision e6f14104e82489f43d4acc1eeda0d50ad5408fb6)
+++ explore/sph_j1c.py	(revision e6f14104e82489f43d4acc1eeda0d50ad5408fb6)
@@ -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: sasmodels/models/fuzzy_sphere.py
===================================================================
--- sasmodels/models/fuzzy_sphere.py	(revision 8dca856265977fe3d9e69f78610587271bb2ed27)
+++ sasmodels/models/fuzzy_sphere.py	(revision 8dca856265977fe3d9e69f78610587271bb2ed27)
@@ -0,0 +1,140 @@
+r"""
+For information about polarised and magnetic scattering, click here_.
+
+.. _here: polar_mag_help.html
+
+Definition
+----------
+
+The scattering intensity *I(q)* is calculated as:
+
+.. math::
+    I(q) = \frac{scale}{V}(\Delta \rho)^2 A^2(q) S(q) +\text{background}
+
+
+where the amplitude *A(q)* is given as the typical sphere scattering convoluted
+with a Gaussian to get a gradual drop-off in the scattering length density:
+
+.. math::
+
+    A(q) = \frac{3\left[\sin(qR) - qR \cos(qR)\right]}{(qR)^3}
+           \exp\left(\frac{-(o_{fuzzy}q)^2}{2}\right)
+
+Here *|A(q)|*:sup:`2`\  is the form factor, *P(q)*. The scale is equivalent to the
+volume fraction of spheres, each of volume, *V*\. Contrast (|drho|) is the
+difference of scattering length densities of the sphere and the surrounding
+solvent.
+
+Poly-dispersion in radius and in fuzziness is provided for.
+
+
+
+From the reference:
+
+  The "fuzziness" of the interface is defined by the parameter
+  |sigma| :sub:`fuzzy`\ . The particle radius *R* represents the radius of the
+  particle where the scattering length density profile decreased to 1/2 of the
+  core density. The |sigma| :sub:`fuzzy`\ is the width of the smeared particle
+  surface; i.e., the standard deviation from the average height of the fuzzy
+  interface. The inner regions of the microgel that display a higher density
+  are described by the radial box profile extending to a radius of
+  approximately *Rbox* ~ *R* - 2\ |sigma|\ . The profile approaches zero as
+  *Rsans* ~ *R* + 2\ |sigma|\ .
+
+For 2D data: The 2D scattering intensity is calculated in the same way as 1D,
+where the *q* vector is defined as
+
+.. math::
+
+    q = \sqrt{{q_x}^2 + {q_y}^2}
+
+
+.. figure:: img/fuzzy_sphere.jpg
+
+     This example dataset is produced by running the FuzzySphereModel,
+     using 200 data points, *qmin* = 0.001 -1,
+     *qmax* = 0.7 |Ang^-1|, background = 0.001 |cm^-1| and the default values.
+
+
+
+References
+----------
+
+M Stieger, J. S Pedersen, P Lindner, W Richtering, *Langmuir*,
+20 (2004) 7283-7292
+"""
+
+from numpy import inf
+
+name = "fuzzy_sphere"
+title = "Scattering from spherical particles with a fuzzy surface."
+description = """\
+scale: scale factor times volume fraction,
+or just volume fraction for absolute scale data
+radius: radius of the solid sphere
+fuzziness = the STD of the height of fuzzy interfacial
+thickness (ie., so-called interfacial roughness)
+sld: the SLD of the sphere
+solvend_sld: the SLD of the solvent
+background: incoherent background
+Note: By definition, this function works only when fuzziness << radius.
+"""
+category = "shape:sphere"
+
+# pylint: disable=bad-whitespace,line-too-long
+# ["name", "units", default, [lower, upper], "type","description"],
+parameters = [["sld",         "1e-6/Ang^2",  1, [-inf, inf], "",       "Layer scattering length density"],
+              ["solvent_sld", "1e-6/Ang^2",  3, [-inf, inf], "",       "Solvent scattering length density"],
+              ["radius",      "Ang",        60, [0, inf],    "volume", "Sphere radius"],
+              ["fuzziness",   "Ang",        10, [0, inf],    "",       "The STD of the height of fuzzy interfacial"],
+             ]
+# pylint: enable=bad-whitespace,line-too-long
+
+source = ["lib/sph_j1c.c"]
+
+# No volume normalization despite having a volume parameter
+# This should perhaps be volume normalized?
+form_volume = """
+    return 1.333333333333333*M_PI*radius*radius*radius;
+    """
+
+Iq = """
+    const double qr = q*radius;
+    const double bes = sph_j1c(qr);
+    const double qf = q*fuzziness;
+    const double fq = bes * (sld - solvent_sld) * form_volume(radius) * exp(-0.5*qf*qf);
+    return 1.0e-4*fq*fq;
+    """
+
+Iqxy = """
+    // never called since no orientation or magnetic parameters.
+    //return -1.0;
+    return Iq(sqrt(qx*qx + qy*qy), sld, solvent_sld, radius, fuzziness);
+    """
+
+def ER(radius):
+    """
+    Return radius
+    """
+    return radius
+
+# VR defaults to 1.0
+
+demo = dict(scale=1, background=0.001,
+            sld=1, solvent_sld=3,
+            radius=60,
+            fuzziness=10,
+            radius_pd=.2, radius_pd_n=45,
+            fuzziness_pd=.2, fuzziness_pd_n=0)
+
+oldname = "FuzzySphereModel"
+oldpars = dict(sld='sldSph', solvent_sld='sldSolv', radius='radius', fuzziness='fuzziness')
+
+
+tests = [
+    # Accuracy tests based on content in test/utest_models_new1_3.py
+    #[{'background': 0.001}, 1.0, 0.001],
+
+    [{}, 0.00301005, 359.2315],
+
+    ]
Index: sasmodels/models/lib/lanczos_gamma.c
===================================================================
--- sasmodels/models/lib/lanczos_gamma.c	(revision 5c1d3414d3ff3bc81b5e9cfefb15dded91ee5d3e)
+++ sasmodels/models/lib/lanczos_gamma.c	(revision e6f14104e82489f43d4acc1eeda0d50ad5408fb6)
@@ -2,5 +2,5 @@
 double lanczos_gamma(double q)
 {
-    // Lanczos approximation to the Gamma function.
+    // Lanczos approximation to the Log Gamma function.
 
     double x,y,tmp,ser;
Index: sasmodels/models/lib/sph_j1c.c
===================================================================
--- sasmodels/models/lib/sph_j1c.c	(revision c437dbb3bdf21e923e4309e102d0e76cc8b331db)
+++ sasmodels/models/lib/sph_j1c.c	(revision e6f14104e82489f43d4acc1eeda0d50ad5408fb6)
@@ -5,66 +5,52 @@
 * Note that the values differ from sasview ~ 5e-12 rather than 5e-14, but
 * in this case it is likely cancellation errors in the original expression
-* using double precision that are the source.  Single precision only
-* requires the first 3 terms.  Double precision requires the 4th term.
-* The fifth term is not needed, and is commented out.
-* Taylor expansion:
-*      1.0 + q2*(-3./30. + q2*(3./840.))+ q2*(-3./45360. + q2*(3./3991680.))))
+* using double precision that are the source.
 */
+
+// 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 **
+//   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.
+//
+// See explore/sph_j1c.py for code to explore these ranges.
+
+// Use 4th order series
+#if FLOAT_SIZE>4
+#define SPH_J1C_CUTOFF 0.1
+#else
+#define SPH_J1C_CUTOFF 0.7
+#endif
 
 double sph_j1c(double q);
 double sph_j1c(double q)
 {
-    const double q2 = q*q;
-    double sin_q, cos_q;
-
-    SINCOS(q, sin_q, cos_q);
-
-#if FLOAT_SIZE>4
-// For double precision, choose a cutoff of 0.18, which is the lower limit
-// for the trig expression to 14 digits.  If only 12 digits are needed, then
-// only 4 terms of the Taylor expansion are needed.
-#define CUTOFF 0.18
-#else
-// For single precision, choose a cutoff halfway between the single precision
-// lower limit for the trig expression of 1.03 and the upper limit of 1.3
-// for the Taylor series expansion with 5 terms (or 9 if you count the zeros).
-// For 5 digits of precision, you can drop two terms of the taylor series
-// and choose a cutoff of 0.3.
-#define CUTOFF 1.15
-#endif
-
-    const double bessel = (q < CUTOFF)
-        ? (1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360. + q2*(3./3991680.)))))
-        : 3.0*(sin_q/q - cos_q)/q2;
-    return bessel;
-
-/*
-    // Code to test various expressions
-    // tested using sympy.mpmath.mp with
-    //    mp.dps = 50 (which is good for q down to 1e-6)
-    //    def j1c(x): return 3*(mp.sin(x)/x - mp.cos(x))/x**2
-    double ret;
-    if (sizeof(q2) > 8) {
-
-        ret = 3.0*(sinl(q)/q - cosl(q))/q2;
-printf("%.15Lf %.15Lg\n", q, ret);
-    //} else if (q < 0.384038453352533) {
-    //} else if (q < 0.18) {
-    } else if (q < 18e-10) {
-        // NB: all are good to 5 digits below 0.18f
-        // good is defined as 14 digits for double and 7 for single
-        //ret = 1.0 + q2*q2*(3./840.) - q2*(3./30.); // good below 0.02d 0.34f
-        ret = square((1. + 3./5600.*q2*q2) - q2/20.); // good below 0.03d 0.08f
-        //ret = 1.0 - q2*(3./30.); // good below 0.02d 0.07f
-        //ret = 1.0 + q2*(-3./30. + q2*(3./840.)); // good below 0.02d 0.34f
-        //ret = 1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360.))); // good below 0.1d 0.8f, 12 digits below 0.18d
-        //ret = 1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360. + q2*(3./3991680.)))); // good below 0.18d 1.3f
-printf("%.15g %.15g\n", q, ret);
+    if (q < SPH_J1C_CUTOFF) {
+        const double q2 = q*q;
+        return (1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360.))));// + q2*(3./3991680.)))));
     } else {
-        // NB: can use a cutoff of 0.1f if the goal is 5 digits rather than 7
-        ret = 3.0*(sin_q/q - cos_q)/q2; // good above 0.18d, 1.03f
-printf("%.15g %.15g\n", q, ret);
+        double sin_q, cos_q;
+        SINCOS(q, sin_q, cos_q);
+        return 3.0*(sin_q/q - cos_q)/(q*q);
     }
-    return ret;
-*/
 }
Index: st/rpa.py
===================================================================
--- test/rpa.py	(revision 82c299fa51ab091ee8ead0a235d4177ccfeecb99)
+++ 	(revision )
@@ -1,202 +1,0 @@
-from collections import namedtuple
-import numpy as np
-from numpy import sqrt, exp, expm1
-
-AVOGADRO = 6.022e23
-
-Polymer = namedtuple("Polymer", "n phi v a b".split())
-def E(q, poly):
-    qrsq = (q*Rg(poly))**2
-    retval = exp(-qrsq)
-    return retval
-
-def F(q, poly):
-    qrsq = (q*Rg(poly))**2
-    retval = -expm1(-qrsq)/qrsq
-    return retval
-
-def P_ii(q, poly):
-    qrsq = (q*Rg(poly))**2
-    retval = 2 * (expm1(-qrsq) + qrsq)/qrsq**2
-    return retval
-
-def P_ij(q, poly_list):
-    i, j = poly_list[0], poly_list[-1]
-    retval = F(q, i) * np.prod([E(q,p) for p in poly_list[1:-1]]) * F(q, j)
-    return retval
-
-def Rg(poly):
-    return sqrt(poly.n/6.)*poly.a
-
-def S0_ii(q, poly):
-    retval = poly.n*poly.phi*poly.v*P_ii(q, poly)
-    return retval
-
-def S0_ij(q, poly_list):
-    i,j = poly_list[0], poly_list[-1]
-    retval = sqrt(i.n*i.phi*i.v*j.n*j.phi*j.v) * P_ij(q, poly_list)
-    return retval
-
-def drho(poly, base):
-    return (poly.b/poly.v -  base.b/base.v)*1e-13*sqrt(AVOGADRO)
-
-def binary_blend(q, C, D, Kcd):
-    """
-    de Gennes, 1979
-    """
-    S0cc = S0_ii(q, C)
-    S0dd = S0_ii(q, D)
-    vcc = 1/S0dd - 2*Kcd #/v0
-    Scc = S0cc/(1 + vcc*S0cc)
-    rhocd = drho(C,D)
-    Iq = rhocd**2 * Scc
-    return Iq
-
-def ternary_blend(q, B, C, D, Kbc, Kbd, Kcd):
-    S0bb = S0_ii(q, B)
-    S0cc = S0_ii(q, C)
-    S0dd = S0_ii(q, D)
-    vbb = 1/S0dd - 2*Kbd
-    vcc = 1/S0dd - 2*Kcd
-    vbc = 1/S0dd + Kbd - Kbc - Kcd
-    rhobd = drho(B,D)
-    rhocd = drho(C,D)
-    det = (1+vbb*S0bb)*(1+vcc*S0cc) - vbc**2*S0bb*S0cc
-    Sbb = S0bb*(1+vcc*S0cc)/det
-    Scc = S0cc*(1+vbb*S0bb)/det
-    Sbc = -S0bb*vbc*S0cc/det
-    Iq = rhobd**2*Sbb + rhocd**2*Scc + 2*rhobd*rhocd*Sbc
-    return Iq
-
-def diblock_copolymer(q, C, D, Kcd):
-    """
-    Leibler, 1980
-    """
-    S0cc = S0_ii(q, C)
-    S0dd = S0_ii(q, D)
-    S0cd = S0_ij(q, [C, D])
-    Scc = (S0cc*S0dd - S0cd**2)/((S0cc+S0dd + 2*S0cd)-2*Kcd*(S0cc+S0dd-2*S0cd))
-    rhocd = drho(C,D)
-    Iq = rhocd**2 * Scc
-    return Iq
-
-def matrix_form(q, case_num, polys,
-                Kab=0., Kac=0., Kad=0., Kbc=0., Kbd=0., Kcd=0.):
-    if case_num < 2:
-        C, D = polys
-        A = B = D
-    elif case_num < 5:
-        B, C, D = polys
-        A = D
-    else:
-        A, B, C, D = polys
-
-    rho = np.matrix([[drho(p, D)] for p in (A,B,C)])
-    S0aa = S0_ii(q, A)
-    S0bb = S0_ii(q, B)
-    S0cc = S0_ii(q, C)
-    S0ab = S0_ij(q, [A, B])
-    S0ac = S0_ij(q, [A, B, C])
-    S0bc = S0_ij(q, [B, C])
-    if case_num == 4: # No a-c interaction in triblock copolymer
-        S0ac *= 0.0
-    elif case_num == 9: # No a-c or a-d interaction in tetrablock copolymer
-        S0ac *= 0.0
-
-    S0 = np.array([[S0aa, S0ab, S0ac], [S0ab, S0bb, S0bc], [S0ac, S0bc, S0cc]])
-    S0dd = S0_ii(q, D)
-    vaa = 1./S0dd - 2*Kad
-    vbb = 1./S0dd - 2*Kbd
-    vcc = 1./S0dd - 2*Kcd
-    vab = 1./S0dd + Kab - Kad - Kbd
-    vac = 1./S0dd + Kac - Kad - Kcd
-    vbc = 1./S0dd + Kbc - Kbd - Kcd
-    v = np.array([[vaa, vab, vac], [vab, vbb, vbc], [vac, vbc, vcc]])
-
-    Iq = np.empty_like(q)
-    for k, qk in enumerate(q):
-        S0_k = S0[:,:,k].M
-        v_k = v[:,:,k].M
-        S_k = np.linalg.inv(1 + S0_k*v_k)*S0_k
-        Iq[k] = rho.T * S_k * rho
-
-def build_pars(case_num, polys, **interactions):
-    def set_poly(x, poly):
-        pars["N"+x] = poly.n
-        pars["Phi"+x] = poly.phi
-        pars["v"+x] = poly.v
-        pars["b"+x] = poly.a
-        pars["L"+x] = poly.b
-    pars = interactions.copy()
-    pars["case_num"] = case_num
-    polys = list(reversed(polys))
-    if len(polys) >= 4: set_poly("a",polys[3])
-    if len(polys) >= 3: set_poly("b",polys[2])
-    if len(polys) >= 2: set_poly("c",polys[1])
-    if len(polys) >= 1: set_poly("d",polys[0])
-    return pars
-
-def sasmodels_rpa(q, pars):
-    from sasmodels.models import rpa
-    from sasmodels.core import load_model
-    from sasmodels.direct_model import DirectModel
-    from sasmodels.data import empty_data1D
-    data = empty_data1D(q, resolution=0.0)
-    model = load_model(rpa, dtype="double", platform="dll")
-    #model = load_model(rpa, dtype="single", platform="ocl")
-    M = DirectModel(data, model)
-    return M(**pars)
-
-def sasview_rpa(q, pars):
-    from sasmodels.models import rpa
-    from sasmodels.compare import eval_sasview
-    from sasmodels.data import empty_data1D
-    data = empty_data1D(q, resolution=0.0)
-    M = eval_sasview(rpa, data)
-    return M(**pars)
-
-def demo():
-    import sys
-    case_num = 0 if len(sys.argv) < 2 else int(sys.argv[1])
-
-    B = Polymer(n=525,phi=0.57,v=97.5,b=-4.99,a=8)
-    C = Polymer(n=525,phi=0.57,v=97.5,b=-4.99,a=8)
-    D = Polymer(n=1105,phi=0.43,v=81.9,b=53.1,a=2)
-    q = np.logspace(-4,-1,400)
-    #q = np.array([0.1])
-    Kab=Kac=Kad=0.0
-    Kcd = 0.0106*0.0035 - 1.84e-5
-    Kbd = Kcd + 2e-4
-    Kbc = (Kcd + 1e-4)*0.5
-    K = dict(Kab=Kab,Kac=Kac,Kad=Kad,Kbc=Kbc,Kbd=Kbd,Kcd=Kcd)
-    if case_num == 0:
-        Iq = binary_blend(q, C, D, Kcd)
-    elif case_num == 1:
-        Iq = diblock_copolymer(q, C, D, Kcd)
-    elif case_num == 2:
-        Iq = ternary_blend(q, B, C, D, Kbc, Kbd, Kcd)
-    else:
-        raise ValueError("Case %d not implmented"%case_num)
-
-    pars = build_pars(case_num, [B, C, D], **K)
-    print "eval sasmodels"
-    Iq_sasmodels = sasmodels_rpa(q, pars)
-    print "eval sasview"
-    Iq_sasview = sasview_rpa(q, pars)
-    print 1./Iq[0], 1./Iq_sasmodels[0], 1./Iq_sasview[0]
-
-    #return
-    import pylab
-    pylab.subplot(121)
-    pylab.loglog(q, Iq, label='direct')
-    pylab.loglog(q, Iq_sasmodels, label='sasmodels')
-    pylab.loglog(q, Iq_sasview, label='sasview')
-    pylab.legend()
-    pylab.subplot(122)
-    pylab.loglog(q, abs(Iq_sasview - Iq_sasmodels)/Iq_sasmodels, label='sasview-sasmodels')
-    pylab.loglog(q, abs(Iq_sasmodels - Iq)/Iq, label='sasmodels-direct')
-    pylab.legend()
-    pylab.show()
-
-if __name__ == "__main__":
-    demo()
