Index: sasmodels/convert.py
===================================================================
--- sasmodels/convert.py	(revision 256dfe1527a1dc32e1e2b9114144fb1c47f2b256)
+++ sasmodels/convert.py	(revision c5ac2b2355ae3705fd04e87490d99bfe4140d545)
@@ -5,4 +5,5 @@
 
 from os.path import join as joinpath, abspath, dirname
+import math
 import warnings
 import json
@@ -34,5 +35,4 @@
     'vesicle',
     'multilayer_vesicle',
-    'core_multi_shell',
 ]
 
@@ -287,3 +287,5 @@
         elif name == 'poly_gauss_coil':
             pars['i_zero'] = 1
-
+        elif name == 'core_multi_shell':
+            pars['n'] = max(math.ceil(pars['n']), 4)
+
Index: sasmodels/modelinfo.py
===================================================================
--- sasmodels/modelinfo.py	(revision 256dfe1527a1dc32e1e2b9114144fb1c47f2b256)
+++ sasmodels/modelinfo.py	(revision c5ac2b2355ae3705fd04e87490d99bfe4140d545)
@@ -526,5 +526,5 @@
         Return the list of parameters for the given data type.
 
-        Vector parameters are expanded as in place.  If multiple parameters
+        Vector parameters are expanded in place.  If multiple parameters
         share the same vector length, then the parameters will be interleaved
         in the result.  The control parameters come first.  For example,
Index: sasmodels/models/core_multi_shell.c
===================================================================
--- sasmodels/models/core_multi_shell.c	(revision abdd01cc0315960ae6e5951909d2866515e37bae)
+++ sasmodels/models/core_multi_shell.c	(revision c5ac2b2355ae3705fd04e87490d99bfe4140d545)
@@ -1,152 +1,2 @@
-FourShell(double dp[], double q)
-{
-    // variables are:
-    //[0] scale factor
-    //[1] radius of core [ï¿œ]
-    //[2] SLD of the core   [ï¿œ-2]
-    //[3] thickness of shell 1 [ï¿œ]
-    //[4] SLD of shell 1
-    //[5] thickness of shell 2 [ï¿œ]
-    //[6] SLD of shell 2
-    //[7] thickness of shell 3
-    //[8] SLD of shell 3
-    //[9] thickness of shell 3
-    //[10] SLD of shell 3
-    //[11] SLD of solvent
-    //[12] background   [cm-1]
-    
-    double x,pi;
-    double scale,rcore,thick1,rhocore,rhoshel1,rhosolv,bkg;     //my local names
-    double bes,f,vol,qr,contr,f2;
-    double rhoshel2,thick2,rhoshel3,thick3,rhoshel4,thick4;
-    
-    pi = 4.0*atan(1.0);
-    x=q;
-    
-    scale = dp[0];
-    rcore = dp[1];
-    rhocore = dp[2];
-    thick1 = dp[3];
-    rhoshel1 = dp[4];
-    thick2 = dp[5];
-    rhoshel2 = dp[6];   
-    thick3 = dp[7];
-    rhoshel3 = dp[8];
-    thick4 = dp[9];
-    rhoshel4 = dp[10];  
-    rhosolv = dp[11];
-    bkg = dp[12];
-    
-        // core first, then add in shells
-    qr=x*rcore;
-    contr = rhocore-rhoshel1;
-    if(qr == 0){
-        bes = 1.0;
-    }else{
-        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
-    }
-    vol = 4.0*pi/3.0*rcore*rcore*rcore;
-    f = vol*bes*contr;
-    //now the shell (1)
-    qr=x*(rcore+thick1);
-    contr = rhoshel1-rhoshel2;
-    if(qr == 0){
-        bes = 1.0;
-    }else{
-        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
-    }
-    vol = 4.0*pi/3.0*(rcore+thick1)*(rcore+thick1)*(rcore+thick1);
-    f += vol*bes*contr;
-    //now the shell (2)
-    qr=x*(rcore+thick1+thick2);
-    contr = rhoshel2-rhoshel3;
-    if(qr == 0){
-        bes = 1.0;
-    }else{
-        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
-    }
-    vol = 4.0*pi/3.0*(rcore+thick1+thick2)*(rcore+thick1+thick2)*(rcore+thick1+thick2);
-    f += vol*bes*contr;
-    //now the shell (3)
-    qr=x*(rcore+thick1+thick2+thick3);
-    contr = rhoshel3-rhoshel4;
-    if(qr == 0){
-        bes = 1.0;
-    }else{
-        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
-    }
-    vol = 4.0*pi/3.0*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3);
-    f += vol*bes*contr;
-    //now the shell (4)
-    qr=x*(rcore+thick1+thick2+thick3+thick4);
-    contr = rhoshel4-rhosolv;
-    if(qr == 0){
-        bes = 1.0;
-    }else{
-        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
-    }
-    vol = 4.0*pi/3.0*(rcore+thick1+thick2+thick3+thick4)*(rcore+thick1+thick2+thick3+thick4)*(rcore+thick1+thick2+thick3+thick4);
-    f += vol*bes*contr;
-    
-        
-    // normalize to particle volume and rescale from [ï¿œ-1] to [cm-1]
-    f2 = f*f/vol*1.0e8;
-    
-    //scale if desired
-    f2 *= scale;
-    // then add in the background
-    f2 += bkg;
-    
-    return(f2);
-}
-
-// above is cut and past with no changes from libigor four shell
-static double
-f_exp(double q, double r, double sld_in, double sld_out,
-    double thickness, double A)
-{
-  const double vol = 4.0/3.0 * M_PI * r * r * r;
-  const double qr = q * r;
-  const double alpha = A * r/thickness;
-  const double bes = sph_j1c(qr);
-  const double B = (sld_out - sld_in)/expm1(A);
-  const double C = sld_in - B;
-  double fun;
-  if (qr == 0.0) {
-    fun = 1.0;
-  } else if (fabs(A) > 0.0) {
-    const double qrsq = qr*qr;
-    const double alphasq = alpha*alpha;
-    const double sumsq = alphasq + qrsq;
-    double sinqr, cosqr;
-    SINCOS(qr, sinqr, cosqr);
-    fun = -3.0(
-            ((alphasq - qrsq)*sinqr/qr - 2.0*alpha*cosqr) / sumsq
-                - (alpha*sinqr/qr - cosqr)
-        ) / sumsq;
-  } else {
-    fun = bes;
-  }
-  return vol * (B*fun + C*bes);
-}
-
-static double
-f_linear(double q, double r, double sld, double slope)
-{
-  const double vol = 4.0/3.0 * M_PI * r * r * r;
-  const double qr = q * r;
-  const double bes = sph_j1c(qr);
-  double fun = 0.0;
-  if (qr > 0.0) {
-    const double qrsq = qr*qr;
-    double sinqr, cosqr;
-    SINCOS(qr, sinqr, cosqr);
-    // Jae-He's code seems to simplify to this
-    //     fun = 3.0 * slope * r * (2.0*qr*sinqr - (qrsq-2.0)*cosqr)/(qrsq*qrsq);
-    // Rederiving the math, we get the following instead:
-    fun = 3.0 * slope * r * (2.0*cosqr + qr*sinqr)/(qrsq*qrsq);
-  }
-  return vol * (sld*bes + fun);
-}
 
 static double
