Index: sasmodels/models/core_shell_cylinder.c
===================================================================
--- sasmodels/models/core_shell_cylinder.c	(revision 994d77f9b713691e41106d1f9af9661dfb2dc419)
+++ sasmodels/models/core_shell_cylinder.c	(revision 3936ad3134459cadbcc9696bac14e0ffc9da869b)
@@ -11,5 +11,5 @@
 double _cyl(double twovd, double besarg, double siarg)
 {
-    const double bj = (besarg == 0.0 ? 0.5 : J1(besarg)/besarg);
+    const double bj = (besarg == 0.0 ? 0.5 : j1(besarg)/besarg);
     const double si = (siarg == 0.0 ? 1.0 : sin(siarg)/siarg);
     return twovd*si*bj;
Index: sasmodels/models/core_shell_cylinder.py
===================================================================
--- sasmodels/models/core_shell_cylinder.py	(revision f0aa7f88c13775739f654942a5ec343bf210d1f5)
+++ sasmodels/models/core_shell_cylinder.py	(revision 3936ad3134459cadbcc9696bac14e0ffc9da869b)
@@ -152,5 +152,5 @@
              ]
 
-source = ["lib/J1.c", "lib/gauss76.c", "core_shell_cylinder.c"]
+source = ["lib/polevl.c", "lib/j1d.c", "lib/gauss76.c", "core_shell_cylinder.c"]
 
 def ER(radius, thickness, length):
