Index: explore/J1.py
===================================================================
--- explore/J1.py	(revision a5af4e1ca0acba8ff274ff823f33faa5a25b0704)
+++ explore/J1.py	(revision a5af4e1ca0acba8ff274ff823f33faa5a25b0704)
@@ -0,0 +1,103 @@
+r"""
+Show numerical precision of $2 J_1(x)/x$.
+"""
+
+import numpy as np
+#from sympy.mpmath import mp
+from 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)
+    #ret = calculator(background=0)
+    #print ret
+    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: sasmodels/models/bessel.py
===================================================================
--- sasmodels/models/bessel.py	(revision cbd37a7e9b2cd2085f3e1d42c1ab4a48c4db56a7)
+++ sasmodels/models/bessel.py	(revision a5af4e1ca0acba8ff274ff823f33faa5a25b0704)
@@ -78,5 +78,5 @@
 
 Iq = """
-    return 2.0*j1(q)/q;
+    return J1(q);
     """
 
Index: sasmodels/models/lib/j1_cephes.c
===================================================================
--- sasmodels/models/lib/j1_cephes.c	(revision e2af2a953458fc0c3f5a1c0c5754c8455115eac8)
+++ sasmodels/models/lib/j1_cephes.c	(revision a5af4e1ca0acba8ff274ff823f33faa5a25b0704)
@@ -39,7 +39,7 @@
 Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
 */
-double j1(double );
+double J1(double );
 
