Index: sasmodels/models/lib/lanczos_gamma.c
===================================================================
--- sasmodels/models/lib/lanczos_gamma.c	(revision 5c1d3414d3ff3bc81b5e9cfefb15dded91ee5d3e)
+++ sasmodels/models/lib/lanczos_gamma.c	(revision 5c1d3414d3ff3bc81b5e9cfefb15dded91ee5d3e)
@@ -0,0 +1,21 @@
+double lanczos_gamma(double q);
+double lanczos_gamma(double q)
+{
+    // Lanczos approximation to the Gamma function.
+
+    double x,y,tmp,ser;
+    double coeff[6]=
+        {76.18009172947146,     -86.50532032941677,
+         24.01409824083091,     -1.231739572450155,
+          0.1208650973866179e-2,-0.5395239384953e-5};
+
+    y=x=q;
+    tmp  = x+5.5;
+    tmp -= (x+0.5)*log(tmp);
+    ser  = 1.000000000190015;
+    for (int j=0; j<=5; j++) {
+        y+=1.0;
+        ser += coeff[j]/y;
+    }
+    return -tmp+log(2.5066282746310005*ser/x);
+}
Index: sasmodels/models/mass_fractal.c
===================================================================
--- sasmodels/models/mass_fractal.c	(revision 5c1d3414d3ff3bc81b5e9cfefb15dded91ee5d3e)
+++ sasmodels/models/mass_fractal.c	(revision 5c1d3414d3ff3bc81b5e9cfefb15dded91ee5d3e)
@@ -0,0 +1,68 @@
+double form_volume(double radius);
+
+double Iq(double q,
+          double radius,
+          double mass_dim,
+          double cutoff_length);
+
+double Iqxy(double qx, double qy,
+          double radius,
+          double mass_dim,
+          double cutoff_length);
+
+
+static double _mass_fractal_kernel(double q,
+          double radius,
+          double mass_dim,
+          double cutoff_length)
+{
+
+    //calculate P(q) for the spherical subunits; not normalized
+    //double pq = pow((3.0*(sin(q*radius) - q*radius*cos(q*radius))/pow((q*radius),3)),2);
+    double pq = J1c(q*radius);
+
+    //pq = 1.0;
+
+    //calculate S(q)
+    double mmo = mass_dim-1.0;
+    double sq = exp(lanczos_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;
+
+    //sq = 1.0;
+
+    //combine and return
+    double result = pq * sq;
+
+    return result;
+}
+double form_volume(double radius){
+
+    return 1.333333333333333*M_PI*radius*radius*radius;
+}
+
+double Iq(double q,
+          double radius,
+          double mass_dim,
+          double cutoff_length)
+{
+    return _mass_fractal_kernel(q,
+           radius,
+           mass_dim,
+           cutoff_length);
+}
+
+// Iqxy is never called since no orientation or magnetic parameters.
+double Iqxy(double qx, double qy,
+          double radius,
+          double mass_dim,
+          double cutoff_length)
+{
+    double q = sqrt(qx*qx + qy*qy);
+    return _mass_fractal_kernel(q,
+           radius,
+           mass_dim,
+           cutoff_length);
+}
+
Index: sasmodels/models/mass_fractal.py
===================================================================
--- sasmodels/models/mass_fractal.py	(revision 5c1d3414d3ff3bc81b5e9cfefb15dded91ee5d3e)
+++ sasmodels/models/mass_fractal.py	(revision 5c1d3414d3ff3bc81b5e9cfefb15dded91ee5d3e)
@@ -0,0 +1,119 @@
+r"""
+Calculates the scattering from fractal-like aggregates based on
+the Mildner reference.
+
+Definition
+----------
+
+The scattering intensity $I(q)$ is calculated as
+
+.. 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(D_m-1)\zeta^{D_m-1}}{\left[1+(q\zeta)^2
+    \right]^{(D_m-1)/2}}
+    \frac{sin\left[(D_m - 1) 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
+
+where $R$ is the radius of the building block, $D_m$ is the **mass** 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.
+
+.. note::
+
+    The mass fractal dimension ( $D_m$ ) is only
+    valid if $0 < mass_dim < 6$. It is also only valid over a limited
+    $q$ range (see the reference for details).
+
+.. figure:: img/mass_fractal_1d.jpg
+
+    1D plot using the default values.
+
+Reference
+---------
+
+D Mildner and P Hall, *J. Phys. D: Appl. Phys.*,  19 (1986) 1535-1545 Equation(9)
+
+
+"""
+
+from numpy import inf
+
+name = "mass_fractal"
+title = "Mass Fractal model"
+description = """
+        The scattering intensity  I(x) = scale*P(x)*S(x) + background, where
+        scale = scale_factor  * V * delta^(2)
+        p(x)=  F(x*radius)^(2)
+        F(x) = 3*[sin(x)-x cos(x)]/x**3
+        S(x) = [(gamma(Dm-1)*colength^(Dm-1)*[1+(x^2*colength^2)]^((1-Dm)/2)
+        * sin[(Dm-1)*arctan(x*colength)])/x]
+        where delta = sldParticle -sldSolv.
+        radius       =  Particle radius
+        mass_dim  =  Mass fractal dimension
+        cutoff_length  =  Cut-off length
+        background   =  background
+        Ref.:Mildner, Hall,J Phys D Appl Phys(1986), 9, 1535-1545
+        Note I: This model is valid for 1<mass_dim<6.
+        Note II: This model is not in absolute scale.
+        """
+category = "shape-independent"
+
+#             ["name", "units", default, [lower, upper], "type","description"],
+parameters = [["radius",        "Ang",  10.0, [0.0, inf], "",
+               "Particle radius"],
+              ["mass_dim",      "",      1.9, [1.0, 6.0],   "",
+               "Mass fractal dimension"],
+              ["cutoff_length", "Ang", 100.0, [0.0, inf], "",
+               "Cut-off length"],
+              ]
+
+
+source = ["lib/J1c.c", "lib/lanczos_gamma.c", "mass_fractal.c"]
+
+demo = dict(scale=1, background=0,
+            radius=10.0,
+            mass_dim=1.9,
+            cutoff_length=2.3)
+
+oldname = 'MassFractalModel'
+oldpars = dict(radius='radius',
+               mass_dim='mass_dim',
+               cutoff_length='co_length')
+
+tests = [[{'radius':        2.0,
+           'mass_dim':      3.3,
+           'cutoff_length': 1.0,
+           }, 0.001, 0.18462699016],
+
+         [{'radius':        1.0,
+           'mass_dim':      1.3,
+           'cutoff_length': 1.0,
+           'background':    0.8,
+           }, 0.001, 1.16539753641],
+
+         [{'radius':        1.0,
+           'mass_dim':      2.3,
+           'cutoff_length': 1.0,
+           'scale':        10.0,
+           }, 0.051, 0.000169548800377],
+         ]
Index: sasmodels/models/surface_fractal.c
===================================================================
--- sasmodels/models/surface_fractal.c	(revision 513efc509d218370e6a8ca9cabf24f087a2de902)
+++ sasmodels/models/surface_fractal.c	(revision 5c1d3414d3ff3bc81b5e9cfefb15dded91ee5d3e)
@@ -10,26 +10,4 @@
             double surface_dim,
             double cutoff_length);
