Index: doc/genmodel.py
===================================================================
--- doc/genmodel.py	(revision a5b8477c53695e7df03ecec7346f0c1bd770d3f7)
+++ doc/genmodel.py	(revision 745b7bb20cb177250d7a71a24058de0dcca609e4)
@@ -6,5 +6,5 @@
 sys.path.insert(0, os.path.abspath('..'))
 from sasmodels import generate, core
-from sasmodels.direct_model import DirectModel
+from sasmodels.direct_model import DirectModel, call_profile
 from sasmodels.data import empty_data1D, empty_data2D
 
@@ -17,4 +17,5 @@
     from sasmodels.kernel import KernelModel
     from sasmodels.modelinfo import ModelInfo
+
 
 def plot_1d(model, opts, ax):
@@ -55,4 +56,19 @@
     ax.set_xlabel(r'$Q_x \/(\AA^{-1})$')
     ax.set_ylabel(r'$Q_y \/(\AA^{-1})$')
+
+def plot_profile_inset(model_info, ax):
+    p = ax.get_position()
+    width, height = 0.4*(p.x1-p.x0), 0.4*(p.y1-p.y0)
+    left, bottom = p.x1-width, p.y1-height
+    inset = plt.gcf().add_axes([left, bottom, width, height])
+    x, y, labels = call_profile(model_info)
+    inset.plot(x, y, '-')
+    inset.locator_params(nbins=4)
+    #inset.set_xlabel(labels[0])
+    #inset.set_ylabel(labels[1])
+    inset.text(0.99, 0.99, "profile",
+               horizontalalignment="right",
+               verticalalignment="top",
+               transform=inset.transAxes)
 
 def figfile(model_info):
@@ -100,4 +116,7 @@
         ax1d = fig.add_axes([ax_left, ax_bottom, ax_width, ax_height])
         plot_1d(model, opts, ax1d)
+
+    if model_info.profile:
+        plot_profile_inset(model_info, ax1d)
 
     # Save image in model/img
Index: sasmodels/direct_model.py
===================================================================
--- sasmodels/direct_model.py	(revision 9eb36329a5763c7ea05bd24c85b4810c08c2c2e2)
+++ sasmodels/direct_model.py	(revision 745b7bb20cb177250d7a71a24058de0dcca609e4)
@@ -73,4 +73,5 @@
     #print("values:", values)
     return calculator(call_details, values, cutoff, is_magnetic)
+
 
 def get_weights(parameter, values):
@@ -95,4 +96,18 @@
         disperser, npts, width, nsigma, value, limits, relative)
     return value, weight / np.sum(weight)
+
+
+def call_profile(model_info, **pars):
+    args = {}
+    for p in model_info.parameters.kernel_parameters:
+        if p.length > 1:
+            value = np.array([pars.get(p.id+str(j), p.default)
+                              for j in range(1, p.length+1)])
+        else:
+            value = pars.get(p.id, p.default)
+        args[p.id] = value
+    x, y = model_info.profile(**args)
+    return x, y, model_info.profile_axes
+
 
 class DataMixin(object):
@@ -290,4 +305,11 @@
         self._set_data(Iq, noise=noise)
 
+    def profile(self, **pars):
+        # type: (**float) -> None
+        """
+        Generate a plottable profile.
+        """
+        return call_profile(self.model.info, **pars)
+
 def main():
     # type: () -> None
