Index: sasmodels/models/spherical_sld.c
===================================================================
--- sasmodels/models/spherical_sld.c	(revision 1bf66d96606bebf9d19bae13c5509e74dddc0477)
+++ sasmodels/models/spherical_sld.c	(revision af0e70c2250b46a7401b4a45e22de4dfc5b2dc8f)
@@ -2,4 +2,5 @@
     int n_shells,
     double radius_core,
+    double thick_inter0,
     double thick_flat[],
     double thick_inter[])
@@ -7,4 +8,5 @@
     int i;
     double r = radius_core;
+    r += thick_inter0;
     for (i=0; i < n_shells; i++) {
         r += thick_inter[i];
@@ -15,152 +17,148 @@
 
 
-static double sphere_sld_kernel(double dp[], double q) {
-  int n = dp[0];
-  int i,j,k;
-
-  double scale = dp[1];
-  double thick_inter_core = dp[2];
-  double sld_core = dp[4];
-  double sld_solv = dp[5];
-  double background = dp[6];
-  double npts = dp[57]; //number of sub_layers in each interface
-  double nsl=npts;//21.0; //nsl = Num_sub_layer:  must be ODD double number
-  int n_s;
-
-  double sld_i,sld_f,dz,bes,fun,f,vol,qr,r,contr,f2;
-  double sign,slope=0.0;
-  double pi;
-
-  double total_thick=0.0;
-
-  int fun_type[12];
-  double sld[12];
-  double thick_inter[12];
-  double thick[12];
-  double fun_coef[12];
-
-  fun_type[0] = dp[3];
-  fun_coef[0] = fabs(dp[58]);
-  for (i =1; i<=n; i++){
-    sld[i] = dp[i+6];
-    thick_inter[i]= dp[i+16];
-    thick[i] = dp[i+26];
-    fun_type[i] = dp[i+36];
-    fun_coef[i] = fabs(dp[i+46]);
-    total_thick += thick[i];
-    total_thick += thick_inter[i];
-  }
-  sld[0] = sld_core;
-  sld[n+1] = sld_solv;
-  thick[0] = dp[59];
-  thick[n+1] = total_thick/5.0;
-  thick_inter[0] = thick_inter_core;
-  thick_inter[n+1] = 0.0;
-  fun_coef[n+1] = 0.0;
-  pi = 4.0*atan(1.0);
-  f = 0.0;
-  r = 0.0;
-  vol = 0.0;
-  //vol_pre = 0.0;
-  //vol_sub = 0.0;
-  sld_f = sld_core;
-
-  //floor_nsl = floor(nsl/2.0);
-
-  dz = 0.0;
-  // iteration for # of shells + core + solvent
-  for (i=0;i<=n+1; i++){
-    //iteration for N sub-layers
-    //if (fabs(thick[i]) <= 1e-24){
-    //   continue;
-    //}
-    // iteration for flat and interface
-    for (j=0;j<2;j++){
-      // iteration for sub_shells in the interface
-      // starts from #1 sub-layer
-      for (n_s=1;n_s<=nsl; n_s++){
-        // for solvent, it doesn't have an interface
-        if (i==n+1 && j==1)
-          break;
-        // for flat layers
-        if (j==0){
-          dz = thick[i];
-          sld_i = sld[i];
-          slope = 0.0;
+static double sphere_sld_kernel(
+    double q,
+    int n_shells,
+    int npts_inter,
+    double radius_core,
+    double sld_core,
+    double sld_solvent,
+    double func_inter_core,
+    double thick_inter_core,
+    double nu_inter_core,
+    double sld_flat[],
+    double thick_flat[],
+    double func_inter[],
+    double thick_inter[],
+    double nu_inter[] ) {
+
+    int i,j,k;
+    int n_s;
+
+    double sld_i,sld_f,dz,bes,fun,f,vol,qr,r,contr,f2;
+    double sign,slope=0.0;
+    double pi;
+
+    double total_thick=0.0;
+
+    //TODO: This part can be further simplified
+    int fun_type[12];
+    double sld[12];
+    double thick_internal[12];
+    double thick[12];
+    double fun_coef[12];
+
+    fun_type[0] = func_inter_core;
+    fun_coef[0] = fabs(nu_inter_core);
+    sld[0] = sld_core;
+    thick[0] = radius_core;
+    thick_internal[0] = thick_inter_core;
+
+    for (i =1; i<=n_shells; i++){
+        sld[i] = sld_flat[i-1];
+        thick_internal[i]= thick_inter[i-1];
+        thick[i] = thick_flat[i-1];
+        fun_type[i] = func_inter[i-1];
+        fun_coef[i] = fabs(nu_inter[i-1]);
+        total_thick += thick[i];
+        total_thick += thick_internal[i]; //doesn't account for core layer
+    }
+
+    sld[n_shells+1] = sld_solvent;
+    thick[n_shells+1] = total_thick/5.0;
+    thick_internal[n_shells+1] = 0.0;
+    fun_coef[n_shells+1] = 0.0;
+    fun_type[n_shells+1] = 0;
+
+    pi = 4.0*atan(1.0);
+    f = 0.0;
+    r = 0.0;
+    vol = 0.0;
+    sld_f = sld_core;
+
+    dz = 0.0;
+    // iteration for # of shells + core + solvent
+    for (i=0;i<=n_shells+1; i++){
+        // iteration for flat and interface
+        for (j=0;j<2;j++){
+            // iteration for sub_shells in the interface
+            // starts from #1 sub-layer
+            for (n_s=1;n_s<=npts_inter; n_s++){
+                // for solvent, it doesn't have an interface
+                if (i==n_shells+1 && j==1)
+                    break;
+                // for flat layers
+                if (j==0){
+                    dz = thick[i];
+                    sld_i = sld[i];
+                    slope = 0.0;
+                }
+                // for interfacial sub_shells
+                else{
+                    dz = thick_internal[i]/npts_inter;
+                    // find sld_i at the outer boundary of sub-layer #n_s
+                    sld_i = intersldfunc(fun_type[i], npts_inter, n_s,
+                            fun_coef[i], sld[i], sld[i+1]);
+                    // calculate slope
+                    slope= (sld_i -sld_f)/dz;
+                }
+                contr = sld_f-slope*r;
+                // iteration for the left and right boundary of the shells
+                for (k=0; k<2; k++){
+                    // At r=0, the contribution to I is zero, so skip it.
+                    if ( i == 0 && j == 0 && k == 0){
+                        continue;
+                    }
+                    // On the top of slovent there is no interface; skip it.
+                    if (i == n_shells+1 && k == 1){
+                        continue;
+                    }
+                    // At the right side (outer) boundary
+                    if ( k == 1){
+                        sign = 1.0;
+                        r += dz;
+                    }
+                    // At the left side (inner) boundary
+                    else{
+                        sign = -1.0;
+                    }
+
+                    qr = q * r;
+                    fun = 0.0;
+
+                    if(qr == 0.0){
+                        bes = sign * 1.0;
+                    }
+                    else{
+                        // for flat sub-layer
+                        //TODO: Single precision calculation fails here
+                        bes = sign *  sph_j1c(qr);
+                        if (fabs(slope) > 0.0 ){
+                            const double qrsq = qr*qr;
+                            double sinqr, cosqr;
+                            SINCOS(qr, sinqr, cosqr);
+                            fun = sign * 3.0 * r *
+                            (2.0*qr*sinqr - (qrsq-2.0)*cosqr)/(qrsq * qrsq);
+                            // In the onioon model Jae-He's formula is rederived
+                            // and gives following:
+                            //fun = 3.0 * sign * r *
+                            //(2.0*cosqr + qr*sinqr)/(qrsq*qrsq);
+                            //But this seems not to be working in this case...
+                        }
+                    }
+
+                    // update total volume
+                    vol = M_4PI_3 * cube(r);
+                    f += vol * (bes * contr + fun * slope);
+                }
+                sld_f = sld_i;
+                // no sub-layer iteration (n_s loop) for the flat layer
+                if (j==0)
+                    break;
+            }
         }
-        // for interfacial sub_shells
-        else{
-          dz = thick_inter[i]/nsl;
-          // find sld_i at the outer boundary of sub-layer #n_s
-          sld_i = intersldfunc(fun_type[i],nsl, n_s, fun_coef[i], sld[i], sld[i+1]);
-          // calculate slope
-          slope= (sld_i -sld_f)/dz;
-        }
-        contr = sld_f-slope*r;
-        // iteration for the left and right boundary of the shells(or sub_shells)
-        for (k=0; k<2; k++){
-          // At r=0, the contribution to I is zero, so skip it.
-          if ( i == 0 && j == 0 && k == 0){
-            continue;
-          }
-          // On the top of slovent there is no interface; skip it.
-          if (i == n+1 && k == 1){
-            continue;
-          }
-          // At the right side (outer) boundary
-          if ( k == 1){
-            sign = 1.0;
-            r += dz;
-          }
-          // At the left side (inner) boundary
-          else{
-            sign = -1.0;
-          }
-          qr = q * r;
-          fun = 0.0;
-
-          if(qr == 0.0){
-            // sigular point
-            bes = sign * 1.0;
-          }
-          else{
-            // for flat sub-layer
-            //TODO: Single precision calculation most likely fails here
-            //bes = sign *  3.0 * (sin(qr) - qr * cos(qr)) / (qr * qr * qr);
-            bes = sign *  sph_j1c(qr);
-            if (fabs(slope) > 0.0 ){
-              //fun = sign * 3.0 * r * (2.0*qr*sin(qr)-((qr*qr)-2.0)*cos(qr))/(qr * qr * qr * qr);
-              fun = sign * r * sph_j1c(qr) + sign * 3.0 * sin(qr)/(qr * qr * q )
-                + sign * 6.0 * cos(qr)/(qr * qr * qr * q);
-            }
-          }
-
-          //Some initial optimization tries
-          /*bes = (qr == 0.0 ? sign * 1.0 : sign *  3.0 * (sin(qr) - qr * cos(qr)) / (qr * qr * qr));
-          //TODO: Will have to chnage this function
-          if (qr!= 0.0 && fabs(slope) > 0.0 ){
-            fun = sign * 3.0 * r * (2.0*qr*sin(qr)-((qr*qr)-2.0)*cos(qr))/(qr * qr * qr * qr);
-          }*/
-
-          // update total volume
-          vol = 4.0 * pi / 3.0 * r * r * r;
-          // we won't do the following volume correction for now.
-          // substrate empty area of volume
-          //if (k == 1 && fabs(sld_in[i]-sld_solv) < 1e-04*fabs(sld_solv) && fun_type[i]==0){
-          //  vol_sub += (vol_pre - vol);
-          //}
-          f += vol * (bes * contr + fun * slope);
-        }
-        // remember this sld as sld_f
-        sld_f = sld_i;
-        // no sub-layer iteration (n_s loop) for the flat layer
-        if (j==0)
-          break;
-      }
-    }
-  }
-  f2 = f * f / vol;
-
-  return (f2);
+    }
+    f2 = f * f / vol;
+    return (f2);
 }
 
@@ -177,4 +175,7 @@
     double sld_core,
     double sld_solvent,
+    double func_inter0,
+    double thick_inter0,
+    double nu_inter0,
     double sld_flat[],
     double thick_flat[],
@@ -186,21 +187,15 @@
     double intensity;
     //TODO: Remove this container at later stage.
-    double dp[60];
+    /*double dp[60];
     dp[0] = n_shells;
     //This is scale will also have to be removed at some stage
     dp[1] = 1.0;
-    dp[2] = thick_inter[0];
-    dp[3] = func_inter[0];
+    dp[2] = thick_inter0;
+    dp[3] = func_inter0;
     dp[4] = sld_core;
     dp[5] = sld_solvent;
     dp[6] = 0.0;
-    dp[7] = sld_flat[0];
-    //TODO: Something is messed up with this data strcucture!
-    dp[17] = thick_inter[0];
-    dp[27] = thick_flat[0];
-    dp[37] = func_inter[0];
-    dp[47] = nu_inter[0];
-
-    for (int i=1; i<=n_shells; i++){
+
+    for (int i=0; i<n_shells; i++){
         dp[i+7] = sld_flat[i];
         dp[i+17] = thick_inter[i];
@@ -211,8 +206,11 @@
 
     dp[57] = npts_inter;
-    dp[58] = nu_inter[0];
+    dp[58] = nu_inter0;
     dp[59] = radius_core;
-
-    intensity = 1.0e-4*sphere_sld_kernel(dp,q);
+    */
+    intensity = sphere_sld_kernel(q, n_shells, npts_inter, radius_core,
+                sld_core, sld_solvent, func_inter0, thick_inter0, nu_inter0,
+                sld_flat, thick_flat, func_inter, thick_inter, nu_inter);
+    intensity *=1.0e-4;
     //printf("%10d\n",intensity);
     return intensity;
Index: sasmodels/models/spherical_sld.py
===================================================================
--- sasmodels/models/spherical_sld.py	(revision 1448db7a455dc72982dab9c6602b6989ff8bfc65)
+++ sasmodels/models/spherical_sld.py	(revision af0e70c2250b46a7401b4a45e22de4dfc5b2dc8f)
@@ -1,11 +1,15 @@
 r"""
-This model calculates an empirical functional form for SAS data using SpericalSLD profile
-
-Similarly to the OnionExpShellModel, this model provides the form factor, P(q), for a multi-shell sphere,
-where the interface between the each neighboring shells can be described by one of a number of functions
-including error, power-law, and exponential functions. This model is to calculate the scattering intensity
-by building a continuous custom SLD profile against the radius of the particle. The SLD profile is composed
-of a flat core, a flat solvent, a number (up to 9 ) flat shells, and the interfacial layers between
-the adjacent flat shells (or core, and solvent) (see below).
+This model calculates an empirical functional form for SAS data using
+SpericalSLD profile
+
+Similarly to the OnionExpShellModel, this model provides the form factor,
+P(q), for a multi-shell sphere, where the interface between the each neighboring
+shells can be described by one of a number of functions including error,
+power-law, and exponential functions.
+This model is to calculate the scattering intensity by building a continuous
+custom SLD profile against the radius of the particle.
+The SLD profile is composed of a flat core, a flat solvent, a number (up to 9 )
+flat shells, and the interfacial layers between the adjacent flat shells
+(or core, and solvent) (see below).
 
 .. figure:: img/spherical_sld_profile.gif
@@ -13,8 +17,9 @@
     Exemplary SLD profile
 
-Unlike the <onion> model (using an analytical integration),
-the interfacial layers here are sub-divided and numerically integrated assuming each of the sub-layers are described
-by a line function. The number of the sub-layer can be given by users by setting the integer values of npts_inter.
-The form factor is normalized by the total volume of the sphere.
+Unlike the <onion> model (using an analytical integration), the interfacial
+layers here are sub-divided and numerically integrated assuming each of the
+sub-layers are described by a line function.
+The number of the sub-layer can be given by users by setting the integer values
+of npts_inter. The form factor is normalized by the total volume of the sphere.
 
 Definition
@@ -29,5 +34,6 @@
     \sum_{\text{flat}_i=0}^N f_{\text{flat}_i} +f_\text{solvent}
 
-For a spherically symmetric particle with a particle density $\rho_x(r)$ the sld function can be defined as:
+For a spherically symmetric particle with a particle density $\rho_x(r)$
+the sld function can be defined as:
 
 .. math::
@@ -39,19 +45,27 @@
 
 .. math::
-    f_\text{core} = 4 \pi \int_{0}^{r_\text{core}} \rho_\text{core} \frac{\sin(qr)} {qr} r^2 dr =
+    f_\text{core} = 4 \pi \int_{0}^{r_\text{core}} \rho_\text{core}
+    \frac{\sin(qr)} {qr} r^2 dr =
     3 \rho_\text{core} V(r_\text{core})
-    \Big[ \frac{\sin(qr_\text{core}) - qr_\text{core} \cos(qr_\text{core})} {qr_\text{core}^3} \Big]
-
-    f_{\text{inter}_i} = 4 \pi \int_{\Delta t_{ \text{inter}_i } } \rho_{ \text{inter}_i } \frac{\sin(qr)} {qr} r^2 dr
-
-    f_{\text{shell}_i} = 4 \pi \int_{\Delta t_{ \text{inter}_i } } \rho_{ \text{flat}_i } \frac{\sin(qr)} {qr} r^2 dr =
-    3 \rho_{ \text{flat}_i } V ( r_{ \text{inter}_i } + \Delta t_{ \text{inter}_i } )
-    \Big[ \frac{\sin(qr_{\text{inter}_i} + \Delta t_{ \text{inter}_i } ) - q (r_{\text{inter}_i} +
-    \Delta t_{ \text{inter}_i }) \cos(q( r_{\text{inter}_i} + \Delta t_{ \text{inter}_i } ) ) }
+    \Big[ \frac{\sin(qr_\text{core}) - qr_\text{core} \cos(qr_\text{core})}
+    {qr_\text{core}^3} \Big]
+
+    f_{\text{inter}_i} = 4 \pi \int_{\Delta t_{ \text{inter}_i } }
+    \rho_{ \text{inter}_i } \frac{\sin(qr)} {qr} r^2 dr
+
+    f_{\text{shell}_i} = 4 \pi \int_{\Delta t_{ \text{inter}_i } }
+    \rho_{ \text{flat}_i } \frac{\sin(qr)} {qr} r^2 dr =
+    3 \rho_{ \text{flat}_i } V ( r_{ \text{inter}_i } +
+    \Delta t_{ \text{inter}_i } )
+    \Big[ \frac{\sin(qr_{\text{inter}_i} + \Delta t_{ \text{inter}_i } )
+    - q (r_{\text{inter}_i} + \Delta t_{ \text{inter}_i })
+    \cos(q( r_{\text{inter}_i} + \Delta t_{ \text{inter}_i } ) ) }
     {q ( r_{\text{inter}_i} + \Delta t_{ \text{inter}_i } )^3 }  \Big]
     -3 \rho_{ \text{flat}_i } V(r_{ \text{inter}_i })
-    \Big[ \frac{\sin(qr_{\text{inter}_i}) - qr_{\text{flat}_i} \cos(qr_{\text{inter}_i}) } {qr_{\text{inter}_i}^3} \Big]
-
-    f_\text{solvent} = 4 \pi \int_{r_N}^{\infty} \rho_\text{solvent} \frac{\sin(qr)} {qr} r^2 dr =
+    \Big[ \frac{\sin(qr_{\text{inter}_i}) - qr_{\text{flat}_i}
+    \cos(qr_{\text{inter}_i}) } {qr_{\text{inter}_i}^3} \Big]
+
+    f_\text{solvent} = 4 \pi \int_{r_N}^{\infty} \rho_\text{solvent}
+    \frac{\sin(qr)} {qr} r^2 dr =
     3 \rho_\text{solvent} V(r_N)
     \Big[ \frac{\sin(qr_N) - qr_N \cos(qr_N)} {qr_N^3} \Big]
@@ -66,6 +80,8 @@
 .. math::
     \rho_{{inter}_i} (r) = \begin{cases}
-    B \exp\Big( \frac {\pm A(r - r_{\text{flat}_i})} {\Delta t_{ \text{inter}_i }} \Big) +C  & \text{for} A \neq 0 \\
-    B \Big( \frac {(r - r_{\text{flat}_i})} {\Delta t_{ \text{inter}_i }} \Big) +C  & \text{for} A = 0 \\
+    B \exp\Big( \frac {\pm A(r - r_{\text{flat}_i})}
+    {\Delta t_{ \text{inter}_i }} \Big) +C  & \text{for} A \neq 0 \\
+    B \Big( \frac {(r - r_{\text{flat}_i})}
+    {\Delta t_{ \text{inter}_i }} \Big) +C  & \text{for} A = 0 \\
     \end{cases}
 
@@ -74,5 +90,6 @@
 .. math::
     \rho_{{inter}_i} (r) = \begin{cases}
-    \pm B \Big( \frac {(r - r_{\text{flat}_i} )} {\Delta t_{ \text{inter}_i }} \Big) ^A  +C  & \text{for} A \neq 0 \\
+    \pm B \Big( \frac {(r - r_{\text{flat}_i} )} {\Delta t_{ \text{inter}_i }}
+    \Big) ^A  +C  & \text{for} A \neq 0 \\
     \rho_{\text{flat}_{i+1}}  & \text{for} A = 0 \\
     \end{cases}
@@ -82,26 +99,35 @@
 .. math::
     \rho_{{inter}_i} (r) = \begin{cases}
-    B \text{erf} \Big( \frac { A(r - r_{\text{flat}_i})} {\sqrt{2} \Delta t_{ \text{inter}_i }} \Big) +C  & \text{for} A \neq 0 \\
-    B \Big( \frac {(r - r_{\text{flat}_i} )} {\Delta t_{ \text{inter}_i }} \Big)  +C  & \text{for} A = 0 \\
+    B \text{erf} \Big( \frac { A(r - r_{\text{flat}_i})}
+    {\sqrt{2} \Delta t_{ \text{inter}_i }} \Big) +C  & \text{for} A \neq 0 \\
+    B \Big( \frac {(r - r_{\text{flat}_i} )} {\Delta t_{ \text{inter}_i }}
+    \Big)  +C  & \text{for} A = 0 \\
     \end{cases}
 
-The functions are normalized so that they vary between 0 and 1, and they are constrained such that the SLD
-is continuous at the boundaries of the interface as well as each sub-layers. Thus B and C are determined.
-
-Once $\rho_{\text{inter}_i}$ is found at the boundary of the sub-layer of the interface, we can find its contribution
-to the form factor $P(q)$
-
-.. math::
-    f_{\text{inter}_i} = 4 \pi \int_{\Delta t_{ \text{inter}_i } } \rho_{ \text{inter}_i } \frac{\sin(qr)} {qr} r^2 dr =
+The functions are normalized so that they vary between 0 and 1, and they are
+constrained such that the SLD is continuous at the boundaries of the interface
+as well as each sub-layers. Thus B and C are determined.
+
+Once $\rho_{\text{inter}_i}$ is found at the boundary of the sub-layer of the
+interface, we can find its contribution to the form factor $P(q)$
+
+.. math::
+    f_{\text{inter}_i} = 4 \pi \int_{\Delta t_{ \text{inter}_i } }
+    \rho_{ \text{inter}_i } \frac{\sin(qr)} {qr} r^2 dr =
     4 \pi \sum_{j=0}^{npts_{\text{inter}_i} -1 }
-    \int_{r_j}^{r_{j+1}} \rho_{ \text{inter}_i } (r_j) \frac{\sin(qr)} {qr} r^2 dr \approx
+    \int_{r_j}^{r_{j+1}} \rho_{ \text{inter}_i } (r_j)
+    \frac{\sin(qr)} {qr} r^2 dr \approx
 
     4 \pi \sum_{j=0}^{npts_{\text{inter}_i} -1 } \Big[
-    3 ( \rho_{ \text{inter}_i } ( r_{j+1} ) - \rho_{ \text{inter}_i } ( r_{j} ) V ( r_{ \text{sublayer}_j } )
-    \Big[ \frac {r_j^2 \beta_\text{out}^2 \sin(\beta_\text{out}) - (\beta_\text{out}^2-2) \cos(\beta_\text{out}) }
+    3 ( \rho_{ \text{inter}_i } ( r_{j+1} ) - \rho_{ \text{inter}_i }
+    ( r_{j} ) V ( r_{ \text{sublayer}_j } )
+    \Big[ \frac {r_j^2 \beta_\text{out}^2 \sin(\beta_\text{out})
+    - (\beta_\text{out}^2-2) \cos(\beta_\text{out}) }
     {\beta_\text{out}^4 } \Big]
 
-    - 3 ( \rho_{ \text{inter}_i } ( r_{j+1} ) - \rho_{ \text{inter}_i } ( r_{j} ) V ( r_{ \text{sublayer}_j-1 } )
-    \Big[ \frac {r_{j-1}^2 \sin(\beta_\text{in}) - (\beta_\text{in}^2-2) \cos(\beta_\text{in}) }
+    - 3 ( \rho_{ \text{inter}_i } ( r_{j+1} ) - \rho_{ \text{inter}_i }
+    ( r_{j} ) V ( r_{ \text{sublayer}_j-1 } )
+    \Big[ \frac {r_{j-1}^2 \sin(\beta_\text{in})
+    - (\beta_\text{in}^2-2) \cos(\beta_\text{in}) }
     {\beta_\text{in}^4 } \Big]
 
@@ -120,10 +146,12 @@
     V(a) = \frac {4\pi}{3}a^3
 
-    a_\text{in} ~ \frac{r_j}{r_{j+1} -r_j} \text{, } a_\text{out} ~ \frac{r_{j+1}}{r_{j+1} -r_j}
+    a_\text{in} ~ \frac{r_j}{r_{j+1} -r_j} \text{, } a_\text{out}
+    ~ \frac{r_{j+1}}{r_{j+1} -r_j}
 
     \beta_\text{in} = qr_j \text{, } \beta_\text{out} = qr_{j+1}
 
 
-We assume the $\rho_{\text{inter}_i} (r)$ can be approximately linear within a sub-layer $j$
+We assume the $\rho_{\text{inter}_i} (r)$ can be approximately linear
+within a sub-layer $j$
 
 Finally form factor can be calculated by
@@ -131,5 +159,6 @@
 .. math::
 
-    P(q) = \frac{[f]^2} {V_\text{particle}} \text{where} V_\text{particle} = V(r_{\text{shell}_N})
+    P(q) = \frac{[f]^2} {V_\text{particle}} \text{where} V_\text{particle}
+    = V(r_{\text{shell}_N})
 
 For 2D data the scattering intensity is calculated in the same way as 1D,
@@ -150,9 +179,11 @@
 
 .. note::
-    The outer most radius is used as the effective radius for S(Q) when $P(Q) * S(Q)$ is applied.
+    The outer most radius is used as the effective radius for S(Q)
+    when $P(Q) * S(Q)$ is applied.
 
 References
 ----------
-L A Feigin and D I Svergun, Structure Analysis by Small-Angle X-Ray and Neutron Scattering, Plenum Press, New York, (1987)
+L A Feigin and D I Svergun, Structure Analysis by Small-Angle X-Ray
+and Neutron Scattering, Plenum Press, New York, (1987)
 
 """
@@ -170,13 +201,16 @@
 # pylint: disable=bad-whitespace, line-too-long
 #            ["name", "units", default, [lower, upper], "type", "description"],
-parameters = [["n_shells",                "",               1,      [0, 9],         "", "number of shells"],
-              ["npts_inter",       "",               35,     [0, 35],        "", "number of points in each sublayer Must be odd number"],
-              ["radius_core",      "Ang",            50.0,   [0, inf],       "", "intern layer thickness"],
-              ["sld_core",         "1e-6/Ang^2",     2.07,   [-inf, inf],    "", "sld function flat"],
-              ["sld_solvent",      "1e-6/Ang^2",     1.0,    [-inf, inf],    "", "sld function solvent"],
+parameters = [["n_shells",          "",               1,      [0, 9],         "volume", "number of shells"],
+              ["npts_inter",        "",               35,     [0, inf],        "", "number of points in each sublayer Must be odd number"],
+              ["radius_core",       "Ang",            50.0,   [0, inf],       "volume", "intern layer thickness"],
+              ["sld_core",          "1e-6/Ang^2",     2.07,   [-inf, inf],    "", "sld function flat"],
+              ["sld_solvent",       "1e-6/Ang^2",     1.0,    [-inf, inf],    "", "sld function solvent"],
+              ["func_inter0",       "",               0,      [0, 4],         "", "Erf:0, RPower:1, LPower:2, RExp:3, LExp:4"],
+              ["thick_inter0",      "Ang",            50.0,   [0, inf],       "volume", "intern layer thickness for core layer"],
+              ["nu_inter0",         "",               2.5,    [-inf, inf],    "", "steepness parameter for core layer"],
               ["sld_flat[n_shells]",      "1e-6/Ang^2",     4.06,   [-inf, inf],    "", "sld function flat"],
-              ["thick_flat[n_shells]",    "Ang",            100.0,  [0, inf],       "", "flat layer_thickness"],
+              ["thick_flat[n_shells]",    "Ang",            100.0,  [0, inf],       "volume", "flat layer_thickness"],
               ["func_inter[n_shells]",    "",               0,      [0, 4],         "", "Erf:0, RPower:1, LPower:2, RExp:3, LExp:4"],
-              ["thick_inter[n_shells]",   "Ang",            50.0,   [0, inf],       "", "intern layer thickness"],
+              ["thick_inter[n_shells]",   "Ang",            50.0,   [0, inf],       "volume", "intern layer thickness"],
               ["nu_inter[n_shells]",      "",               2.5,    [-inf, inf],    "", "steepness parameter"],
               ]
@@ -184,4 +218,5 @@
 source = ["lib/librefl.c",  "lib/sph_j1c.c", "spherical_sld.c"]
 
+profile_axes = ['Radius (A)', 'SLD (1e-6/A^2)']
 def profile(n_shells, radius_core,  sld_core,  sld_solvent, sld_flat,
             thick_flat, func_inter, thick_inter, nu_inter, npts_inter):
@@ -256,14 +291,17 @@
 
 demo = {
-    "n_shells":4,
-    "npts_inter":35.0,
-    "radius_core":50.0,
-    "sld_core":2.07,
+    "n_shells": 4,
+    "npts_inter": 35.0,
+    "radius_core": 50.0,
+    "sld_core": 2.07,
     "sld_solvent": 1.0,
-    "sld_flat":[4.0,3.5,4.0,3.5,4.0],
-    "thick_flat":[100.0,100.0,100.0,100.0,100.0],
-    "func_inter":[0,0,0,0,0],
-    "thick_inter":[50.0,50.0,50.0,50.0,50.0],
-    "nu_inter":[2.5,2.5,2.5,2.5,2.5]
+    "thick_inter0": 50.0,
+    "func_inter0": 0,
+    "nu_inter0": 2.5,
+    "sld_flat":[4.0,3.5,4.0,3.5],
+    "thick_flat":[100.0,100.0,100.0,100.0],
+    "func_inter":[0,0,0,0],
+    "thick_inter":[50.0,50.0,50.0,50.0],
+    "nu_inter":[2.5,2.5,2.5,2.5],
     }
 
@@ -271,17 +309,14 @@
 tests = [
     # Accuracy tests based on content in test/utest_extra_models.py
-    [{'npts_iter':35,
-        'sld_solv':1,
-        'radius_core':50.0,
-        'sld_core':2.07,
-        'func_inter2':0.0,
-        'thick_inter2':50,
-        'nu_inter2':2.5,
-        'sld_flat2':4,
-        'thick_flat2':100,
-        'func_inter1':0.0,
-        'thick_inter1':50,
-        'nu_inter1':2.5,
-        'background': 0.0,
+    [{"n_shells":4,
+        'npts_inter':35,
+        "radius_core":50.0,
+        "sld_core":2.07,
+        "sld_solvent": 1.0,
+        "sld_flat":[4.0,3.5,4.0,3.5],
+        "thick_flat":[100.0,100.0,100.0,100.0],
+        "func_inter":[0,0,0,0],
+        "thick_inter":[50.0,50.0,50.0,50.0],
+        "nu_inter":[2.5,2.5,2.5,2.5]
     }, 0.001, 0.001],
 ]
