Index: sasmodels/models/lib/j0d.c
===================================================================
--- sasmodels/models/lib/j0d.c	(revision 3936ad3134459cadbcc9696bac14e0ffc9da869b)
+++ sasmodels/models/lib/j0d.c	(revision 0a9d21936ceb8c9316bfd9de49f989fecffc9413)
@@ -96,6 +96,8 @@
 double j0( double );
 
-double j0(double x)
-{
+double j0(double x) {
+
+//Cephes single precission
+if ( FLOAT_SIZE>4 ) {
     double w, z, p, q, xn;
 
@@ -103,4 +105,7 @@
     const double SQ2OPI = 7.9788456080286535587989E-1;
     const double PIO4 = 7.85398163397448309616E-1;
+
+    const double DR1 = 5.78318596294678452118E0;
+    const double DR2 = 3.04712623436620863991E1;
 
     const double PP[8] = {
@@ -193,11 +198,4 @@
     };
 
-
-    /*  5.783185962946784521175995758455807035071 */
-    const double DR1 = 5.78318596294678452118E0;
-    /* 30.47126234366208639907816317502275584842 */
-    const double DR2 = 3.04712623436620863991E1;
-
-
     if( x < 0 )
 	    x = -x;
@@ -225,4 +223,82 @@
     return( p * SQ2OPI / sqrt(x) );
 }
-
-
+//Cephes single precission
+else {
+    double xx, w, z, p, q, xn;
+
+    const double YZ1 =  0.43221455686510834878;
+    const double YZ2 = 22.401876406482861405;
+    const double YZ3 = 64.130620282338755553;
+    const double DR1 =  5.78318596294678452118;
+    const double PIO4F = 0.7853981633974483096;
+
+    double MO[8] = {
+        -6.838999669318810E-002,
+        1.864949361379502E-001,
+        -2.145007480346739E-001,
+        1.197549369473540E-001,
+        -3.560281861530129E-003,
+        -4.969382655296620E-002,
+        -3.355424622293709E-006,
+        7.978845717621440E-001
+    };
+
+    double PH[8] = {
+        3.242077816988247E+001,
+        -3.630592630518434E+001,
+        1.756221482109099E+001,
+        -4.974978466280903E+000,
+        1.001973420681837E+000,
+        -1.939906941791308E-001,
+        6.490598792654666E-002,
+        -1.249992184872738E-001
+    };
+
+    double YP[8] = {
+        9.454583683980369E-008,
+        -9.413212653797057E-006,
+        5.344486707214273E-004,
+        -1.584289289821316E-002,
+        1.707584643733568E-001,
+        0.0,
+        0.0,
+        0.0
+    };
+
+    double JP[8] = {
+        -6.068350350393235E-008,
+        6.388945720783375E-006,
+        -3.969646342510940E-004,
+        1.332913422519003E-002,
+        -1.729150680240724E-001,
+        0.0,
+        0.0,
+        0.0
+    };
+
+    if( x < 0 )
+	    xx = -x;
+    else
+	    xx = x;
+
+    if( x <= 2.0 ) {
+	    z = xx * xx;
+	    if( x < 1.0e-3 )
+		    return( 1.0 - 0.25*z );
+
+	    p = (z-DR1) * polevl( z, JP, 4);
+	    return( p );
+	}
+
+    q = 1.0/x;
+    w = sqrtf(q);
+
+    p = w * polevl( q, MO, 7);
+    w = q*q;
+    xn = q * polevl( w, PH, 7) - PIO4F;
+    p = p * cosf(xn + xx);
+    return(p);
+}
+}
+
+
Index: sasmodels/models/lib/j1d.c
===================================================================
--- sasmodels/models/lib/j1d.c	(revision 3936ad3134459cadbcc9696bac14e0ffc9da869b)
+++ sasmodels/models/lib/j1d.c	(revision 0a9d21936ceb8c9316bfd9de49f989fecffc9413)
@@ -79,4 +79,7 @@
 
 double j1(double x) {
+
+//Cephes double pression function
+if (FLOAT_SIZE>4) {
 
     double w, z, p, q, xn;
@@ -156,7 +159,4 @@
     };
 
-    //const double Z1 = 1.46819706421238932572E1;
-    //const double Z2 = 4.92184563216946036703E1;
-
     w = x;
     if( x < 0 )
@@ -182,3 +182,68 @@
 
     return( p * SQ2OPI / sqrt(x) );
+
 }
+//Single precission version of cephes
+else {
+    double xx, w, z, p, q, xn;
+
+    const double Z1 = 1.46819706421238932572E1;
+    const double THPIO4F =  2.35619449019234492885;    /* 3*pi/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 )
+	    xx = -x;
+
+    if( xx <= 2.0 )
+	{
+	    z = xx * xx;
+	    p = (z-Z1) * xx * polevl( z, JP, 4 );
+	    return( p );
+	}
+
+    q = 1.0/x;
+    w = sqrt(q);
+
+    p = w * polevl( q, MO1, 7);
+    w = q*q;
+    xn = q * polevl( w, PH1, 7) - THPIO4F;
+    p = p * cos(xn + xx);
+
+    return(p);
+}
+}
+
Index: sasmodels/models/lib/jnd.c
===================================================================
--- sasmodels/models/lib/jnd.c	(revision 3936ad3134459cadbcc9696bac14e0ffc9da869b)
+++ sasmodels/models/lib/jnd.c	(revision 0a9d21936ceb8c9316bfd9de49f989fecffc9413)
@@ -49,73 +49,103 @@
 
 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;
+double jn( int n, double x ) {
 
-if( n < 0 )
-	{
-	n = -n;
-	if( (n & 1) == 0 )	/* -1**n */
-		sign = 1;
-	else
-		sign = -1;
+    const double MACHEP = 1.11022302462515654042E-16;
+    double pkm2, pkm1, pk, xk, r, ans, xinv;
+    int k, sign;
+
+    if( n < 0 ) {
+	    n = -n;
+	    if( (n & 1) == 0 )	/* -1**n */
+		    sign = 1;
+	    else
+		    sign = -1;
 	}
-else
-	sign = 1;
+    else
+	    sign = 1;
 
-if( x < 0.0 )
-	{
-	if( n & 1 )
-		sign = -sign;
-	x = -x;
+    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( 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 );
+    if( x < MACHEP )
+	    return( 0.0 );
 
+    if (FLOAT_SIZE > 4)
+        k = 53;
+    else
+        k = 24;
 
-pk = 2 * (n + k);
-ans = pk;
-xk = x * x;
+    pk = 2 * (n + k);
+    ans = pk;
+    xk = x * x;
 
-do
-	{
-	pk -= 2.0;
-	ans = pk - (xk/ans);
+    do {
+	    pk -= 2.0;
+	    ans = pk - (xk/ans);
+	} while( --k > 0 );
+
+    /* backward recurrence */
+
+    pk = 1.0;
+
+    if (FLOAT_SIZE > 4) {
+        ans = x/ans;
+        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 );
 	}
-while( --k > 0 );
-ans = x/ans;
+    else {
+        xinv = 1.0/x;
+        pkm1 = ans * xinv;
+        k = n-1;
+        r = (float )(2 * k);
 
-/* backward recurrence */
+        do {
+	        pkm2 = (pkm1 * r  -  pk * x) * xinv;
+	        pk = pkm1;
+	        pkm1 = pkm2;
+	        r -= 2.0;
+	    }
+        while( --k > 0 );
 
-pk = 1.0;
-pkm1 = 1.0/ans;
-k = n-1;
-r = 2 * k;
+        r = pk;
+        if( r < 0 )
+	        r = -r;
+        ans = pkm1;
+        if( ans < 0 )
+	        ans = -ans;
 
-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 );
+        if( r > ans )  /* if( fabs(pk) > fabs(pkm1) ) */
+	        ans = sign * j1(x)/pk;
+        else
+	        ans = sign * j0(x)/pkm1;
+        return( ans );
+    }
 }
 