@@ -154,47 +4,38 @@
 {
   const double bes = sph_j1c(q * r);
-  const double vol = 4.0/3.0 * M_PI * r * r * r;
+  const double vol = M_4PI_3 * cube(r);
   return sld * vol * bes;
 }
 
-static double
+double
+form_volume(double core_radius, double n, double thickness[]);
+double
 form_volume(double core_radius, double n, double thickness[])
 {
-  int i;
   double r = core_radius;
-  for (i=0; i < n; i++) {
+  for (int i=0; i < n; i++) {
     r += thickness[i];
   }
-  return 4.0/3.0 * M_PI * r * r * r;
+  return M_4PI_3 * cube(r);
 }
 
-static double
-Iq(double q, double core_sld, double core_radius, double solvent_sld,
-    double n, double in_sld[], double out_sld[], double thickness[],
-    double A[])
+double
+Iq(double q, double core_sld, double core_radius,
+   double solvent_sld, double num_shells, double sld[], double thickness[]);
+double
+Iq(double q, double core_sld, double core_radius,
+   double solvent_sld, double num_shells, double sld[], double thickness[])
 {
-  int i;
+  const int n = (int)ceil(num_shells);
+  double f, r, last_sld;
   r = core_radius;
-  f = f_constant(q, r, core_sld);
-  for (i=0; i<n; i++){
-    const double r0 = r;
+  last_sld = core_sld;
+  f = 0.;
+  for (int i=0; i<n; i++) {
+    f += M_4PI_3 * cube(r) * (sld[i] - last_sld) * sph_j1c(q*r);
+    last_sld = sld[i];
     r += thickness[i];
-    if (r == r0) {
-      // no thickness, so nothing to add
-    } else if (fabs(A[i]) < 1.0e-16 || sld_out[i] == sld_in[i]) {
-      f -= f_constant(q, r0, sld_in[i]);
-      f += f_constant(q, r, sld_in[i]);
-    } else if (fabs(A[i]) < 1.0e-4) {
-      const double slope = (sld_out[i] - sld_in[i])/thickness[i];
-      f -= f_linear(q, r0, sld_in[i], slope);
-      f += f_linear(q, r, sld_out[i], slope);
-    } else {
-      f -= f_exp(q, r0, sld_in[i], sld_out[i], thickness[i], A[i]);
-      f += f_exp(q, r, sld_in[i], sld_out[i], thickness[i], A[i]);
-    }
   }
-  f -= f_constant(q, r, solvent_sld);
-  f2 = f * f * 1.0e-4;
-
-  return f2;
+  f += M_4PI_3 * cube(r) * (solvent_sld - last_sld) * sph_j1c(q*r);
+  return f * f * 1.0e-4;
 }
Index: sasmodels/models/core_multi_shell.py
===================================================================
--- sasmodels/models/core_multi_shell.py	(revision 46ed76021e194fbd75cd7aeac6434e865025002f)
+++ sasmodels/models/core_multi_shell.py	(revision c5ac2b2355ae3705fd04e87490d99bfe4140d545)
@@ -87,7 +87,5 @@
 
 #             ["name", "units", default, [lower, upper], "type","description"],
-parameters = [["volfraction", "", 0.05, [0,1],"",
-               "volume fraction of spheres"],
-              ["sld_core", "1e-6/Ang^2", 1.0, [-inf, inf], "",
+parameters = [["sld_core", "1e-6/Ang^2", 1.0, [-inf, inf], "",
                "Core scattering length density"],
               ["radius", "Ang", 200., [0, inf], "volume",
@@ -95,5 +93,5 @@
               ["sld_solvent", "1e-6/Ang^2", 6.4, [-inf, inf], "",
                "Solvent scattering length density"],
-              ["n", "", 1, [0, 4], "volume",
+              ["n", "", 1, [0, 10], "volume",
                "number of shells"],
               ["sld[n]", "1e-6/Ang^2", 1.7, [-inf, inf], "",
@@ -103,14 +101,7 @@
               ]
 
-#source = ["lib/sph_j1c.c", "onion.c"]
+source = ["lib/sph_j1c.c", "core_multi_shell.c"]
 
-def Iq(q, *args, **kw):
-    return q
-
-def Iqxy(qx, *args, **kw):
-    return qx
-
-
-def profile(volfraction, sld_core, radius, sld_solvent, n, sld, thickness):
+def profile(sld_core, radius, sld_solvent, n, sld, thickness):
     """
     SLD profile
@@ -120,37 +111,4 @@
 
     """
-#        r = []
-#        beta = []
-#        # for core at r=0
-#        r.append(0)
-#        beta.append(self.params['sld_core0'])
-#        # for core at r=rad_core
-#        r.append(self.params['rad_core0'])
-#        beta.append(self.params['sld_core0'])
-#
-#        # for shells
-#        for n in range(1, self.n_shells+1):
-#            # Left side of each shells
-#            r0 = r[len(r)-1]
-#            r.append(r0)
-#            exec "beta.append(self.params['sld_shell%s'% str(n)])"
-#
-#            # Right side of each shells
-#            exec "r0 += self.params['thick_shell%s'% str(n)]"
-#            r.append(r0)
-#            exec "beta.append(self.params['sld_shell%s'% str(n)])"
-#
-#        # for solvent
-#        r0 = r[len(r)-1]            
-#        r.append(r0)
-#        beta.append(self.params['sld_solv'])
-#        r_solv = 5*r0/4
-#        r.append(r_solv)
-#        beta.append(self.params['sld_solv'])
-#
-#        return r, beta
-# above is directly from old code -- below is alotered from Woitek's first
-# cut an the onion.  Seems like we should be consistant?
-
     total_radius = 1.25*(sum(thickness[:n]) + radius + 1)
 
@@ -172,5 +130,5 @@
         r.append(r0 + thickness[k])
         beta.append(sld[k])
-   # add in the solvent
+    # add in the solvent
     r.append(r[-1])
     beta.append(sld_solvent)
@@ -187,9 +145,13 @@
     return 1.0, 1.0
 
-demo = dict(volfraction = 1.0,
-            sld_core = 6.4,
+demo = dict(sld_core = 6.4,
             radius = 60,
             sld_solvent = 6.4,
-            n = 1,
-            sld = [2.0],
-            thickness = [10])
+            n = 2,
+            sld = [2.0, 3.0],
+            thickness = 20,
+            thickness1_pd = 0.3,
+            thickness2_pd = 0.3,
+            thickness1_pd_n = 10,
+            thickness2_pd_n = 10,
+            )
