Index: sasmodels/models/fractal.c
===================================================================
--- sasmodels/models/fractal.c	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
+++ sasmodels/models/fractal.c	(revision 6d96b66e1e58166fbc135a8df1ddf8d37f64b54a)
@@ -1,54 +1,45 @@
-double Iq(double q,
-          double volfraction,
-          double radius,
-          double fractal_dim,
-          double cor_length,
-          double sld_block,
-          double sld_solvent);
+#define INVALID(p) (p.fractal_dim <= 0.0)
 
-double Iq(double q,
-          double volfraction,
-          double radius,
-          double fractal_dim,
-          double cor_length,
-          double sld_block,
-          double sld_solvent)
+static double
+Iq(double q,
+   double volfraction,
+   double radius,
+   double fractal_dim,
+   double cor_length,
+   double sld_block,
+   double sld_solvent)
 {
-    double qr,r0,Df,corr,phi,sldp,sldm;
-    double pq,sq,inten;
+    //calculate P(q) for the spherical subunits
+    const double pq = M_4PI_3*cube(radius) * square((sld_block-sld_solvent)*sph_j1c(q*radius));
     
-     // Actively check the argument - needed for mass fractal - is it needie
-     //here?
-    if (fractal_dim <= 0.0){
-       return 0.0;
+    //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;
     }
-   
-    phi = volfraction;        // volume fraction of building block spheres...
-    r0 = radius;     //  radius of building block
-    Df = fractal_dim;     //  fractal dimension
-    corr = cor_length;       //  correlation length of fractal-like aggregates
-    sldp = sld_block;       // SLD of building block
-    sldm = sld_solvent;       // SLD of matrix or solution
- 
-     qr=q*r0;
-    
-    //calculate P(q) for the spherical subunits
-    pq = phi*M_4PI_3*r0*r0*r0*(sldp-sldm)*(sldp-sldm)*sph_j1c(qr)*sph_j1c(qr);
-    
-    //calculate S(q)
-    sq = Df*sas_gamma(Df-1.0)*sin((Df-1.0)*atan(q*corr));
-    sq /= pow(qr,Df) * pow((1.0 + 1.0/(q*corr)/(q*corr)),((Df-1.0)/2.0));
-    sq += 1.0;
-    
-    //combine, scale to units cm-1 sr-1 (assuming data on absolute scale)
-    //and return
-    inten = pq*sq;
-    // convert I(1/A) to (1/cm)
-    inten *= 1.0e8;
-    //convert rho^2 in 10^-6 1/A to 1/A
-    inten *= 1.0e-12;    
-    
-    
-    return(inten);
+
+    // scale to units cm-1 sr-1 (assuming data on absolute scale)
+    //    convert I(1/A) to (1/cm)  => 1e8 * I(q)
+    //    convert rho^2 in 10^-6 1/A to 1/A  => 1e-12 * I(q)
+    //    combined: 1e-4 * I(q)
+
+    return 1.e-4 * volfraction * pq * sq;
 }
 
Index: sasmodels/models/fractal.py
===================================================================
--- sasmodels/models/fractal.py	(revision 785cbec0749457756eec925df00538ade654c842)
+++ sasmodels/models/fractal.py	(revision 6d96b66e1e58166fbc135a8df1ddf8d37f64b54a)
@@ -84,5 +84,5 @@
               ["radius",    "Ang",  5.0, [0.0, inf], "",
                "radius of particles"],
-              ["fractal_dim",      "",  2.0, [0.0, 6.0], "",
+              ["fractal_dim",      "",  2.0, [1.0, 6.0], "",
                "fractal dimension"],
               ["cor_length", "Ang", 100.0, [0.0, inf], "",
Index: sasmodels/models/fractal_core_shell.c
===================================================================
--- sasmodels/models/fractal_core_shell.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
+++ sasmodels/models/fractal_core_shell.c	(revision 6d96b66e1e58166fbc135a8df1ddf8d37f64b54a)
@@ -27,21 +27,32 @@
 
 
-   double intensity = core_shell_kernel(q,
-                              radius,
-                              thickness,
-                              core_sld,
-                              shell_sld,
-                              solvent_sld);
+    const double pq = core_shell_kernel(q, radius, thickness,
+                                        core_sld, shell_sld, solvent_sld);
+
+
     //calculate S(q)
-    double frac_1 = fractal_dim-1.0;
-    double qr = q*radius;
+    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;
+    }
 
-    double t1 = fractal_dim*sas_gamma(frac_1)*sin(frac_1*atan(q*cor_length));
-    double t2 = (1.0 + 1.0/(q*cor_length)/(q*cor_length));
-    double t3 = pow(qr, fractal_dim) * pow(t2, (frac_1/2.0));
-    double sq = t1/t3;
-    sq += 1.0;
-
-    return sq*intensity*volfraction;
+    // Note: core_shell_kernel already performs the 1e-4 unit conversion
+    return volfraction * sq * pq;
 }
 