Index: smodels/models/lib/librefl.c
===================================================================
--- sasmodels/models/lib/librefl.c	(revision ee5d14de635651e6d83bc978943c11ee44a6b250)
+++ 	(revision )
@@ -1,202 +1,0 @@
-double err_mod_func(double n_sub, double ind, double nu);
-double linearfunc(double n_sub, double ind, double nu);
-double power_r(double n_sub, double ind, double nu);
-double power_l(double n_sub, double ind, double nu);
-double exp_r(double n_sub, double ind, double nu);
-double exp_l(double n_sub, double ind, double nu);
-double intersldfunc(int fun_type, double n_sub, double i, double nu, double sld_l, double sld_r);
-double interfunc(int fun_type, double n_sub, double i, double sld_l, double sld_r);
-
-// normalized and modified erf
-//   |
-// 1 +                __  - - - -
-//   |             _
-//	 |            _
-//   |        __
-// 0 + - - -
-//   |-------------+------------+--
-//   0           center       n_sub    --->
-//                                     ind
-//
-// n_sub = total no. of bins(or sublayers)
-// ind = x position: 0 to max
-// nu = max x to integration
-double err_mod_func(double n_sub, double ind, double nu)
-{
-  double center, func;
-  if (nu == 0.0)
-		nu = 1.0e-14;
-	if (n_sub == 0.0)
-		n_sub = 1.0;
-
-
-	//ind = (n_sub-1.0)/2.0-1.0 +ind;
-	center = n_sub/2.0;
-	// transform it so that min(ind) = 0
-	ind -= center;
-	// normalize by max limit
-	ind /= center;
-	// divide by sqrt(2) to get Gaussian func
-	nu /= sqrt(2.0);
-	ind *= nu;
-	// re-scale and normalize it so that max(erf)=1, min(erf)=0
-	func = sas_erf(ind)/sas_erf(nu)/2.0;
-	// shift it by +0.5 in y-direction so that min(erf) = 0
-	func += 0.5;
-
-	return func;
-}
-double linearfunc(double n_sub, double ind, double nu)
-{
-  double bin_size, func;
-	if (n_sub == 0.0)
-		n_sub = 1.0;
-
-	bin_size = 1.0/n_sub;  //size of each sub-layer
-	// rescale
-	ind *= bin_size;
-	func = ind;
-
-	return func;
-}
-// use the right hand side from the center of power func
-double power_r(double n_sub, double ind, double nu)
-{
-  double bin_size,func;
-	if (nu == 0.0)
-		nu = 1.0e-14;
-	if (n_sub == 0.0)
-		n_sub = 1.0;
-
-	bin_size = 1.0/n_sub;  //size of each sub-layer
-	// rescale
-	ind *= bin_size;
-	func = pow(ind, nu);
-
-	return func;
-}
-// use the left hand side from the center of power func
-double power_l(double n_sub, double ind, double nu)
-{
-  double bin_size, func;
-	if (nu == 0.0)
-		nu = 1.0e-14;
-	if (n_sub == 0.0)
-		n_sub = 1.0;
-
-	bin_size = 1.0/n_sub;  //size of each sub-layer
-	// rescale
-	ind *= bin_size;
-	func = 1.0-pow((1.0-ind),nu);
-
-	return func;
-}
-// use 1-exp func from x=0 to x=1
-double exp_r(double n_sub, double ind, double nu)
-{
-  double bin_size, func;
-	if (nu == 0.0)
-		nu = 1.0e-14;
-	if (n_sub == 0.0)
-		n_sub = 1.0;
-
-	bin_size = 1.0/n_sub;  //size of each sub-layer
-	// rescale
-	ind *= bin_size;
-	// modify func so that func(0) =0 and func(max)=1
-	func = 1.0-exp(-nu*ind);
-	// normalize by its max
-	func /= (1.0-exp(-nu));
-
-	return func;
-}
-
-// use the left hand side mirror image of exp func
-double exp_l(double n_sub, double ind, double nu)
-{
-  double bin_size, func;
-	if (nu == 0.0)
-		nu = 1.0e-14;
-	if (n_sub == 0.0)
-		n_sub = 1.0;
-
-	bin_size = 1.0/n_sub;  //size of each sub-layer
-	// rescale
-	ind *= bin_size;
-	// modify func
-	func = exp(-nu*(1.0-ind))-exp(-nu);
-	// normalize by its max
-	func /= (1.0-exp(-nu));
-
-	return func;
-}
-
-// To select function called
-// At nu = 0 (singular point), call line function
-double intersldfunc(int fun_type, double n_sub, double i, double nu, double sld_l, double sld_r)
-{
-	double sld_i, func;
-	// this condition protects an error from the singular point
-	if (nu == 0.0){
-		nu = 1.0e-13;
-	}
-	// select func
-	switch(fun_type){
-		case 1 :
-			func = power_r(n_sub, i, nu);
-			break;
-		case 2 :
-			func = power_l(n_sub, i, nu);
-			break;
-		case 3 :
-			func = exp_r(n_sub, i, nu);
-			break;
-		case 4 :
-			func = exp_l(n_sub, i, nu);
-			break;
-		case 5 :
-			func = linearfunc(n_sub, i, nu);
-			break;
-		default:
-			func = err_mod_func(n_sub, i, nu);
-			break;
-	}
-	// compute sld
-	if (sld_r>sld_l){
-		sld_i = (sld_r-sld_l)*func+sld_l; //sld_cal(sld[i],sld[i+1],n_sub,dz,thick);
-	}
-	else if (sld_r<sld_l){
-		func = 1.0-func;
-		sld_i = (sld_l-sld_r)*func+sld_r; //sld_cal(sld[i],sld[i+1],n_sub,dz,thick);
-	}
-	else{
-		sld_i = sld_r;
-	}
-	return sld_i;
-}
-
-
-// used by refl.c
-double interfunc(int fun_type, double n_sub, double i, double sld_l, double sld_r)
-{
-	double sld_i, func;
-	switch(fun_type){
-		case 0 :
-			func = err_mod_func(n_sub, i, 2.5);
-			break;
-		default:
-			func = linearfunc(n_sub, i, 1.0);
-			break;
-	}
-	if (sld_r>sld_l){
-		sld_i = (sld_r-sld_l)*func+sld_l; //sld_cal(sld[i],sld[i+1],n_sub,dz,thick);
-	}
-	else if (sld_r<sld_l){
-		func = 1.0-func;
-		sld_i = (sld_l-sld_r)*func+sld_r; //sld_cal(sld[i],sld[i+1],n_sub,dz,thick);
-	}
-	else{
-		sld_i = sld_r;
-	}
-	return sld_i;
-}
Index: sasmodels/models/spherical_sld.py
===================================================================
--- sasmodels/models/spherical_sld.py	(revision 54bcd4abeaba4faf03e24319c050286fc10969da)
+++ sasmodels/models/spherical_sld.py	(revision 745b7bb20cb177250d7a71a24058de0dcca609e4)
@@ -1,25 +1,18 @@
 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).