-
-static double _gamln(double q)
-{
-    // Lanczos approximation to the Gamma function.
-    // Should be refactored out to lib/, if used elsewhere.
-    double x,y,tmp,ser;
-    double coeff[6]=
-        {76.18009172947146,     -86.50532032941677,
-         24.01409824083091,     -1.231739572450155,
-          0.1208650973866179e-2,-0.5395239384953e-5};
-    int j;
-
-    y=x=q;
-    tmp  = x+5.5;
-    tmp -= (x+0.5)*log(tmp);
-    ser  = 1.000000000190015;
-    for (j=0; j<=5; j++) {
-        y+=1.0;
-        ser += coeff[j]/y;
-    }
-    return -tmp+log(2.5066282746310005*ser/x);
-}
 
 static double _surface_fractal_kernel(double q,
@@ -48,5 +26,5 @@
     //calculate S(q)
     mmo = 5.0 - surface_dim;
-    sq  = exp(_gamln(mmo))*sin(-(mmo)*atan(q*cutoff_length));
+    sq  = exp(lanczos_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));
Index: sasmodels/models/surface_fractal.py
===================================================================
--- sasmodels/models/surface_fractal.py	(revision 513efc509d218370e6a8ca9cabf24f087a2de902)
+++ sasmodels/models/surface_fractal.py	(revision 5c1d3414d3ff3bc81b5e9cfefb15dded91ee5d3e)
@@ -90,5 +90,5 @@
 
 
-source = ["lib/J1c.c", "surface_fractal.c"]
+source = ["lib/J1c.c", "lib/lanczos_gamma.c", "surface_fractal.c"]
 
 demo = dict(scale=1, background=0,
