Index: sasmodels/models/surface_fractal.c
===================================================================
--- sasmodels/models/surface_fractal.c	(revision a8072067ae8a224e45d23ba71f663aadbe1ba4a9)
+++ sasmodels/models/surface_fractal.c	(revision b716cc63882c24230b143debfbf45e545ba6bac2)
@@ -1,2 +1,5 @@
+// Don't need invalid test since fractal_dim_surf is not polydisperse
+// #define INVALID(v) (v.fractal_dim_surf <= 1.0 || v.fractal_dim_surf >= 3.0)
+
 double form_volume(double radius);
 
@@ -11,27 +14,28 @@
     double cutoff_length)
 {
-    double pq, sq, mmo, result;
+    // calculate P(q)
+    const double pq = square(sph_j1c(q*radius));
 
-    //Replaced the original formula with Taylor expansion near zero.
-    //pq = pow((3.0*(sin(q*radius) - q*radius*cos(q*radius))/pow((q*radius),3)),2);
+    // calculate S(q)
+    // Note: lim q->0 S(q) = -gamma(mmo) cutoff_length^mmo (mmo cutoff_length)
+    // however, the surface fractal formula is invalid outside the range
+    const double mmo = 5.0 - fractal_dim_surf;
+    const double sq = sas_gamma(mmo) * pow(cutoff_length, mmo)
+           * pow(1.0 + square(q*cutoff_length), -0.5*mmo)
+           * sin(-mmo * atan(q*cutoff_length)) / q;
 
-    pq = sph_j1c(q*radius);
-    pq = pq*pq;
+    // Empirically determined that the results are valid within this range.
+    // Above 1/r, the form starts to oscillate;  below
+    //const double result = (q > 5./(3-fractal_dim_surf)/cutoff_length) && q < 1./radius
+    //                      ? pq * sq : 0.);
 
-    //calculate S(q)
-    mmo = 5.0 - fractal_dim_surf;
-    sq  = sas_gamma(mmo)*sin(-(mmo)*atan(q*cutoff_length));
-    sq *= pow(cutoff_length, mmo);
-    sq /= pow((1.0 + (q*cutoff_length)*(q*cutoff_length)),(mmo/2.0));
-    sq /= q;
+    double result = pq * sq;
 
-    //combine and return
-    result = pq * sq;
-
-    return result;
+    // exclude negative results
+    return result > 0. ? result : 0.;
 }
-double form_volume(double radius){
-
-    return 1.333333333333333*M_PI*radius*radius*radius;
+double form_volume(double radius)
+{
+    return M_4PI_3*cube(radius);
 }
 
Index: sasmodels/models/surface_fractal.py
===================================================================
--- sasmodels/models/surface_fractal.py	(revision a8072067ae8a224e45d23ba71f663aadbe1ba4a9)
+++ sasmodels/models/surface_fractal.py	(revision b716cc63882c24230b143debfbf45e545ba6bac2)
@@ -10,37 +10,22 @@
 .. math::
 
-    I(q) = scale \times P(q)S(q) + background
-
-.. math::
-
-    P(q) = F(qR)^2
-
-.. math::
-
-    F(x) = \frac{3\left[sin(x)-xcos(x)\right]}{x^3}
-
-.. math::
-
-    S(q) = \frac{\Gamma(5-D_S)\zeta^{5-D_S}}{\left[1+(q\zeta)^2
-    \right]^{(5-D_S)/2}}
-    \frac{sin\left[(D_S - 5) tan^{-1}(q\zeta) \right]}{q}
-
-.. math::
-
-    scale = scale\_factor \times NV^2(\rho_{particle} - \rho_{solvent})^2
-
-.. math::
-
-    V = \frac{4}{3}\pi R^3
+    I(q) &= \text{scale} \times P(q)S(q) + \text{background} \\
+    P(q) &= F(qR)^2 \\
+    F(x) &= \frac{3\left[\sin(x)-x\cos(x)\right]}{x^3} \\
+    S(q) &= \Gamma(5-D_S)\xi^{\,5-D_S}\left[1+(q\xi)^2 \right]^{-(5-D_S)/2}
+            \sin\left[-(5-D_S) \tan^{-1}(q\xi) \right] q^{-1} \\
+    \text{scale} &= \phi N V^2(\rho_\text{particle} - \rho_\text{solvent})^2 \\
+    V &= \frac{4}{3}\pi R^3
 
 where $R$ is the radius of the building block, $D_S$ is the **surface** fractal
-dimension,| \zeta\|  is the cut-off length, $\rho_{solvent}$ is the scattering
-length density of the solvent,
-and $\rho_{particle}$ is the scattering length density of particles.
+dimension, $\xi$ is the cut-off length, $\rho_\text{solvent}$ is the scattering
+length density of the solvent, $\rho_\text{particle}$ is the scattering
+length density of particles and $\phi$ is the volume fraction of the particles.
 
 .. note::
-    The surface fractal dimension $D_s$ is only valid if $1<surface\_dim<3$.
-    It is also only valid over a limited $q$ range (see the reference for
-    details)
+
+    The surface fractal dimension is only valid if $1<D_S<3$. The result is
+    only valid over a limited $q$ range, $\tfrac{5}{3-D_S}\xi^{\,-1} < q < R^{-1}$.
+    See the reference for details.
 
 
@@ -89,5 +74,5 @@
 source = ["lib/sph_j1c.c", "lib/sas_gamma.c", "surface_fractal.c"]
 
-demo = dict(scale=1, background=0,
+demo = dict(scale=1, background=1e-5,
             radius=10, fractal_dim_surf=2.0, cutoff_length=500)
 