-
-.. figure:: img/spherical_sld_profile.gif
-
-    Exemplary SLD profile
+Similarly to the onion, 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 the error function, power-law, or exponential
+functions.  The scattering intensity is computed by building a continuous
+custom SLD profile along the radius of the particle. The SLD profile is
+composed of a number of uniform shells with interfacial shells between them.
+
+.. figure:: img/spherical_sld_profile.png
+
+    Example 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.
+shells here are sub-divided and numerically integrated assuming each
+sub-shell is described by a line function, with *n_steps* sub-shells per
+interface. The form factor is normalized by the total volume of the sphere.
 
 Definition
@@ -42,5 +35,5 @@
 
 
-so that individual terms can be calcualted as follows:
+so that individual terms can be calculated as follows:
 
 .. math::
@@ -72,5 +65,5 @@
 
 
-Here we assumed that the SLDs of the core and solvent are constant against $r$.
+Here we assumed that the SLDs of the core and solvent are constant in $r$.
 The SLD at the interface between shells, $\rho_{\text {inter}_i}$
 is calculated with a function chosen by an user, where the functions are
@@ -107,7 +100,7 @@
 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
+as well as each sub-shell. Thus B and C are determined.
+
+Once $\rho_{\text{inter}_i}$ is found at the boundary of the sub-shell of the
 interface, we can find its contribution to the form factor $P(q)$
 
