Index: sasmodels/models/lib/sas_J0.c
===================================================================
--- sasmodels/models/lib/sas_J0.c	(revision 0278e3fec71d4a97458c09903fd5e3bd354166c1)
+++ sasmodels/models/lib/sas_J0.c	(revision 3f8584a2835abe8ddfd9dfbdb33f0edb8715277c)
@@ -53,4 +53,6 @@
  * except YP, YQ which are designed for absolute error. */
 
+#if FLOAT_SIZE>4
+//Cephes double precission
 
  constant double PPJ0[8] = {
@@ -144,5 +146,45 @@
   };
 
- constant double MOJ0[8] = {
+double j0(double x);
+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 DR1 = 5.78318596294678452118E0;
+    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, RPJ0, 3)/p1evl( z, RQJ0, 8 );
+        return( p );
+    }
+
+    w = 5.0/x;
+    q = 25.0/(x*x);
+    p = polevl( q, PPJ0, 6)/polevl( q, PQJ0, 6 );
+    q = polevl( q, QPJ0, 7)/p1evl( q, QQJ0, 7 );
+    xn = x - PIO4;
+
+    double sn, cn;
+    SINCOS(xn, sn, cn);
+    p = p * cn - w * q * sn;
+
+    return( p * SQ2OPI / sqrt(x) );
+}
+#else
+
+ constant float MOJ0[8] = {
         -6.838999669318810E-002,
         1.864949361379502E-001,
@@ -155,5 +197,5 @@
   };
 
- constant double PHJ0[8] = {
+ constant float PHJ0[8] = {
         3.242077816988247E+001,
         -3.630592630518434E+001,
@@ -166,5 +208,5 @@
   };
 
- constant double JPJ0[8] = {
+ constant float JPJ0[8] = {
         -6.068350350393235E-008,
         6.388945720783375E-006,
@@ -177,67 +219,29 @@
  };
 
-double sas_J0(double x);
-double sas_J0(double x)
+//Cephes single precission
+float j0f(float x);
+float j0f(float x)
 {
-
-//Cephes single precission
-#if FLOAT_SIZE>4
-    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 DR1 = 5.78318596294678452118E0;
-    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, RPJ0, 3)/p1evl( z, RQJ0, 8 );
-	    return( p );
-	}
-
-    w = 5.0/x;
-    q = 25.0/(x*x);
-    p = polevl( q, PPJ0, 6)/polevl( q, PQJ0, 6 );
-    q = polevl( q, QPJ0, 7)/p1evl( q, QQJ0, 7 );
-    xn = x - PIO4;
-
-    double sn, cn;
-    SINCOS(xn, sn, cn);
-    p = p * cn - w * q * sn;
-
-    return( p * SQ2OPI / sqrt(x) );
-//Cephes single precission
-#else
-    double xx, w, z, p, q, xn;
+    float 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;
+    const float DR1 =  5.78318596294678452118;
+    const float PIO4F = 0.7853981633974483096;
 
     if( x < 0 )
-	    xx = -x;
+        xx = -x;
     else
-	    xx = x;
+        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, JPJ0, 4);
-	    return( p );
-	}
+        z = xx * xx;
+        if( x < 1.0e-3 )
+            return( 1.0 - 0.25*z );
+
+        p = (z-DR1) * polevl( z, JPJ0, 4);
+        return( p );
+    }
 
     q = 1.0/x;
@@ -249,7 +253,10 @@
     p = p * cos(xn + xx);
     return(p);
+}
 #endif
 
-}
-
-
+#if FLOAT_SIZE>4
+#define sas_J0 j0
+#else
+#define sas_J0 j0f
+#endif
Index: sasmodels/models/lib/sas_J1.c
===================================================================
--- sasmodels/models/lib/sas_J1.c	(revision 0278e3fec71d4a97458c09903fd5e3bd354166c1)
+++ sasmodels/models/lib/sas_J1.c	(revision 3f8584a2835abe8ddfd9dfbdb33f0edb8715277c)
@@ -40,4 +40,7 @@
 */
 
+#if FLOAT_SIZE>4
+
+//Cephes double pression function
 constant double RPJ1[8] = {
     -8.99971225705559398224E8,
@@ -102,5 +105,44 @@
     0.0 };
 
