Index: sasmodels/models/fractal.c
===================================================================
--- sasmodels/models/fractal.c	(revision 6d96b66e1e58166fbc135a8df1ddf8d37f64b54a)
+++ sasmodels/models/fractal.c	(revision 217590bce5ffcf2c4de704a4770a9cab003dd287)
@@ -1,3 +1,3 @@
-#define INVALID(p) (p.fractal_dim <= 0.0)
+#define INVALID(p) (p.fractal_dim < 0.0)
 
 static double
@@ -10,29 +10,9 @@
    double sld_solvent)
 {
+    const double sq = fractal_sq(q, radius, fractal_dim, cor_length);
+
     //calculate P(q) for the spherical subunits
-    const double pq = M_4PI_3*cube(radius) * square((sld_block-sld_solvent)*sph_j1c(q*radius));
-    
-    //calculate S(q),  using Teixeira, Eq(15)
-    double sq;
-    if (q > 0. && fractal_dim > 1.) {
-        // q>0, D>0
-        const double D = fractal_dim;
-        const double Dm1 = fractal_dim - 1.0;
-        // Note: for large Dm1, sin(Dm1*atan(q*cor_length) can go negative
-        const double t1 = D*sas_gamma(Dm1)*sin(Dm1*atan(q*cor_length));
-        const double t2 = pow(q*radius, -D);
-        const double t3 = pow(1.0 + 1.0/square(q*cor_length), -0.5*Dm1);
-        sq = 1.0 + t1 * t2 * t3;
-    } else if (q > 0.) {
-        // q>0, D=1
-        sq = 1.0 + atan(q*cor_length) / (q*radius);
-    } else if (fractal_dim > 1.) {
-        // q=0, D>1
-        const double D = fractal_dim;
-        sq = 1.0 + pow(cor_length/radius, D)*sas_gamma(D+1.0);
-    } else {
-        // q=0, D=1
-        sq = 1.0 + cor_length/radius;
-    }
+    const double V = M_4PI_3*cube(radius);
+    const double pq = V * square((sld_block-sld_solvent)*sph_j1c(q*radius));
 
     // scale to units cm-1 sr-1 (assuming data on absolute scale)
@@ -41,5 +21,5 @@
     //    combined: 1e-4 * I(q)
 
-    return 1.e-4 * volfraction * pq * sq;
+    return 1.e-4 * volfraction * sq * pq;
 }
 
Index: sasmodels/models/fractal.py
===================================================================
--- sasmodels/models/fractal.py	(revision 6d96b66e1e58166fbc135a8df1ddf8d37f64b54a)
+++ sasmodels/models/fractal.py	(revision 217590bce5ffcf2c4de704a4770a9cab003dd287)
@@ -84,5 +84,5 @@
               ["radius",    "Ang",  5.0, [0.0, inf], "",
                "radius of particles"],