@@ -121,5 +114,5 @@
     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 } )
+    ( r_{j} ) V ( r_{ \text{subshell}_j } )
     \Big[ \frac {r_j^2 \beta_\text{out}^2 \sin(\beta_\text{out})
     - (\beta_\text{out}^2-2) \cos(\beta_\text{out}) }
@@ -127,5 +120,5 @@
 
     - 3 ( \rho_{ \text{inter}_i } ( r_{j+1} ) - \rho_{ \text{inter}_i }
-    ( r_{j} ) V ( r_{ \text{sublayer}_j-1 } )
+    ( r_{j} ) V ( r_{ \text{subshell}_j-1 } )
     \Big[ \frac {r_{j-1}^2 \sin(\beta_\text{in})
     - (\beta_\text{in}^2-2) \cos(\beta_\text{in}) }
@@ -152,8 +145,8 @@
 
 
-We assume the $\rho_{\text{inter}_i} (r)$ can be approximately linear
-within a sub-layer $j$
-
-Finally form factor can be calculated by
+We assume $\rho_{\text{inter}_j} (r)$ is approximately linear
+within the sub-shell $j$.
+
+Finally the form factor can be calculated by
 
 .. math::
@@ -169,14 +162,6 @@
     q = \sqrt{q_x^2 + q_y^2}
 
-
-.. figure:: img/spherical_sld_1d.jpg
-
-    1D plot using the default values (w/400 data point).
-
-.. figure:: img/spherical_sld_default_profile.jpg
-
-    SLD profile from the default values.
-
 .. note::
+
     The outer most radius is used as the effective radius for S(Q)
     when $P(Q) * S(Q)$ is applied.
@@ -206,5 +191,5 @@
 # pylint: disable=bad-whitespace, line-too-long
 #            ["name", "units", default, [lower, upper], "type", "description"],
-parameters = [["n_shells",             "",           1,      [1, 11],        "volume", "number of shells"],
+parameters = [["n_shells",             "",           1,      [1, 10],        "volume", "number of shells"],
               ["sld_solvent",          "1e-6/Ang^2", 1.0,    [-inf, inf],    "sld", "solvent sld"],
               ["sld[n_shells]",        "1e-6/Ang^2", 4.06,   [-inf, inf],    "sld", "sld of the shell"],
@@ -258,5 +243,5 @@
     rho.append(sld_solvent)
     # return sld profile (r, beta)
-    return np.asarray(z), np.asarray(rho)*1e-6
+    return np.asarray(z), np.asarray(rho)
 
 
Index: sasmodels/sasview_model.py
===================================================================
--- sasmodels/sasview_model.py	(revision 4edec6f1bf80f86d3d885e9135d676090bf2a652)
+++ sasmodels/sasview_model.py	(revision 745b7bb20cb177250d7a71a24058de0dcca609e4)
@@ -323,14 +323,16 @@
                 beta is a list of the corresponding SLD values
         """
-        args = [] # type: List[Union[float, np.ndarray]]
+        args = {} # type: Dict[str, Any]
         for p in self._model_info.parameters.kernel_parameters:
             if p.id == self.multiplicity_info.control:
-                args.append(float(self.multiplicity))
+                value = float(self.multiplicity)
             elif p.length == 1:
-                args.append(self.params.get(p.id, np.NaN))
+                value = self.params.get(p.id, np.NaN)
             else:
-                args.append([self.params.get(p.id+str(k), np.NaN)
-                             for k in range(1,p.length+1)])
-        return self._model_info.profile(*args)
+                value = np.array([self.params.get(p.id+str(k), np.NaN)
+                                  for k in range(1,p.length+1)]
+            args[p.id] = value
+
+        return self._model_info.profile(**args)
 
     def setParam(self, name, value):
