Index: sasmodels/models/lib/j0_cephes.c
===================================================================
--- sasmodels/models/lib/j0_cephes.c	(revision bfef528b72ff0998652ebfa1704aa03b745dc7af)
+++ sasmodels/models/lib/j0_cephes.c	(revision bfef528b72ff0998652ebfa1704aa03b745dc7af)
@@ -0,0 +1,304 @@
+/*							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) {
+
+//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;
+
+    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,
+    };
+
+    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) );
+//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);
+#endif
+
+}
+
+
Index: smodels/models/lib/j0d.c
===================================================================
--- sasmodels/models/lib/j0d.c	(revision 07142f3db170174ebf3c5f6b28931ae502435967)
+++ 	(revision )
@@ -1,304 +1,0 @@
-/*							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) {
-
-//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;
-
-    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,
-    };
-
-    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) );
-//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);
-#endif
-
-}
-
-
Index: sasmodels/models/lib/j1_cephes.c
===================================================================
--- sasmodels/models/lib/j1_cephes.c	(revision bfef528b72ff0998652ebfa1704aa03b745dc7af)
+++ sasmodels/models/lib/j1_cephes.c	(revision bfef528b72ff0998652ebfa1704aa03b745dc7af)
@@ -0,0 +1,249 @@
+/*							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) {
+
+//Cephes double pression function
+#if FLOAT_SIZE>4
+
+    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;
+    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) );
+
+
+//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);
+#endif
+}
+
Index: smodels/models/lib/j1d.c
===================================================================
--- sasmodels/models/lib/j1d.c	(revision 07142f3db170174ebf3c5f6b28931ae502435967)
+++ 	(revision )
@@ -1,249 +1,0 @@
-/*							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) {
-
-//Cephes double pression function
-#if FLOAT_SIZE>4
-
-    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;
-    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) );
-
-
-//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);
-#endif
-}
-
Index: sasmodels/models/lib/jn_cephes.c
===================================================================
--- sasmodels/models/lib/jn_cephes.c	(revision bfef528b72ff0998652ebfa1704aa03b745dc7af)
+++ sasmodels/models/lib/jn_cephes.c	(revision bfef528b72ff0998652ebfa1704aa03b745dc7af)
@@ -0,0 +1,152 @@
+/*							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 );
+
+double jn( int n, double x ) {
+
+    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;
+
+    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 );
+
+    #if FLOAT_SIZE > 4
+        k = 53;
+    #else
+        k = 24;
+    #endif
+
+    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;
+
+    #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 );
+
+    #else
+        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 * j1(x)/pk;
+        else
+	        ans = sign * j0(x)/pkm1;
+        return( ans );
+    #endif
+}
+
Index: smodels/models/lib/jnd.c
===================================================================
--- sasmodels/models/lib/jnd.c	(revision 07142f3db170174ebf3c5f6b28931ae502435967)
+++ 	(revision )
@@ -1,152 +1,0 @@
-/*							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 );
-
-double jn( int n, double x ) {
-
-    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;
-
-    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 );
-
-    #if FLOAT_SIZE > 4
-        k = 53;
-    #else
-        k = 24;
-    #endif
-
-    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;
-
-    #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 );
-
-    #else
-        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 * j1(x)/pk;
-        else
-	        ans = sign * j0(x)/pkm1;
-        return( ans );
-    #endif
-}
-
Index: smodels/models/lib/lgamma.c
===================================================================
--- sasmodels/models/lib/lgamma.c	(revision cae6ce63a3ff4a4d02996fecf8be2ac218727beb)
+++ 	(revision )
@@ -1,330 +1,0 @@
-/*							lgamma.c
- *
- *	Log Gamma function
- *
- */
-
-/*							lgamma()
- *
- *	Natural logarithm of gamma function
- *
- *
- *
- * SYNOPSIS:
- *
- * double x, y, lgamma();
- * extern int sgngam;
- *
- * y = lgamma( x );
- *
- *
- *
- * DESCRIPTION:
- *
- * Returns the base e (2.718...) logarithm of the absolute
- * value of the gamma function of the argument.
- * The sign (+1 or -1) of the gamma function is returned in a
- * global (extern) variable named sgngam.
- *
- * For arguments greater than 13, the logarithm of the gamma
- * function is approximated by the logarithmic version of
- * Stirling's formula using a polynomial approximation of
- * degree 4. Arguments between -33 and +33 are reduced by
- * recurrence to the interval [2,3] of a rational approximation.
- * The cosecant reflection formula is employed for arguments
- * less than -33.
- *
- * Arguments greater than MAXLGM return MAXNUM and an error
- * message.  MAXLGM = 2.035093e36 for DEC
- * arithmetic or 2.556348e305 for IEEE arithmetic.
- *
- *
- *
- * ACCURACY:
- *
- *
- * arithmetic      domain        # trials     peak         rms
- *    DEC     0, 3                  7000     5.2e-17     1.3e-17
- *    DEC     2.718, 2.035e36       5000     3.9e-17     9.9e-18
- *    IEEE    0, 3                 28000     5.4e-16     1.1e-16
- *    IEEE    2.718, 2.556e305     40000     3.5e-16     8.3e-17
- * The error criterion was relative when the function magnitude
- * was greater than one but absolute when it was less than one.
- *
- * The following test used the relative error criterion, though
- * at certain points the relative error could be much higher than
- * indicated.
- *    IEEE    -200, -4             10000     4.8e-16     1.3e-16
- *
- */
-
-
-/*
-Cephes Math Library Release 2.8:  June, 2000
-Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
-*/
-
-
-double lgamma( double );
-
-double lgamma( double x) {
-
-#if FLOAT_SIZE > 4
-    double p, q, u, w, z;
-    int i;
-    int sgngam = 1;
-
-    const double LS2PI  =  0.91893853320467274178;
-    const double MAXLGM = 2.556348e305;
-    const double MAXNUM =  1.79769313486231570815E308;
-    const double LOGPI = 1.14472988584940017414;
-    const double PI = M_PI;
-
-    double A[8] = {
-        8.11614167470508450300E-4,
-        -5.95061904284301438324E-4,
-        7.93650340457716943945E-4,
-        -2.77777777730099687205E-3,
-        8.33333333333331927722E-2,
-        0.0,
-        0.0,
-        0.0
-    };
-    double B[8] = {
-        -1.37825152569120859100E3,
-        -3.88016315134637840924E4,
-        -3.31612992738871184744E5,
-        -1.16237097492762307383E6,
-        -1.72173700820839662146E6,
-        -8.53555664245765465627E5,
-        0.0,
-        0.0
-    };
-
-    double C[8] = {
-        /* 1.00000000000000000000E0, */
-        -3.51815701436523470549E2,
-        -1.70642106651881159223E4,
-        -2.20528590553854454839E5,
-        -1.13933444367982507207E6,
-        -2.53252307177582951285E6,
-        -2.01889141433532773231E6,
-        0.0,
-        0.0
-    };
-
-    sgngam = 1;
-
-    if( x < -34.0 ) {
-	    q = -x;
-	    w = lanczos_gamma(q); /* note this modifies sgngam! */
-	    p = floor(q);
-	    if( p == q ) {
-            lgsing:
-        		goto loverf;
-		}
-	    i = p;
-	    if( (i & 1) == 0 )
-		    sgngam = -1;
-	    else
-		    sgngam = 1;
-	    z = q - p;
-	    if( z > 0.5 ) {
-		    p += 1.0;
-		    z = p - q;
-		}
-	    z = q * sin( PI * z );
-	    if( z == 0.0 )
-		    goto lgsing;
-	    z = LOGPI - log( z ) - w;
-	    return( z );
-	}
-
-    if( x < 13.0 ) {
-	    z = 1.0;
-	    p = 0.0;
-	    u = x;
-	    while( u >= 3.0 ) {
-		    p -= 1.0;
-		    u = x + p;
-		    z *= u;
-		}
-	    while( u < 2.0 ) {
-		    if( u == 0.0 )
-			    goto lgsing;
-		    z /= u;
-		    p += 1.0;
-		    u = x + p;
-		}
-	    if( z < 0.0 ) {
-		    sgngam = -1;
-		    z = -z;
-		}
-	    else
-		    sgngam = 1;
-	    if( u == 2.0 )
-		    return( log(z) );
-	    p -= 2.0;
-	    x = x + p;
-	    p = x * polevl( x, B, 5 ) / p1evl( x, C, 6);
-	    return( log(z) + p );
-	}
-
-    if( x > MAXLGM ) {
-        loverf:
-	        return( sgngam * MAXNUM );
-	}
-
-    q = ( x - 0.5 ) * log(x) - x + LS2PI;
-    if( x > 1.0e8 )
-	    return( q );
-
-    p = 1.0/(x*x);
-    if( x >= 1000.0 )
-	    q += ((   7.9365079365079365079365e-4 * p
-		    - 2.7777777777777777777778e-3) *p
-		    + 0.0833333333333333333333) / x;
-    else
-	    q += polevl( p, A, 4 ) / x;
-    return( q );
-#else
-
-    double p, q, w, z, xx;
-    double nx, tx;
-    int i, direction;
-    int sgngamf = 1;
-
-     double B[8] = {
-        6.055172732649237E-004,
-        -1.311620815545743E-003,
-        2.863437556468661E-003,
-        -7.366775108654962E-003,
-        2.058355474821512E-002,
-        -6.735323259371034E-002,
-        3.224669577325661E-001,
-        4.227843421859038E-001
-    };
-
-    /* log gamma(x+1), -.25 < x < .25 */
-    double C[8] = {
-        1.369488127325832E-001,
-        -1.590086327657347E-001,
-        1.692415923504637E-001,
-        -2.067882815621965E-001,
-        2.705806208275915E-001,
-        -4.006931650563372E-001,
-        8.224670749082976E-001,
-        -5.772156501719101E-001
-    };
-
-    /* log( sqrt( 2*pi ) ) */
-    const double LS2PI  =  0.91893853320467274178;
-    const double MAXLGM = 2.035093e36;
-    const double PIINV =  0.318309886183790671538;
-    const double MAXNUMF = 3.4028234663852885981170418348451692544e38;
-    const double PIF = 3.141592653589793238;
-    /* Logarithm of gamma function */
-
-    sgngamf = 1;
-
-    xx = x;
-    if( xx < 0.0 ) {
-	    q = -xx;
-	    w = lgamma(q); /* note this modifies sgngam! */
-	    p = floor(q);
-	    if( p == q )
-		    goto loverf;
-	    i = p;
-
-	    if( (i & 1) == 0 )
-		    sgngamf = -1;
-	    else
-		    sgngamf = 1;
-
-	    z = q - p;
-
-	    if( z > 0.5 ) {
-		    p += 1.0;
-		    z = p - q;
-	    }
-	    z = q * sin( PIF * z );
-	    if( z == 0.0 )
-		    goto loverf;
-	    z = -log( PIINV*z ) - w;
-	    return( z );
-	}
-
-    if( x < 6.5 ) {
-	    direction = 0;
-	    z = 1.0;
-	    tx = x;
-	    nx = 0.0;
-	    if( x >= 1.5 ) {
-		    while( tx > 2.5 ) {
-			    nx -= 1.0;
-			    tx = x + nx;
-			    z *=tx;
-			}
-		    x += nx - 2.0;
-            iv1r5:
-		    p = x * polevl( x, B, 7 );
-		    goto cont;
-	}
-	if( x >= 1.25 ) {
-		z *= x;
-		x -= 1.0; /* x + 1 - 2 */
-		direction = 1;
-		goto iv1r5;
-	}
-	if( x >= 0.75 ) {
-		x -= 1.0;
-		p = x * polevl( x, C, 7 );
-		q = 0.0;
-		goto contz;
-	}
-	while( tx < 1.5 ) {
-		if( tx == 0.0 )
-			goto loverf;
-		z *=tx;
-		nx += 1.0;
-		tx = x + nx;
-	}
-	direction = 1;
-	x += nx - 2.0;
-	p = x * polevl( x, B, 7 );
-
-    cont:
-	    if( z < 0.0 ) {
-		    sgngamf = -1;
-		    z = -z;
-        }
-	    else {
-		    sgngamf = 1;
-		}
-	    q = log(z);
-	    if( direction )
-		    q = -q;
-    contz:
-	    return( p + q );
-	}
-
-    if( x > MAXLGM ) {
-        loverf:
-	    return( sgngamf * MAXNUMF );
-	}
-
-    /* Note, though an asymptotic formula could be used for x >= 3,
-    * there is cancellation error in the following if x < 6.5.  */
-    q = LS2PI - x;
-    q += ( x - 0.5 ) * log(x);
-
-    if( x <= 1.0e4 ) {
-	    z = 1.0/x;
-	    p = z * z;
-	    q += ((    6.789774945028216E-004 * p
-		 - 2.769887652139868E-003 ) * p
-		+  8.333316229807355E-002 ) * z;
-	}
-    return( q );
-#endif
-}
Index: sasmodels/models/lib/lgamma_cephes.c
===================================================================
--- sasmodels/models/lib/lgamma_cephes.c	(revision bfef528b72ff0998652ebfa1704aa03b745dc7af)
+++ sasmodels/models/lib/lgamma_cephes.c	(revision bfef528b72ff0998652ebfa1704aa03b745dc7af)
@@ -0,0 +1,330 @@
+/*							lgamma.c
+ *
+ *	Log Gamma function
+ *
+ */
+
+/*							lgamma()
+ *
+ *	Natural logarithm of gamma function
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, lgamma();
+ * extern int sgngam;
+ *
+ * y = lgamma( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns the base e (2.718...) logarithm of the absolute
+ * value of the gamma function of the argument.
+ * The sign (+1 or -1) of the gamma function is returned in a
+ * global (extern) variable named sgngam.
+ *
+ * For arguments greater than 13, the logarithm of the gamma
+ * function is approximated by the logarithmic version of
+ * Stirling's formula using a polynomial approximation of
+ * degree 4. Arguments between -33 and +33 are reduced by
+ * recurrence to the interval [2,3] of a rational approximation.
+ * The cosecant reflection formula is employed for arguments
+ * less than -33.
+ *
+ * Arguments greater than MAXLGM return MAXNUM and an error
+ * message.  MAXLGM = 2.035093e36 for DEC
+ * arithmetic or 2.556348e305 for IEEE arithmetic.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *
+ * arithmetic      domain        # trials     peak         rms
+ *    DEC     0, 3                  7000     5.2e-17     1.3e-17
+ *    DEC     2.718, 2.035e36       5000     3.9e-17     9.9e-18
+ *    IEEE    0, 3                 28000     5.4e-16     1.1e-16
+ *    IEEE    2.718, 2.556e305     40000     3.5e-16     8.3e-17
+ * The error criterion was relative when the function magnitude
+ * was greater than one but absolute when it was less than one.
+ *
+ * The following test used the relative error criterion, though
+ * at certain points the relative error could be much higher than
+ * indicated.
+ *    IEEE    -200, -4             10000     4.8e-16     1.3e-16
+ *
+ */
+
+
+/*
+Cephes Math Library Release 2.8:  June, 2000
+Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
+*/
+
+
+double lgamma( double );
+
+double lgamma( double x) {
+
+#if FLOAT_SIZE > 4
+    double p, q, u, w, z;
+    int i;
+    int sgngam = 1;
+
+    const double LS2PI  =  0.91893853320467274178;
+    const double MAXLGM = 2.556348e305;
+    const double MAXNUM =  1.79769313486231570815E308;
+    const double LOGPI = 1.14472988584940017414;
+    const double PI = M_PI;
+
+    double A[8] = {
+        8.11614167470508450300E-4,
+        -5.95061904284301438324E-4,
+        7.93650340457716943945E-4,
+        -2.77777777730099687205E-3,
+        8.33333333333331927722E-2,
+        0.0,
+        0.0,
+        0.0
+    };
+    double B[8] = {
+        -1.37825152569120859100E3,
+        -3.88016315134637840924E4,
+        -3.31612992738871184744E5,
+        -1.16237097492762307383E6,
+        -1.72173700820839662146E6,
+        -8.53555664245765465627E5,
+        0.0,
+        0.0
+    };
+
+    double C[8] = {
+        /* 1.00000000000000000000E0, */
+        -3.51815701436523470549E2,
+        -1.70642106651881159223E4,
+        -2.20528590553854454839E5,
+        -1.13933444367982507207E6,
+        -2.53252307177582951285E6,
+        -2.01889141433532773231E6,
+        0.0,
+        0.0
+    };
+
+    sgngam = 1;
+
+    if( x < -34.0 ) {
+	    q = -x;
+	    w = lanczos_gamma(q); /* note this modifies sgngam! */
+	    p = floor(q);
+	    if( p == q ) {
+            lgsing:
+        		goto loverf;
+		}
+	    i = p;
+	    if( (i & 1) == 0 )
+		    sgngam = -1;
+	    else
+		    sgngam = 1;
+	    z = q - p;
+	    if( z > 0.5 ) {
+		    p += 1.0;
+		    z = p - q;
+		}
+	    z = q * sin( PI * z );
+	    if( z == 0.0 )
+		    goto lgsing;
+	    z = LOGPI - log( z ) - w;
+	    return( z );
+	}
+
+    if( x < 13.0 ) {
+	    z = 1.0;
+	    p = 0.0;
+	    u = x;
+	    while( u >= 3.0 ) {
+		    p -= 1.0;
+		    u = x + p;
+		    z *= u;
+		}
+	    while( u < 2.0 ) {
+		    if( u == 0.0 )
+			    goto lgsing;
+		    z /= u;
+		    p += 1.0;
+		    u = x + p;
+		}
+	    if( z < 0.0 ) {
+		    sgngam = -1;
+		    z = -z;
+		}
+	    else
+		    sgngam = 1;
+	    if( u == 2.0 )
+		    return( log(z) );
+	    p -= 2.0;
+	    x = x + p;
+	    p = x * polevl( x, B, 5 ) / p1evl( x, C, 6);
+	    return( log(z) + p );
+	}
+
+    if( x > MAXLGM ) {
+        loverf:
+	        return( sgngam * MAXNUM );
+	}
+
+    q = ( x - 0.5 ) * log(x) - x + LS2PI;
+    if( x > 1.0e8 )
+	    return( q );
+
+    p = 1.0/(x*x);
+    if( x >= 1000.0 )
+	    q += ((   7.9365079365079365079365e-4 * p
+		    - 2.7777777777777777777778e-3) *p
+		    + 0.0833333333333333333333) / x;
+    else
+	    q += polevl( p, A, 4 ) / x;
+    return( q );
+#else
+
+    double p, q, w, z, xx;
+    double nx, tx;
+    int i, direction;
+    int sgngamf = 1;
+
+     double B[8] = {
+        6.055172732649237E-004,
+        -1.311620815545743E-003,
+        2.863437556468661E-003,
+        -7.366775108654962E-003,
+        2.058355474821512E-002,
+        -6.735323259371034E-002,
+        3.224669577325661E-001,
+        4.227843421859038E-001
+    };
+
+    /* log gamma(x+1), -.25 < x < .25 */
+    double C[8] = {
+        1.369488127325832E-001,
+        -1.590086327657347E-001,
+        1.692415923504637E-001,
+        -2.067882815621965E-001,
+        2.705806208275915E-001,
+        -4.006931650563372E-001,
+        8.224670749082976E-001,
+        -5.772156501719101E-001
+    };
+
+    /* log( sqrt( 2*pi ) ) */
+    const double LS2PI  =  0.91893853320467274178;
+    const double MAXLGM = 2.035093e36;
+    const double PIINV =  0.318309886183790671538;
+    const double MAXNUMF = 3.4028234663852885981170418348451692544e38;
+    const double PIF = 3.141592653589793238;
+    /* Logarithm of gamma function */
+
+    sgngamf = 1;
+
+    xx = x;
+    if( xx < 0.0 ) {
+	    q = -xx;
+	    w = lgamma(q); /* note this modifies sgngam! */
+	    p = floor(q);
+	    if( p == q )
+		    goto loverf;
+	    i = p;
+
+	    if( (i & 1) == 0 )
+		    sgngamf = -1;
+	    else
+		    sgngamf = 1;
+
+	    z = q - p;
+
+	    if( z > 0.5 ) {
+		    p += 1.0;
+		    z = p - q;
+	    }
+	    z = q * sin( PIF * z );
+	    if( z == 0.0 )
+		    goto loverf;
+	    z = -log( PIINV*z ) - w;
+	    return( z );
+	}
+
+    if( x < 6.5 ) {
+	    direction = 0;
+	    z = 1.0;
+	    tx = x;
+	    nx = 0.0;
+	    if( x >= 1.5 ) {
+		    while( tx > 2.5 ) {
+			    nx -= 1.0;
+			    tx = x + nx;
+			    z *=tx;
+			}
+		    x += nx - 2.0;
+            iv1r5:
+		    p = x * polevl( x, B, 7 );
+		    goto cont;
+	}
+	if( x >= 1.25 ) {
+		z *= x;
+		x -= 1.0; /* x + 1 - 2 */
+		direction = 1;
+		goto iv1r5;
+	}
+	if( x >= 0.75 ) {
+		x -= 1.0;
+		p = x * polevl( x, C, 7 );
+		q = 0.0;
+		goto contz;
+	}
+	while( tx < 1.5 ) {
+		if( tx == 0.0 )
+			goto loverf;
+		z *=tx;
+		nx += 1.0;
+		tx = x + nx;
+	}
+	direction = 1;
+	x += nx - 2.0;
+	p = x * polevl( x, B, 7 );
+
+    cont:
+	    if( z < 0.0 ) {
+		    sgngamf = -1;
+		    z = -z;
+        }
+	    else {
+		    sgngamf = 1;
+		}
+	    q = log(z);
+	    if( direction )
+		    q = -q;
+    contz:
+	    return( p + q );
+	}
+
+    if( x > MAXLGM ) {
+        loverf:
+	    return( sgngamf * MAXNUMF );
+	}
+
+    /* Note, though an asymptotic formula could be used for x >= 3,
+    * there is cancellation error in the following if x < 6.5.  */
+    q = LS2PI - x;
+    q += ( x - 0.5 ) * log(x);
+
+    if( x <= 1.0e4 ) {
+	    z = 1.0/x;
+	    p = z * z;
+	    q += ((    6.789774945028216E-004 * p
+		 - 2.769887652139868E-003 ) * p
+		+  8.333316229807355E-002 ) * z;
+	}
+    return( q );
+#endif
+}