-              ["fractal_dim",      "",  2.0, [1.0, 6.0], "",
+              ["fractal_dim",      "",  2.0, [0.0, 6.0], "",
                "fractal dimension"],
               ["cor_length", "Ang", 100.0, [0.0, inf], "",
@@ -95,5 +95,5 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/sph_j1c.c", "lib/sas_gamma.c", "fractal.c"]
+source = ["lib/sph_j1c.c", "lib/sas_gamma.c", "lib/fractal_sq.c", "fractal.c"]
 
 demo = dict(volfraction=0.05,
Index: sasmodels/models/fractal_core_shell.c
===================================================================
--- sasmodels/models/fractal_core_shell.c	(revision 6d96b66e1e58166fbc135a8df1ddf8d37f64b54a)
+++ sasmodels/models/fractal_core_shell.c	(revision 217590bce5ffcf2c4de704a4770a9cab003dd287)
@@ -1,55 +1,22 @@
-double form_volume(double radius, double thickness);
-
-double Iq(double q,
-          double radius,
-          double thickness,
-          double core_sld,
-          double shell_sld,
-          double solvent_sld,
-          double volfraction,
-          double fractal_dim,
-          double cor_length);
-
-double form_volume(double radius, double thickness)
+static double
+form_volume(double radius, double thickness)
 {
     return M_4PI_3 * cube(radius + thickness);
 }
 
-double Iq(double q,
-          double radius,
-          double thickness,
-          double core_sld,
-          double shell_sld,
-          double solvent_sld,
-          double volfraction,
-          double fractal_dim,
-          double cor_length) {
-
-
+static double
+Iq(double q,
+   double radius,
+   double thickness,
+   double core_sld,
+   double shell_sld,
+   double solvent_sld,
+   double volfraction,
+   double fractal_dim,
+   double cor_length)
+{
+    const double sq = fractal_sq(q, radius, fractal_dim, cor_length);
     const double pq = core_shell_kernel(q, radius, thickness,
                                         core_sld, shell_sld, solvent_sld);
-
-
-    //calculate S(q)
-    double sq;
-    if (q > 0. && fractal_dim > 1.) {
-        // q>0, D>0
-        const double D = fractal_dim;
-        const double Dm1 = fractal_dim - 1.0;
-        const double t1 = D*sas_gamma(Dm1)*sin((Dm1)*atan(q*cor_length));
-        const double t2 = pow(q*radius, -D);
-        const double t3 = pow(1.0 + 1.0/square(q*cor_length), -0.5*Dm1);
-        sq = 1.0 + t1 * t2 * t3;
-    } else if (q > 0.) {
-        // q>0, D=1
-        sq = 1.0 + atan(q*cor_length) / (q*radius);
-    } else if (fractal_dim > 1.) {
-        // q=0, D>1
-        const double D = fractal_dim;
-        sq = 1.0 + pow(cor_length/radius, D)*sas_gamma(D+1.0);
-    } else {
-        // q=0, D=1
-        sq = 1.0 + cor_length/radius;
-    }
 
     // Note: core_shell_kernel already performs the 1e-4 unit conversion
Index: sasmodels/models/fractal_core_shell.py
===================================================================
--- sasmodels/models/fractal_core_shell.py	(revision 6d96b66e1e58166fbc135a8df1ddf8d37f64b54a)
+++ sasmodels/models/fractal_core_shell.py	(revision 217590bce5ffcf2c4de704a4770a9cab003dd287)
@@ -69,10 +69,11 @@
     ["sld_solvent", "1e-6/Ang^2", 3.0,  [-inf, inf], "sld",    "Solvent scattering length density"],
     ["volfraction", "",           1.0,  [0.0, inf],  "",       "Volume fraction of building block spheres"],
-    ["fractal_dim",    "",        2.0,  [1.0, 6.0],  "",       "Fractal dimension"],
+    ["fractal_dim",    "",        2.0,  [0.0, 6.0],  "",       "Fractal dimension"],
     ["cor_length",  "Ang",      100.0,  [0.0, inf],  "",       "Correlation length of fractal-like aggregates"],
 ]
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/sph_j1c.c", "lib/sas_gamma.c", "lib/core_shell.c", "fractal_core_shell.c"]
+source = ["lib/sph_j1c.c", "lib/sas_gamma.c", "lib/core_shell.c",
+          "lib/fractal_sq.c", "fractal_core_shell.c"]
 
 demo = dict(scale=0.05,
Index: sasmodels/models/lib/fractal_sq.c
===================================================================
--- sasmodels/models/lib/fractal_sq.c	(revision 217590bce5ffcf2c4de704a4770a9cab003dd287)
+++ sasmodels/models/lib/fractal_sq.c	(revision 217590bce5ffcf2c4de704a4770a9cab003dd287)
@@ -0,0 +1,28 @@
+static double
+fractal_sq(double q, double radius, double fractal_dim, double cor_length)
+{
+    //calculate S(q),  using Teixeira, Eq(15)
+    // mathematica query to check limiting conditions:
+    //    lim x->0 of [ x gamma(x-1) sin(arctan(q c (x-1))) (q r)^(-x) (1 + 1/(q c)^2)^((1-x)/2) ]
+    // Note: gamma(x) may be unreliable for x<0, so the gamma(D-1) is risky.
+    // We instead transform D*gamma(D-1) into gamma(D+1)/(D-1).
+    double term;
+    if (q == 0.) {
+        const double D = fractal_dim;
+        term = pow(cor_length/radius, D)*sas_gamma(D+1.);
+    } else if (fractal_dim == 0.) {
+        term = 1.0;
+    } else if (fractal_dim == 1.) {
+        term = atan(q*cor_length)/(q*radius);
+    } else {
+        // q>0, D>0
+        const double D = fractal_dim;
+        const double Dm1 = fractal_dim - 1.0;
+        // Note: for large Dm1, sin(Dm1*atan(q*cor_length) can go negative
+        const double t1 = sas_gamma(D+1.)/Dm1*sin(Dm1*atan(q*cor_length));
+        const double t2 = pow(q*radius, -D);
+        const double t3 = pow(1.0 + 1.0/square(q*cor_length), -0.5*Dm1);
+        term = t1 * t2 * t3;
+    }
+    return 1.0 + term;
+}
Index: sasmodels/models/lib/sas_gamma.c
===================================================================
--- sasmodels/models/lib/sas_gamma.c	(revision 1596de320c8f962693aa46eb70d0b5bec0e6c601)
+++ sasmodels/models/lib/sas_gamma.c	(revision 217590bce5ffcf2c4de704a4770a9cab003dd287)
@@ -137,3 +137,16 @@
 
 
-inline double sas_gamma( double x) { return tgamma(x+1)/x; }
+inline double sas_gamma(double x) {
+#if 1
+    // Fast reliable gamma in [-n,171], where n is a small integer
+    double norm = 1.;
+    while (x < 1.) {
+        norm *= x;
+        x += 1.0;
+    }
+    return tgamma(x)/norm;
+#else
+    // Fast reliable gamma in [0,170]
+    return tgamma(x+1)/x;
+#endif
+}