Index: sasmodels/models/lib/j0d.c
===================================================================
--- sasmodels/models/lib/j0d.c	(revision 3936ad3134459cadbcc9696bac14e0ffc9da869b)
+++ sasmodels/models/lib/j0d.c	(revision 3936ad3134459cadbcc9696bac14e0ffc9da869b)
@@ -0,0 +1,228 @@
+/*							j0.c
+ *
+ *	Bessel function of order zero
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, j0();
+ *
+ * y = j0( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns Bessel function of order zero of the argument.
+ *
+ * The domain is divided into the intervals [0, 5] and
+ * (5, infinity). In the first interval the following rational
+ * approximation is used:
+ *
+ *
+ *        2         2
+ * (w - r  ) (w - r  ) P (w) / Q (w)
+ *       1         2    3       8
+ *
+ *            2
+ * where w = x  and the two r's are zeros of the function.
+ *
+ * In the second interval, the Hankel asymptotic expansion
+ * is employed with two rational functions of degree 6/6
+ * and 7/7.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Absolute error:
+ * arithmetic   domain     # trials      peak         rms
+ *    DEC       0, 30       10000       4.4e-17     6.3e-18
+ *    IEEE      0, 30       60000       4.2e-16     1.1e-16
+ *
+ */
+
+/*							y0.c
+ *
+ *	Bessel function of the second kind, order zero
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, y0();
+ *
+ * y = y0( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns Bessel function of the second kind, of order
+ * zero, of the argument.
+ *
+ * The domain is divided into the intervals [0, 5] and
+ * (5, infinity). In the first interval a rational approximation
+ * R(x) is employed to compute
+ *   y0(x)  = R(x)  +   2 * log(x) * j0(x) / PI.
+ * Thus a call to j0() is required.
+ *
+ * In the second interval, the Hankel asymptotic expansion
+ * is employed with two rational functions of degree 6/6
+ * and 7/7.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *  Absolute error, when y0(x) < 1; else relative error:
+ *
+ * arithmetic   domain     # trials      peak         rms
+ *    DEC       0, 30        9400       7.0e-17     7.9e-18
+ *    IEEE      0, 30       30000       1.3e-15     1.6e-16
+ *
+ */
+
+
+/*
+Cephes Math Library Release 2.8:  June, 2000
+Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
+*/
+
+/* Note: all coefficients satisfy the relative error criterion
+ * except YP, YQ which are designed for absolute error. */
+
+double j0( double );
+
+double j0(double x)
+{
+    double w, z, p, q, xn;
+
+    const double TWOOPI = 6.36619772367581343075535E-1;
+    const double SQ2OPI = 7.9788456080286535587989E-1;
+    const double PIO4 = 7.85398163397448309616E-1;
+
+    const double PP[8] = {
+        7.96936729297347051624E-4,
+        8.28352392107440799803E-2,
+        1.23953371646414299388E0,
+        5.44725003058768775090E0,
+        8.74716500199817011941E0,
+        5.30324038235394892183E0,
+        9.99999999999999997821E-1,
+        0.0
+    };
+
+    const double PQ[8] = {
+        9.24408810558863637013E-4,
+        8.56288474354474431428E-2,
+        1.25352743901058953537E0,
+        5.47097740330417105182E0,
+        8.76190883237069594232E0,
+        5.30605288235394617618E0,
+        1.00000000000000000218E0,
+        0.0
+    };
+
+    const double QP[8] = {
+        -1.13663838898469149931E-2,
+        -1.28252718670509318512E0,
+        -1.95539544257735972385E1,
+        -9.32060152123768231369E1,
+        -1.77681167980488050595E2,
+        -1.47077505154951170175E2,
+        -5.14105326766599330220E1,
+        -6.05014350600728481186E0,
+    };
+
+    const double QQ[8] = {
+        /*  1.00000000000000000000E0,*/
+        6.43178256118178023184E1,
+        8.56430025976980587198E2,
+        3.88240183605401609683E3,
+        7.24046774195652478189E3,
+        5.93072701187316984827E3,
+        2.06209331660327847417E3,
+        2.42005740240291393179E2,
+    };
+
+    const double YP[8] = {
+        1.55924367855235737965E4,
+        -1.46639295903971606143E7,
+        5.43526477051876500413E9,
+        -9.82136065717911466409E11,
+        8.75906394395366999549E13,
+        -3.46628303384729719441E15,
+        4.42733268572569800351E16,
+        -1.84950800436986690637E16,
+    };
+
+    const double YQ[7] = {
+        /* 1.00000000000000000000E0,*/
+        1.04128353664259848412E3,
+        6.26107330137134956842E5,
+        2.68919633393814121987E8,
+        8.64002487103935000337E10,
+        2.02979612750105546709E13,
+        3.17157752842975028269E15,
+        2.50596256172653059228E17,
+    };
+
+    const double RP[8] = {
+        -4.79443220978201773821E9,
+        1.95617491946556577543E12,
+        -2.49248344360967716204E14,
+        9.70862251047306323952E15,
+        0.0,
+        0.0,
+        0.0,
+        0.0
+    };
+
+    const double RQ[8] = {
+        /* 1.00000000000000000000E0,*/
+        4.99563147152651017219E2,
+        1.73785401676374683123E5,
+        4.84409658339962045305E7,
+        1.11855537045356834862E10,
+        2.11277520115489217587E12,
+        3.10518229857422583814E14,
+        3.18121955943204943306E16,
+        1.71086294081043136091E18,
+    };
+
+
+    /*  5.783185962946784521175995758455807035071 */
+    const double DR1 = 5.78318596294678452118E0;
+    /* 30.47126234366208639907816317502275584842 */
+    const double DR2 = 3.04712623436620863991E1;
+
+
+    if( x < 0 )
+	    x = -x;
+
+    if( x <= 5.0 ) {
+	    z = x * x;
+	    if( x < 1.0e-5 )
+		    return( 1.0 - z/4.0 );
+
+	    p = (z - DR1) * (z - DR2);
+	    p = p * polevl( z, RP, 3)/p1evl( z, RQ, 8 );
+	    return( p );
+	}
+
+    w = 5.0/x;
+    q = 25.0/(x*x);
+    p = polevl( q, PP, 6)/polevl( q, PQ, 6 );
+    q = polevl( q, QP, 7)/p1evl( q, QQ, 7 );
+    xn = x - PIO4;
+
+    double sn, cn;
+    SINCOS(xn, sn, cn);
+    p = p * cn - w * q * sn;
+
+    return( p * SQ2OPI / sqrt(x) );
+}
+
+
Index: sasmodels/models/lib/j1d.c
===================================================================
--- sasmodels/models/lib/j1d.c	(revision 3936ad3134459cadbcc9696bac14e0ffc9da869b)
+++ sasmodels/models/lib/j1d.c	(revision 3936ad3134459cadbcc9696bac14e0ffc9da869b)
@@ -0,0 +1,184 @@
+/*							j1.c
+ *
+ *	Bessel function of order one
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, j1();
+ *
+ * y = j1( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns Bessel function of order one of the argument.
+ *
+ * The domain is divided into the intervals [0, 8] and
+ * (8, infinity). In the first interval a 24 term Chebyshev
+ * expansion is used. In the second, the asymptotic
+ * trigonometric representation is employed using two
+ * rational functions of degree 5/5.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Absolute error:
+ * arithmetic   domain      # trials      peak         rms
+ *    DEC       0, 30       10000       4.0e-17     1.1e-17
+ *    IEEE      0, 30       30000       2.6e-16     1.1e-16
+ *
+ *
+ */
+/*							y1.c
+ *
+ *	Bessel function of second kind of order one
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, y1();
+ *
+ * y = y1( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns Bessel function of the second kind of order one
+ * of the argument.
+ *
+ * The domain is divided into the intervals [0, 8] and
+ * (8, infinity). In the first interval a 25 term Chebyshev
+ * expansion is used, and a call to j1() is required.
+ * In the second, the asymptotic trigonometric representation
+ * is employed using two rational functions of degree 5/5.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Absolute error:
+ * arithmetic   domain      # trials      peak         rms
+ *    DEC       0, 30       10000       8.6e-17     1.3e-17
+ *    IEEE      0, 30       30000       1.0e-15     1.3e-16
+ *
+ * (error criterion relative when |y1| > 1).
+ *
+ */
+
+/*
+Cephes Math Library Release 2.8:  June, 2000
+Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
+*/
+double j1(double );
+
+double j1(double x) {
+
+    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,
+    };
+
+    //const double Z1 = 1.46819706421238932572E1;
+    //const double Z2 = 4.92184563216946036703E1;
+
+    w = x;
+    if( x < 0 )
+	    w = -x;
+
+    if( w <= 5.0 )
+	{
+	    z = x * x;
+	    w = polevl( z, RP, 3 ) / p1evl( z, RQ, 8 );
+	    w = w * x * (z - Z1) * (z - Z2);
+	    return( w );
+	}
+
+    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 );
+    xn = x - THPIO4;
+
+    double sn, cn;
+    SINCOS(xn, sn, cn);
+    p = p * cn - w * q * sn;
+
+    return( p * SQ2OPI / sqrt(x) );
+}
Index: sasmodels/models/lib/jnd.c
===================================================================
--- sasmodels/models/lib/jnd.c	(revision 3936ad3134459cadbcc9696bac14e0ffc9da869b)
+++ sasmodels/models/lib/jnd.c	(revision 3936ad3134459cadbcc9696bac14e0ffc9da869b)
@@ -0,0 +1,121 @@
+/*							jn.c
+ *
+ *	Bessel function of integer order
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * int n;
+ * double x, y, jn();
+ *
+ * y = jn( n, x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns Bessel function of order n, where n is a
+ * (possibly negative) integer.
+ *
+ * The ratio of jn(x) to j0(x) is computed by backward
+ * recurrence.  First the ratio jn/jn-1 is found by a
+ * continued fraction expansion.  Then the recurrence
+ * relating successive orders is applied until j0 or j1 is
+ * reached.
+ *
+ * If n = 0 or 1 the routine for j0 or j1 is called
+ * directly.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Absolute error:
+ * arithmetic   range      # trials      peak         rms
+ *    DEC       0, 30        5500       6.9e-17     9.3e-18
+ *    IEEE      0, 30        5000       4.4e-16     7.9e-17
+ *
+ *
+ * Not suitable for large n or x. Use jv() instead.
+ *
+ */
+
+
+/*							jn.c
+Cephes Math Library Release 2.8:  June, 2000
+Copyright 1984, 1987, 2000 by Stephen L. Moshier
+*/
+
+double jn( int n, double x );
+#define MACHEP 1.11022302462515654042E-16
+
+double jn( int n, double x )
+{
+double pkm2, pkm1, pk, xk, r, ans;
+int k, sign;
+
+if( n < 0 )
+	{
+	n = -n;
+	if( (n & 1) == 0 )	/* -1**n */
+		sign = 1;
+	else
+		sign = -1;
+	}
+else
+	sign = 1;
+
+if( x < 0.0 )
+	{
+	if( n & 1 )
+		sign = -sign;
+	x = -x;
+	}
+
+if( n == 0 )
+	return( sign * j0(x) );
+if( n == 1 )
+	return( sign * j1(x) );
+if( n == 2 )
+	return( sign * (2.0 * j1(x) / x  -  j0(x)) );
+
+if( x < MACHEP )
+	return( 0.0 );
+
+
+pk = 2 * (n + k);
+ans = pk;
+xk = x * x;
+
+do
+	{
+	pk -= 2.0;
+	ans = pk - (xk/ans);
+	}
+while( --k > 0 );
+ans = x/ans;
+
+/* backward recurrence */
+
+pk = 1.0;
+pkm1 = 1.0/ans;
+k = n-1;
+r = 2 * k;
+
+do
+	{
+	pkm2 = (pkm1 * r  -  pk * x) / x;
+	pk = pkm1;
+	pkm1 = pkm2;
+	r -= 2.0;
+	}
+while( --k > 0 );
+
+if( fabs(pk) > fabs(pkm1) )
+	ans = j1(x)/pk;
+else
+	ans = j0(x)/pkm1;
+return( sign * ans );
+}
+
Index: sasmodels/models/lib/polevl.c
===================================================================
--- sasmodels/models/lib/polevl.c	(revision 3936ad3134459cadbcc9696bac14e0ffc9da869b)
+++ sasmodels/models/lib/polevl.c	(revision 3936ad3134459cadbcc9696bac14e0ffc9da869b)
@@ -0,0 +1,86 @@
+/*							polevl.c
+ *							p1evl.c
+ *
+ *	Evaluate polynomial
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * int N;
+ * double x, y, coef[N+1], polevl[];
+ *
+ * y = polevl( x, coef, N );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Evaluates polynomial of degree N:
+ *
+ *                     2          N
+ * y  =  C  + C x + C x  +...+ C x
+ *        0    1     2          N
+ *
+ * Coefficients are stored in reverse order:
+ *
+ * coef[0] = C  , ..., coef[N] = C  .
+ *            N                   0
+ *
+ *  The function p1evl() assumes that coef[N] = 1.0 and is
+ * omitted from the array.  Its calling arguments are
+ * otherwise the same as polevl().
+ *
+ *
+ * SPEED:
+ *
+ * In the interest of speed, there are no checks for out
+ * of bounds arithmetic.  This routine is used by most of
+ * the functions in the library.  Depending on available
+ * equipment features, the user may wish to rewrite the
+ * program in microcode or assembly language.
+ *
+ */
+
+
+
+/*
+Cephes Math Library Release 2.1:  December, 1988
+Copyright 1984, 1987, 1988 by Stephen L. Moshier
+Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+*/
+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 ;
+
+}
+
+/*							p1evl()	*/
+/*                                          N
+ * Evaluate polynomial when coefficient of x  is 1.0.
+ * Otherwise same as polevl.
+ */
+
+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 );
+
+}
