Index: sasmodels/models/core_multi_shell.c
===================================================================
--- sasmodels/models/core_multi_shell.c	(revision f7930be4820f6e24bcbf1db3df75a77f3f369c7a)
+++ sasmodels/models/core_multi_shell.c	(revision f7930be4820f6e24bcbf1db3df75a77f3f369c7a)
@@ -0,0 +1,200 @@
+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
+f_constant(double q, double r, double sld)
+{
+  const double bes = sph_j1c(q * r);
+  const double vol = 4.0/3.0 * M_PI * r * r * r;
+  return sld * vol * bes;
+}
+
+static double
+form_volume(double core_radius, double n, double thickness[])
+{
+  int i;
+  double r = core_radius;
+  for (i=0; i < n; i++) {
+    r += thickness[i];
+  }
+  return 4.0/3.0 * M_PI * r * r * 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[])
+{
+  int i;
+  r = core_radius;
+  f = f_constant(q, r, core_sld);
+  for (i=0; i<n; i++){
+    const double r0 = r;
+    r += thickness[i];
+    if (r == r0) {
+      // no thickness, so nothing to add
+    } else if (fabs(A[i]) < 1e-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]) < 1e-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;
+}
Index: sasmodels/models/core_multi_shell.py
===================================================================
--- sasmodels/models/core_multi_shell.py	(revision f7930be4820f6e24bcbf1db3df75a77f3f369c7a)
+++ sasmodels/models/core_multi_shell.py	(revision f7930be4820f6e24bcbf1db3df75a77f3f369c7a)
@@ -0,0 +1,221 @@
+r"""
+Definition
+----------
+
+This model is a trivial extension of the CoreShell function to a larger number
+of shells. The scattering length density profile for the default sld values 
+(w/ 4 shells).
+
+.. figure:: img/core_multi_shell_sld_default_profile.jpg
+
+    SLD profile of the core_multi_shell object from the center of sphere out
+    for the default SLDs.*
+
+The 2D scattering intensity is the same as *P(q)* above, regardless of the
+orientation of the *q* vector which is defined as
+
+.. math::
+
+    q = \sqrt{q_x^2 + q_y^2}
+
+.. note:: **Be careful!** The SLDs and scale can be highly correlated. Hold as
+         many of these parameters fixed as possible.
+
+.. note:: The outer most radius (= *radius* + *thickness*) is used as the
+          effective radius for *S(Q)* when *P(Q)* \* *S(Q)* is applied.
+
+For information about polarised and magnetic scattering, click here_.
+
+Our model uses the form factor calculations implemented in a c-library provided
+by the NIST Center for Neutron Research (Kline, 2006).
+
+References
+----------
+See the CoreShellModel_ documentation.
+
+L A Feigin and D I Svergun,
+*Structure Analysis by Small-Angle X-Ray and Neutron Scattering*,
+Plenum Press, New York, 1987.
+
+**Author:** NIST IGOR/DANSE **on:** pre 2010
+
+**Last Modified by:** in progress **on:** March 20, 2016
+
+**Last Reviewed by:** in progress **on:** March 20, 2016
+"""
+
+
+
+from __future__ import division
+
+import numpy as np
+from numpy import inf, nan
+from math import fabs, exp, expm1
+
+name = "core_multi_shell"
+title = "This model provides the scattering from a spherical core with 1 to 4 \
+ concentric shell structures. The SLDs of the core and each shell are \
+ individually specified."
+
+description = """\
+Form factor for a core muti-shell (up to 4) sphere normalized by the volume.
+Each shell can have a unique thickness and sld.
+
+	background:background,
+	rad_core0: radius of sphere(core)
+	thick_shell#:the thickness of the shell#
+	sld_core0: the SLD of the sphere
+	sld_solv: the SLD of the solvent
+	sld_shell: the SLD of the shell#
+	A_shell#: the coefficient in the exponential function
+	
+	
+    scale: 1.0 if data is on absolute scale
+    volfraction: volume fraction of spheres
+    radius: the radius of the core
+    sld: the SLD of the core
+    thick_shelli: the thickness of the i'th shell from the core
+    sld_shelli: the SLD of the i'th shell from the core
+    sld_solvent: the SLD of the solvent
+    background: incoherent background
+
+"""
+
+category = "shape:sphere"
+
+
+#             ["name", "units", default, [lower, upper], "type","description"],
+parameters = [["volfraction", "", 0.05, [0,1],"",
+               "volume fraction of spheres"],
+              ["sld", "1e-6/Ang^2", 1.0, [-inf, inf], "",
+               "Core scattering length density"],
+              ["radius", "Ang", 200., [0, inf], "",
+               "Radius of the core"],
+              ["sld_solvent", "1e-6/Ang^2", 6.4, [-inf, inf], "",
+               "Solvent scattering length density"],
+              ["n", "", 1, [0, 10], "volume",
+               "number of shells"],
+              ["sld_shell[n]", "1e-6/Ang^2", 1.7, [-inf, inf], "",
+               "scattering length density of shell k"],
+              ["thick_shell[n]", "Ang", 40., [0, inf], "volume",
+               "Thickness of shell k"],
+              ]
+
+#source = ["lib/sph_j1c.c", "onion.c"]
+
+def Iq(q, *args, **kw):
+    return q
+
+def Iqxy(qx, *args, **kw):
+    return qx
+
+
+def shape(core_sld, core_radius, solvent_sld, n, in_sld, out_sld, thickness, A):
+    """
+    SLD profile
+    
+    :return: (r, beta) where r is a list of radius of the transition points\
+      and beta is a list of the corresponding SLD values.
+
+    """
+#        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]) + core_radius + 1)
+    dr = total_radius/400  # 400 points for a smooth plot
+
+    r = []
+    beta = []
+
+    # add in the core
+    r.append(0)
+    beta.append(sld)
+    r.append(radius)
+    beta.append(sld)
+
+    # add in the shells
+    for k in range(n):
+        # Left side of each shells
+        r0 = r[-1]
+        r.append(r0)
+        beta.append(sld_shell[k])
+        r.append(r0 + thickness[k])
+        beta.append(sld_shell[k])
+   # add in the solvent
+    r.append(r[-1])
+    beta.append(solvent_sld)
+    r.append(total_radius)
+    beta.append(solvent_sld)
+
+    return np.asarray(r), np.asarray(beta)
+
+def ER(radius, n, thick_shell):
+    return np.sum(thick_shell[:n], axis=0) + core_radius
+
+def VR(radius, n, thick_shell):
+    return 1.0
+
+parameters = [["volfraction", "", 0.05, [0,1],"",
+               "volume fraction of spheres"],
+              ["sld", "1e-6/Ang^2", 1.0, [-inf, inf], "",
+               "Core scattering length density"],
+              ["radius", "Ang", 200., [0, inf], "",
+               "Radius of the core"],
+              ["sld_solvent", "1e-6/Ang^2", 6.4, [-inf, inf], "",
+               "Solvent scattering length density"],
+              ["n", "", 1, [0, 10], "volume",
+               "number of shells"],
+              ["sld_shell[n]", "1e-6/Ang^2", 1.7, [-inf, inf], "",
+               "scattering length density of shell k"],
+              ["thick_shell[n]", "Ang", 40., [0, inf], "volume",
+               "Thickness of shell k"],
+               ]
+
+demo = dict(volfraction = 1.0,
+            sld = 6.4,
+            radius = 60,
+            sld_solvent = 6.4,
+            n = 1,
+            sld_shell = [2.0],
+            thick_shell = [10])
+
+oldname = "CoreMultiShellModel"
+oldpars = dict(
+    sld="sld_core0",
+    radius="rad_core0",
+    sld_solvent="sld_solv",
+    n="n_shells",
+    sld_shell="sld_in_shell",
+    thick_shell="thick_shell",
+    # func_shell is always 2 in the user interface, so isn't included
+    )