Index: sasmodels/models/fractal_core_shell.py
===================================================================
--- sasmodels/models/fractal_core_shell.py	(revision 4962519cc1910422a467e4f643dba223623e0e09)
+++ sasmodels/models/fractal_core_shell.py	(revision 6d96b66e1e58166fbc135a8df1ddf8d37f64b54a)
@@ -63,12 +63,13 @@
 #   ["name", "units", default, [lower, upper], "type","description"],
 parameters = [
-    ["radius",      "Ang",        60.0, [0, inf],    "volume", "Sphere core radius"],
-    ["thickness",   "Ang",        10.0, [0, inf],    "volume", "Sphere shell thickness"],
+    ["radius",      "Ang",        60.0, [0.0, inf],  "volume", "Sphere core radius"],
+    ["thickness",   "Ang",        10.0, [0.0, inf],  "volume", "Sphere shell thickness"],
     ["sld_core",    "1e-6/Ang^2", 1.0,  [-inf, inf], "sld",    "Sphere core scattering length density"],
     ["sld_shell",   "1e-6/Ang^2", 2.0,  [-inf, inf], "sld",    "Sphere shell scattering length density"],
     ["sld_solvent", "1e-6/Ang^2", 3.0,  [-inf, inf], "sld",    "Solvent scattering length density"],
-    ["volfraction", "",           1.0,  [0, inf],    "",       "Volume fraction of building block spheres"],
-    ["fractal_dim",    "",           2.0,  [-inf, inf], "",       "Fractal dimension"],
-    ["cor_length",  "Ang",      100.0,  [0, inf],    "",       "Correlation length of fractal-like aggregates"]]
+    ["volfraction", "",           1.0,  [0.0, inf],  "",       "Volume fraction of building block spheres"],
+    ["fractal_dim",    "",        2.0,  [1.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
 
Index: sasmodels/models/mass_fractal.c
===================================================================
--- sasmodels/models/mass_fractal.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
+++ sasmodels/models/mass_fractal.c	(revision 6d96b66e1e58166fbc135a8df1ddf8d37f64b54a)
@@ -1,48 +1,31 @@
-double form_volume(double radius);
+#define INVALID(p) (p.fractal_dim_mass < 1.0)
 
-double Iq(double q,
-          double radius,
-          double fractal_dim_mass,
-          double cutoff_length);
+static double
+Iq(double q, double radius, double fractal_dim_mass, double cutoff_length)
+{
+    //calculate P(q)
+    const double pq = square(sph_j1c(q*radius));
 
-static double _mass_fractal_kernel(double q,
-          double radius,
-          double fractal_dim_mass,
-          double cutoff_length)
-{
-    // Actively check the argument.
-    if (fractal_dim_mass <= 1.0){
-       return 0.0;
+    //calculate S(q)
+    // S(q) = gamma(D-1) sin((D-1)atan(q c))/q c^(D-1) (1+(q c)^2)^(-(D-1)/2)
+    // lim D->1 [gamma(D-1) sin((D-1) a)] = a
+    // lim q->0 [sin((D-1) atan(q c))/q] = (D-1) c
+    // lim q,D->0,1 [gamma(D-1) sin((D-1)atan(q c))/q] = c
+    double sq;
+    if (q > 0. && fractal_dim_mass > 1.) {
+        const double Dm1 = fractal_dim_mass - 1.0;
+        const double t1 = sas_gamma(Dm1)*sin(Dm1*atan(q*cutoff_length));
+        const double t2 = pow(cutoff_length, Dm1);
+        const double t3 = pow(1.0 + square(q*cutoff_length), -0.5*Dm1);
+        sq = t1 * t2 * t3 / q;
+    } else if (q > 0.) {
+        sq = atan(q*cutoff_length)/q;
+    } else if (fractal_dim_mass > 1.) {
+        const double D = fractal_dim_mass;
+        sq = pow(cutoff_length, D) * sas_gamma(D);
+    } else {
+        sq = cutoff_length;
     }
 
-    //calculate P(q)
-    double pq = sph_j1c(q*radius);
-    pq = pq*pq;
-
-    //calculate S(q)
-    double mmo = fractal_dim_mass-1.0;
-    double 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;
-
-    //combine and return
-    double result = pq * sq;
-
-    return result;
+    return pq * sq;
 }
-double form_volume(double radius){
-
-    return M_4PI_3*radius*radius*radius;
-}
-
-double Iq(double q,
-          double radius,
-          double fractal_dim_mass,
-          double cutoff_length)
-{
-    return _mass_fractal_kernel(q,
-           radius,
-           fractal_dim_mass,
-           cutoff_length);
-}
Index: sasmodels/models/mass_fractal.py
===================================================================
--- sasmodels/models/mass_fractal.py	(revision a8072067ae8a224e45d23ba71f663aadbe1ba4a9)
+++ sasmodels/models/mass_fractal.py	(revision 6d96b66e1e58166fbc135a8df1ddf8d37f64b54a)
@@ -78,9 +78,10 @@
 
 # pylint: disable=bad-whitespace, line-too-long
-#             ["name", "units", default, [lower, upper], "type","description"],
-parameters = [["radius",        "Ang",  10.0, [0.0, inf], "", "Particle radius"],
-              ["fractal_dim_mass",      "",      1.9, [1.0, 6.0], "", "Mass fractal dimension"],
-              ["cutoff_length", "Ang", 100.0, [0.0, inf], "", "Cut-off length"],
-             ]
+#   ["name", "units", default, [lower, upper], "type","description"],
+parameters = [
+    ["radius",           "Ang",  10.0, [0.0, inf], "", "Particle radius"],
+    ["fractal_dim_mass", "",      1.9, [1.0, 6.0], "", "Mass fractal dimension"],
+    ["cutoff_length",    "Ang", 100.0, [0.0, inf], "", "Cut-off length"],
+]
 # pylint: enable=bad-whitespace, line-too-long
 
Index: sasmodels/models/mass_surface_fractal.c
===================================================================
--- sasmodels/models/mass_surface_fractal.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
+++ sasmodels/models/mass_surface_fractal.c	(revision 6d96b66e1e58166fbc135a8df1ddf8d37f64b54a)
@@ -1,11 +1,5 @@
-double form_volume(double radius);
-
-double Iq(double q,
-          double fractal_dim_mass,
-          double fractal_dim_surf,
-          double rg_cluster,
-          double rg_primary);
-
-static double _mass_surface_fractal_kernel(double q,
+#define INVALID(p) (p.fractal_dim_mass + p.fractal_dim_surf > 6.)
+static double
+Iq(double q,
           double fractal_dim_mass,
           double fractal_dim_surf,
@@ -14,40 +8,11 @@
 {
      //computation
-    double tot_dim = 6.0 - fractal_dim_surf - fractal_dim_mass;
-    fractal_dim_mass /= 2.0;
-    tot_dim /= 2.0;
+    const double Dm = 0.5*fractal_dim_mass;
+    const double Dt = 0.5*(6.0 - (fractal_dim_mass + fractal_dim_surf));
 
-    double rc_norm = rg_cluster * rg_cluster / (3.0 * fractal_dim_mass);
-    double rp_norm = rg_primary * rg_primary / (3.0 * tot_dim);
+    const double t1 = Dm==0. ? 1.0 : pow(1.0 + square(q*rg_cluster)/(3.0*Dm), -Dm);
+    const double t2 = Dt==0. ? 1.0 : pow(1.0 + square(q*rg_primary)/(3.0*Dt), -Dt);
+    const double form = t1*t2;
 
-    //x for P
-    double x_val1 = 1.0 +  q * q * rc_norm;
-    double x_val2 = 1.0 +  q * q * rp_norm;
-
-    double inv_form = pow(x_val1, fractal_dim_mass) * pow(x_val2, tot_dim);
-
-    //another singular
-    if (inv_form == 0.0) return 0.0;
-
-    double form_factor = 1.0;
-    form_factor /= inv_form;
-
-    return (form_factor);
+    return form;
 }
-double form_volume(double radius){
-
-    return M_4PI_3*radius*radius*radius;
-}
-
-double Iq(double q,
-          double fractal_dim_mass,
-          double fractal_dim_surf,
-          double rg_cluster,
-          double rg_primary)
-{
-    return _mass_surface_fractal_kernel(q,
-            fractal_dim_mass,
-            fractal_dim_surf,
-            rg_cluster,
-            rg_primary);
-}
Index: sasmodels/models/mass_surface_fractal.py
===================================================================
--- sasmodels/models/mass_surface_fractal.py	(revision 30fbe2eda4d6f97acc840bd5e605de3829caf3c7)
+++ sasmodels/models/mass_surface_fractal.py	(revision 6d96b66e1e58166fbc135a8df1ddf8d37f64b54a)
@@ -78,14 +78,11 @@
 
 # pylint: disable=bad-whitespace, line-too-long
-#             ["name", "units", default, [lower, upper], "type","description"],
-parameters = [["fractal_dim_mass",      "",    1.8, [1e-16, 6.0], "",
-               "Mass fractal dimension"],
-              ["fractal_dim_surf",   "",    2.3, [1e-16, 6.0], "",
-               "Surface fractal dimension"],
-              ["rg_cluster", "Ang",   86.7, [0.0, inf], "",
-               "Cluster radius of gyration"],
-              ["rg_primary", "Ang", 4000.,  [0.0, inf], "",
-               "Primary particle radius of gyration"],
-             ]
+#   ["name", "units", default, [lower, upper], "type","description"],
+parameters = [
+    ["fractal_dim_mass", "",      1.8, [0.0, 6.0], "", "Mass fractal dimension"],
+    ["fractal_dim_surf", "",      2.3, [0.0, 6.0], "", "Surface fractal dimension"],
+    ["rg_cluster",       "Ang",  86.7, [0.0, inf], "", "Cluster radius of gyration"],
+    ["rg_primary",       "Ang", 4000., [0.0, inf], "", "Primary particle radius of gyration"],
+]
 # pylint: enable=bad-whitespace, line-too-long
 