-constant double JPJ1[8] = {
+double j1(double x);
+double j1(double x)
+{
+
+    double w, z, p, q, xn;
+
+    const double Z1 = 1.46819706421238932572E1;
+    const double Z2 = 4.92184563216946036703E1;
+    const double THPIO4 =  2.35619449019234492885;
+    const double SQ2OPI = 0.79788456080286535588;
+
+    w = x;
+    if( x < 0 )
+        w = -x;
+
+    if( w <= 5.0 ) {
+        z = x * x;
+        w = polevl( z, RPJ1, 3 ) / p1evl( z, RQJ1, 8 );
+        w = w * x * (z - Z1) * (z - Z2);
+        return( w );
+    }
+
+    w = 5.0/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;
+
+    double sn, cn;
+    SINCOS(xn, sn, cn);
+    p = p * cn - w * q * sn;
+
+    return( p * SQ2OPI / sqrt(x) );
+}
+
+#else
+//Single precission version of cephes
+constant float JPJ1[8] = {
     -4.878788132172128E-009,
     6.009061827883699E-007,
@@ -113,5 +155,5 @@
     };
 
-constant double MO1J1[8] = {
+constant float MO1J1[8] = {
     6.913942741265801E-002,
     -2.284801500053359E-001,
@@ -124,5 +166,5 @@
     };
 
-constant double PH1J1[8] = {
+constant float PH1J1[8] = {
     -4.497014141919556E+001,
     5.073465654089319E+001,
@@ -135,63 +177,23 @@
     };
 
-double sas_J1(double x);
-double sas_J1(double x)
+float j1f(float x);
+float j1f(float x)
 {
 
-//Cephes double pression function
-#if FLOAT_SIZE>4
-
-    double w, z, p, q, xn;
-
-    const double Z1 = 1.46819706421238932572E1;
-    const double Z2 = 4.92184563216946036703E1;
-    const double THPIO4 =  2.35619449019234492885;
-    const double SQ2OPI = 0.79788456080286535588;
-
-    w = x;
-    if( x < 0 )
-	    w = -x;
-
-    if( w <= 5.0 )
-	{
-	    z = x * x;
-	    w = polevl( z, RPJ1, 3 ) / p1evl( z, RQJ1, 8 );
-	    w = w * x * (z - Z1) * (z - Z2);
-	    return( w );
-	}
-
-    w = 5.0/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;
-
-    double sn, cn;
-    SINCOS(xn, sn, cn);
-    p = p * cn - w * q * sn;
-
-    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 */
+    float xx, w, z, p, q, xn;
+
+    const float Z1 = 1.46819706421238932572E1;
+    const float THPIO4F =  2.35619449019234492885;    /* 3*pi/4 */
 
 
     xx = x;
     if( xx < 0 )
-	    xx = -x;
-
-    if( xx <= 2.0 )
-	{
-	    z = xx * xx;
-	    p = (z-Z1) * xx * polevl( z, JPJ1, 4 );
-	    return( p );
-	}
+        xx = -x;
+
+    if( xx <= 2.0 ) {
+        z = xx * xx;
+        p = (z-Z1) * xx * polevl( z, JPJ1, 4 );
+        return( p );
+    }
 
     q = 1.0/x;
@@ -204,6 +206,12 @@
 
     return(p);
+}
 #endif
-}
+
+#if FLOAT_SIZE>4
+#define sas_J1 j1
+#else
+#define sas_J1 j1f
+#endif
 
 //Finally J1c function that equals 2*J1(x)/x
Index: sasmodels/models/lib/sas_JN.c
===================================================================
--- sasmodels/models/lib/sas_JN.c	(revision e6408d049d3598925e345d8adec1cdb0bb0c6adb)
+++ sasmodels/models/lib/sas_JN.c	(revision 3f8584a2835abe8ddfd9dfbdb33f0edb8715277c)
@@ -48,8 +48,10 @@
 */
 
-double sas_JN( int n, double x );
-
-double sas_JN( int n, double x ) {
-
+#if FLOAT_SIZE > 4
+
+double jn( int n, double x );
+double jn( int n, double x ) {
+
+    // PAK: seems to be machine epsilon/2
     const double MACHEP = 1.11022302462515654042E-16;
     double pkm2, pkm1, pk, xk, r, ans;
@@ -57,35 +59,30 @@
 
     if( n < 0 ) {
-	    n = -n;
-	    if( (n & 1) == 0 )	/* -1**n */
-		    sign = 1;
-	    else
-		    sign = -1;
-	}
-    else
-	    sign = 1;
+        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 & 1 )
+            sign = -sign;
+    x = -x;
+    }
 
     if( n == 0 )
-	    return( sign * sas_J0(x) );
+        return( sign * j0(x) );
     if( n == 1 )
-	    return( sign * sas_J1(x) );
+        return( sign * j1(x) );
     if( n == 2 )
-	    return( sign * (2.0 * sas_J1(x) / x  -  sas_J0(x)) );
+        return( sign * (2.0 * j1(x) / x  -  j0(x)) );
 
     if( x < MACHEP )
-	    return( 0.0 );
-
-    #if FLOAT_SIZE > 4
-        k = 53;
-    #else
-        k = 24;
-    #endif
-
+        return( 0.0 );
+
+    k = 53;
     pk = 2 * (n + k);
     ans = pk;
@@ -93,7 +90,7 @@
 
     do {
-	    pk -= 2.0;
-	    ans = pk - (xk/ans);
-	} while( --k > 0 );
+        pk -= 2.0;
+        ans = pk - (xk/ans);
+    } while( --k > 0 );
 
     /* backward recurrence */
@@ -101,52 +98,104 @@
     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 = sas_J1(x)/pk;
+    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 );
+}
+
+#else
+
+float jnf(int n, float x)
+{
+    // PAK: seems to be machine epsilon/2
+    const double MACHEP = 5.9604645e-08;
+    float pkm2, pkm1, pk, xk, r, ans;
+    int k, sign;
+
+    if( n < 0 ) {
+        n = -n;
+        if( (n & 1) == 0 ) /* -1**n */
+            sign = 1;
         else
-	        ans = sas_J0(x)/pkm1;
-
-	    return( sign * ans );
-
-    #else
-        const double xinv = 1.0/x;
-        pkm1 = ans * xinv;
-        k = n-1;
-        r = (float )(2 * k);
-
-        do {
-	        pkm2 = (pkm1 * r  -  pk * x) * xinv;
-	        pk = pkm1;
-	        pkm1 = pkm2;
-	        r -= 2.0;
-	    }
-        while( --k > 0 );
-
-        r = pk;
-        if( r < 0 )
-	        r = -r;
-        ans = pkm1;
-        if( ans < 0 )
-	        ans = -ans;
-
-        if( r > ans )  /* if( fabs(pk) > fabs(pkm1) ) */
-	        ans = sign * sas_J1(x)/pk;
-        else
-	        ans = sign * sas_J0(x)/pkm1;
-        return( ans );
-    #endif
+            sign = -1;
+    } else {
+        sign = 1;
+    }
+
+    if( x < 0.0 ) {
+        if( n & 1 )
+            sign = -sign;
+        x = -x;
+    }
+
+    if( n == 0 )
+        return( sign * j0f(x) );
+    if( n == 1 )
+        return( sign * j1f(x) );
+    if( n == 2 )
+        return( sign * (2.0 * j1f(x) / x  -  j0f(x)) );
+
+    if( x < MACHEP )
+        return( 0.0 );
+
+    k = 24;
+    pk = 2 * (n + k);
+    ans = pk;
+    xk = x * x;
+
+    do {
+        pk -= 2.0;
+        ans = pk - (xk/ans);
+    } while( --k > 0 );
+
+    /* backward recurrence */
+
+    pk = 1.0;
+
+    const float xinv = 1.0/x;
+    pkm1 = ans * xinv;
+    k = n-1;
+    r = (float )(2 * k);
+
+    do {
+        pkm2 = (pkm1 * r  -  pk * x) * xinv;
+        pk = pkm1;
+        pkm1 = pkm2;
+        r -= 2.0;
+    } while( --k > 0 );
+
+    r = pk;
+    if( r < 0 )
+        r = -r;
+    ans = pkm1;
+    if( ans < 0 )
+        ans = -ans;
+
+    if( r > ans )  /* if( fabs(pk) > fabs(pkm1) ) */
+        ans = sign * j1f(x)/pk;
+    else
+        ans = sign * j0f(x)/pkm1;
+    return( ans );
 }
-
+#endif
+
+#if FLOAT_SIZE>4
+#define sas_JN jn
+#else
+#define sas_JN jnf
+#endif