-double j1(double x) {
+double J1(double x) {
 
 //Cephes double pression function
@@ -48,77 +48,8 @@
     double w, z, p, q, xn;
 
-    const double DR1 = 5.78318596294678452118E0;
-    const double DR2 = 3.04712623436620863991E1;
     const double Z1 = 1.46819706421238932572E1;
     const double Z2 = 4.92184563216946036703E1;
     const double THPIO4 =  2.35619449019234492885;
     const double SQ2OPI = 0.79788456080286535588;
-
-    double RP[8] = {
-    -8.99971225705559398224E8,
-    4.52228297998194034323E11,
-    -7.27494245221818276015E13,
-    3.68295732863852883286E15,
-    0.0,
-    0.0,
-    0.0,
-    0.0
-    };
-
-    double RQ[8] = {
-    /* 1.00000000000000000000E0,*/
-    6.20836478118054335476E2,
-    2.56987256757748830383E5,
-    8.35146791431949253037E7,
-    2.21511595479792499675E10,
-    4.74914122079991414898E12,
-    7.84369607876235854894E14,
-    8.95222336184627338078E16,
-    5.32278620332680085395E18,
-    };
-
-    double PP[8] = {
-    7.62125616208173112003E-4,
-    7.31397056940917570436E-2,
-    1.12719608129684925192E0,
-    5.11207951146807644818E0,
-    8.42404590141772420927E0,
-    5.21451598682361504063E0,
-    1.00000000000000000254E0,
-    0.0,
-    };
-    double PQ[8] = {
-    5.71323128072548699714E-4,
-    6.88455908754495404082E-2,
-    1.10514232634061696926E0,
-    5.07386386128601488557E0,
-    8.39985554327604159757E0,
-    5.20982848682361821619E0,
-    9.99999999999999997461E-1,
-    0.0,
-    };
-
-    double QP[8] = {
-    5.10862594750176621635E-2,
-    4.98213872951233449420E0,
-    7.58238284132545283818E1,
-    3.66779609360150777800E2,
-    7.10856304998926107277E2,
-    5.97489612400613639965E2,
-    2.11688757100572135698E2,
-    2.52070205858023719784E1,
-    };
-
-    double QQ[8] = {
-    /* 1.00000000000000000000E0,*/
-    7.42373277035675149943E1,
-    1.05644886038262816351E3,
-    4.98641058337653607651E3,
-    9.56231892404756170795E3,
-    7.99704160447350683650E3,
-    2.82619278517639096600E3,
-    3.36093607810698293419E2,
-    0.0,
-    };
 
     w = x;
@@ -129,5 +60,5 @@
 	{
 	    z = x * x;
-	    w = polevl( z, RP, 3 ) / p1evl( z, RQ, 8 );
+	    w = polevlRP( z, 3 ) / p1evlRQ( z, 8 );
 	    w = w * x * (z - Z1) * (z - Z2);
 	    return( w );
@@ -136,6 +67,8 @@
     w = 5.0/x;
     z = w * w;
-    p = polevl( z, PP, 6)/polevl( z, PQ, 6 );
-    q = polevl( z, QP, 7)/p1evl( z, QQ, 7 );
+
+    p = polevlPP( z, 6)/polevlPQ( z, 6 );
+    q = polevlQP( z, 7)/p1evlQQ( z, 7 );
+
     xn = x - THPIO4;
 
@@ -155,37 +88,4 @@
 
 
-    double JP[8] = {
-        -4.878788132172128E-009,
-        6.009061827883699E-007,
-        -4.541343896997497E-005,
-        1.937383947804541E-003,
-        -3.405537384615824E-002,
-        0.0,
-        0.0,
-        0.0
-    };
-
-    double MO1[8] = {
-        6.913942741265801E-002,
-        -2.284801500053359E-001,
-        3.138238455499697E-001,
-        -2.102302420403875E-001,
-        5.435364690523026E-003,
-        1.493389585089498E-001,
-        4.976029650847191E-006,
-        7.978845453073848E-001
-    };
-
-    double PH1[8] = {
-        -4.497014141919556E+001,
-        5.073465654089319E+001,
-        -2.485774108720340E+001,
-        7.222973196770240E+000,
-        -1.544842782180211E+000,
-        3.503787691653334E-001,
-        -1.637986776941202E-001,
-        3.749989509080821E-001
-    };
-
     xx = x;
     if( xx < 0 )
@@ -195,5 +95,5 @@
 	{
 	    z = xx * xx;
-	    p = (z-Z1) * xx * polevl( z, JP, 4 );
+	    p = (z-Z1) * xx * polevlJP( z, 4 );
 	    return( p );
 	}
@@ -202,7 +102,7 @@
     w = sqrt(q);
 
-    p = w * polevl( q, MO1, 7);
+    p = w * polevlMO1( q, 7);
     w = q*q;
-    xn = q * polevl( w, PH1, 7) - THPIO4F;
+    xn = q * polevlPH1( w, 7) - THPIO4F;
     p = p * cos(xn + xx);
 
@@ -210,3 +110,2 @@
 #endif
 }
-
Index: sasmodels/models/lib/polevl.c
===================================================================
--- sasmodels/models/lib/polevl.c	(revision e2af2a953458fc0c3f5a1c0c5754c8455115eac8)
+++ sasmodels/models/lib/polevl.c	(revision a5af4e1ca0acba8ff274ff823f33faa5a25b0704)
@@ -51,20 +51,218 @@
 */
 
-double polevl( double x, double coef[8], int N );
-double p1evl( double x, double coef[8], int N );
-
-double polevl( double x, double coef[8], int N ) {
-
-    int i = 0;
-    double ans = coef[i];
-
-    while (i < N) {
-        i++;
-        ans = ans * x + coef[i];
-    }
-
-    return ans ;
-
-}
+double polevlRP(double x, int N ) {
+
+    double coef[8] = {
+    -8.99971225705559398224E8,
+    4.52228297998194034323E11,
+    -7.27494245221818276015E13,
+    3.68295732863852883286E15,
+    0.0,
+    0.0,
+    0.0,
+    0.0 };
+
+    int i = 0;
+    double ans = coef[i];
+
+    while (i < N) {
+        i++;
+        ans = ans * x + coef[i];
+    }
+    return ans ;
+}
+
+double polevlRQ(double x, int N ) {
+
+    double coef[8] = {
+        6.20836478118054335476E2,
+        2.56987256757748830383E5,
+        8.35146791431949253037E7,
+        2.21511595479792499675E10,
+        4.74914122079991414898E12,
+        7.84369607876235854894E14,
+        8.95222336184627338078E16,
+        5.32278620332680085395E18
+    };
+
+    int i = 0;
+    double ans = coef[i];
+
+    while (i < N) {
+        i++;
+        ans = ans * x + coef[i];
+    }
+    return ans ;
+}
+
+double polevlPP(double x, int N ) {
+
+    double coef[8] = {
+    7.62125616208173112003E-4,
+    7.31397056940917570436E-2,
+    1.12719608129684925192E0,
+    5.11207951146807644818E0,
+    8.42404590141772420927E0,
+    5.21451598682361504063E0,
+    1.00000000000000000254E0,
+    0.0} ;
+
+    int i = 0;
+    double ans = coef[i];
+
+    while (i < N) {
+        i++;
+        ans = ans * x + coef[i];
+    }
+    return ans ;
+}
+
+double polevlPQ(double x, int N ) {
+//4: PQ
+    double coef[8] = {
+    5.71323128072548699714E-4,
+    6.88455908754495404082E-2,
+    1.10514232634061696926E0,
+    5.07386386128601488557E0,
+    8.39985554327604159757E0,
+    5.20982848682361821619E0,
+    9.99999999999999997461E-1,
+    0.0 };
+
+    int i = 0;
+    double ans = coef[i];
+
+    while (i < N) {
+        i++;
+        ans = ans * x + coef[i];
+    }
+    return ans ;
+}
+
+double polevlQP(double x, int N ) {
+    double coef[8] = {
+    5.10862594750176621635E-2,
+    4.98213872951233449420E0,
+    7.58238284132545283818E1,
+    3.66779609360150777800E2,
+    7.10856304998926107277E2,
+    5.97489612400613639965E2,
+    2.11688757100572135698E2,
+    2.52070205858023719784E1 };
+
+    int i = 0;
+    double ans = coef[i];
+
+    while (i < N) {
+        i++;
+        ans = ans * x + coef[i];
+    }
+    return ans ;
+
+}
+
+double polevlQQ(double x, int N ) {
+    double coef[8] = {
+    7.42373277035675149943E1,
+    1.05644886038262816351E3,
+    4.98641058337653607651E3,
+    9.56231892404756170795E3,
+    7.99704160447350683650E3,
+    2.82619278517639096600E3,
+    3.36093607810698293419E2,
+    0.0 };
+
+    int i = 0;
+    double ans = coef[i];
+
+    while (i < N) {
+        i++;
+        ans = ans * x + coef[i];
+    }
+    return ans ;
+
+}
+
+double polevlJP( double x, int N ) {
+    double coef[8] = {
+        -4.878788132172128E-009,
+        6.009061827883699E-007,
+        -4.541343896997497E-005,
+        1.937383947804541E-003,
+        -3.405537384615824E-002,
+        0.0,
+        0.0,
+        0.0
+    };
+
+    int i = 0;
+    double ans = coef[i];
+
+    while (i < N) {
+        i++;
+        ans = ans * x + coef[i];
+    }
+    return ans ;
+
+}
+
+double polevlMO1( double x, int N ) {
+    double coef[8] = {
+        6.913942741265801E-002,
+        -2.284801500053359E-001,
+        3.138238455499697E-001,
+        -2.102302420403875E-001,
+        5.435364690523026E-003,
+        1.493389585089498E-001,
+        4.976029650847191E-006,
+        7.978845453073848E-001
+    };
+
+    int i = 0;
+    double ans = coef[i];
+
+    while (i < N) {
+        i++;
+        ans = ans * x + coef[i];
+    }
+    return ans ;
+}
+
+double polevlPH1( double x, int N ) {
+
+    double coef[8] = {
+        -4.497014141919556E+001,
+        5.073465654089319E+001,
+        -2.485774108720340E+001,
+        7.222973196770240E+000,
+        -1.544842782180211E+000,
+        3.503787691653334E-001,
+        -1.637986776941202E-001,
+        3.749989509080821E-001
+    };
+
+    int i = 0;
+    double ans = coef[i];
+
+    while (i < N) {
+        i++;
+        ans = ans * x + coef[i];
+    }
+    return ans ;
+}
+
+/*double polevl( double x, double coef[8], int N ) {
+
+    int i = 0;
+    double ans = coef[i];
+
+    while (i < N) {
+        i++;
+        ans = ans * x + coef[i];
+    }
+
+    return ans ;
+
+}*/
 
 /*							p1evl()	*/
@@ -74,14 +272,218 @@
  */
 
-double p1evl( double x, double coef[8], int N ) {
-    int i=0;
-    double ans = x+coef[i];
-
-    while (i < N-1) {
-        i++;
-        ans = ans*x + coef[i];
-    }
-
-    return( ans );
-
-}
+double p1evlRP( double x, int N ) {
+    double coef[8] = {
+    -8.99971225705559398224E8,
+    4.52228297998194034323E11,
+    -7.27494245221818276015E13,
+    3.68295732863852883286E15,
+    0.0,
+    0.0,
+    0.0,
+    0.0 };
+
+    int i=0;
+    double ans = x+coef[i];
+
+    while (i < N-1) {
+        i++;
+        ans = ans*x + coef[i];
+    }
+
+    return( ans );
+
+}
+
+double p1evlRQ( double x, int N ) {
+//1: RQ
+    double coef[8] = {
+    6.20836478118054335476E2,
+    2.56987256757748830383E5,
+    8.35146791431949253037E7,
+    2.21511595479792499675E10,
+    4.74914122079991414898E12,
+    7.84369607876235854894E14,
+    8.95222336184627338078E16,
+    5.32278620332680085395E18};
+
+    int i=0;
+    double ans = x+coef[i];
+
+    while (i < N-1) {
+        i++;
+        ans = ans*x + coef[i];
+    }
+
+    return( ans );
+}
+
+double p1evlPP( double x, int N ) {
+//3 : PP
+    double coef[8] = {
+    7.62125616208173112003E-4,
+    7.31397056940917570436E-2,
+    1.12719608129684925192E0,
+    5.11207951146807644818E0,
+    8.42404590141772420927E0,
+    5.21451598682361504063E0,
+    1.00000000000000000254E0,
+    0.0};
+
+    int i=0;
+    double ans = x+coef[i];
+
+    while (i < N-1) {
+        i++;
+        ans = ans*x + coef[i];
+    }
+
+    return( ans );
+}
+
+double p1evlPQ( double x, int N ) {
+//4: PQ
+    double coef[8] = {
+    5.71323128072548699714E-4,
+    6.88455908754495404082E-2,
+    1.10514232634061696926E0,
+    5.07386386128601488557E0,
+    8.39985554327604159757E0,
+    5.20982848682361821619E0,
+    9.99999999999999997461E-1,
+    0.0};
+
+    int i=0;
+    double ans = x+coef[i];
+
+    while (i < N-1) {
+        i++;
+        ans = ans*x + coef[i];
+    }
+
+    return( ans );
+}
+
+double p1evlQP( double x, int N ) {
+//5: QP
+    double coef[8] = {
+    5.10862594750176621635E-2,
+    4.98213872951233449420E0,
+    7.58238284132545283818E1,
+    3.66779609360150777800E2,
+    7.10856304998926107277E2,
+    5.97489612400613639965E2,
+    2.11688757100572135698E2,
+    2.52070205858023719784E1 };
+
+    int i=0;
+    double ans = x+coef[i];
+
+    while (i < N-1) {
+        i++;
+        ans = ans*x + coef[i];
+    }
+
+    return( ans );
+}
+
+double p1evlQQ( double x, int N ) {
+    double coef[8] = {
+    7.42373277035675149943E1,
+    1.05644886038262816351E3,
+    4.98641058337653607651E3,
+    9.56231892404756170795E3,
+    7.99704160447350683650E3,
+    2.82619278517639096600E3,
+    3.36093607810698293419E2,
+    0.0};
+
+    int i=0;
+    double ans = x+coef[i];
+
+    while (i < N-1) {
+        i++;
+        ans = ans*x + coef[i];
+    }
+
+    return( ans );
+
+}
+
+double p1evlJP( double x, int N ) {
+    double coef[8] = {
+        -4.878788132172128E-009,
+        6.009061827883699E-007,
+        -4.541343896997497E-005,
+        1.937383947804541E-003,
+        -3.405537384615824E-002,
+        0.0,
+        0.0,
+        0.0};
+
+    int i=0;
+    double ans = x+coef[i];
+
+    while (i < N-1) {
+        i++;
+        ans = ans*x + coef[i];
+    }
+
+    return( ans );
+}
+
+double p1evlMO1( double x, int N ) {
+    double coef[8] = {
+        6.913942741265801E-002,
+        -2.284801500053359E-001,
+        3.138238455499697E-001,
+        -2.102302420403875E-001,
+        5.435364690523026E-003,
+        1.493389585089498E-001,
+        4.976029650847191E-006,
+        7.978845453073848E-001};
+
+    int i=0;
+    double ans = x+coef[i];
+
+    while (i < N-1) {
+        i++;
+        ans = ans*x + coef[i];
+    }
+
+    return( ans );
+
+}
+
+double p1evlPH1( double x, int N ) {
+    double coef[8] = {
+        -4.497014141919556E+001,
+        5.073465654089319E+001,
+        -2.485774108720340E+001,
+        7.222973196770240E+000,
+        -1.544842782180211E+000,
+        3.503787691653334E-001,
+        -1.637986776941202E-001,
+        3.749989509080821E-001};
+    int i=0;
+    double ans = x+coef[i];
+
+    while (i < N-1) {
+        i++;
+        ans = ans*x + coef[i];
+    }
+
+    return( ans );
+}
+
+/*double p1evl( double x, double coef[8], int N ) {
+    int i=0;
+    double ans = x+coef[i];
+
+    while (i < N-1) {
+        i++;
+        ans = ans*x + coef[i];
+    }
+
+    return( ans );
+
+}*/
