Index: sasmodels/models/lib/librefl.c
===================================================================
--- sasmodels/models/lib/librefl.c	(revision eb97b11aadaf88044ec058d0176aacaaca44500d)
+++ sasmodels/models/lib/librefl.c	(revision eb97b11aadaf88044ec058d0176aacaaca44500d)
@@ -0,0 +1,202 @@
+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 = erf(ind)/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.c
===================================================================
--- sasmodels/models/spherical_sld.c	(revision 299dcce4251e2668af3f4a5bab4afdc4c69576a7)
+++ sasmodels/models/spherical_sld.c	(revision eb97b11aadaf88044ec058d0176aacaaca44500d)
@@ -5,5 +5,5 @@
     int n_shells);
 
-double Iq(double q,
+static double Iq(double q,
     int n_shells,
     double thick_inter[],
@@ -17,5 +17,5 @@
     double core_radius);
 
-double Iqxy(double qx, double qy,
+static double Iqxy(double qx, double qy,
     int n_shells,
     double thick_inter[],
@@ -56,5 +56,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 ODD number in double //no other number works now
+  double nsl=npts;//21.0; //nsl = Num_sub_layer:  must be ODD double number
   int n_s;
 
@@ -63,19 +63,6 @@
   double pi;
 
-  //int* fun_type;
-  //double* sld;
-  //double* thick_inter;
-  //double* thick;
-  //double* fun_coef;
-
   double total_thick=0.0;
 
-  //fun_type = (int*)malloc((n+2)*sizeof(int));
-  //sld = (double*)malloc((n+2)*sizeof(double));
-  //thick_inter = (double*)malloc((n+2)*sizeof(double));
-  //thick = (double*)malloc((n+2)*sizeof(double));
-  //fun_coef = (double*)malloc((n+2)*sizeof(double));
-
-  //TODO: Solution to avoid mallocs but probablyu can be done better
   int fun_type[12];
   double sld[12];
@@ -175,5 +162,6 @@
             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);
+              fun = sign * r * sph_j1c(qr) + sign * 3.0 * sin(qr)/(qr * qr * q )
+                + sign * 6.0 * cos(qr)/(qr * qr * qr * q);
             }
           }
@@ -285,19 +273,7 @@
 
     double q = sqrt(qx*qx + qy*qy);
-    return Iq(q, n_shells, thick_inter_0, func_inter_0, core0_sld, solvent_sld,
-    flat1_sld, flat2_sld, flat3_sld, flat4_sld, flat5_sld,
-    flat6_sld, flat7_sld, flat8_sld, flat9_sld, flat10_sld,
-    thick_inter_1, thick_inter_2, thick_inter_3, thick_inter_4, thick_inter_5,
-    thick_inter_6, thick_inter_7, thick_inter_8, thick_inter_9, thick_inter_10,
-    thick_flat_1, thick_flat_2, thick_flat_3, thick_flat_4, thick_flat_5,
-    thick_flat_6, thick_flat_7, thick_flat_8, thick_flat_9, thick_flat_10,
-    func_inter_1, func_inter_2, func_inter_3, func_inter_4, func_inter_5,
-    func_inter_6, func_inter_7, func_inter_8, func_inter_9, func_inter_10,
-    nu_inter_1, nu_inter_2, nu_inter_3, nu_inter_4, nu_inter_5,
-    nu_inter_6, nu_inter_7, nu_inter_8, nu_inter_9, nu_inter_10,
-    npts_inter, nu_inter_0, rad_core_0);
-
-    //TODO: Check if evalute rphi is not needed?
-
-}
-
+    return Iq(q, n_shells, thick_inter[], func_inter[], sld_core, sld_solvent,
+    sld_flat[], thick_flat[],nu_inter[], npts_inter, core_radius)
+
+}
+
Index: sasmodels/models/spherical_sld.py
===================================================================
--- sasmodels/models/spherical_sld.py	(revision d2bb6049fb1ca2c1d60bcd3853b1d2ae889866b2)
+++ sasmodels/models/spherical_sld.py	(revision eb97b11aadaf88044ec058d0176aacaaca44500d)
@@ -170,23 +170,97 @@
 # pylint: disable=bad-whitespace, line-too-long
 #            ["name", "units", default, [lower, upper], "type", "description"],
-parameters = [["n",                "",               1,      [0, 9],         "", "number of shells"],
+parameters = [["n_shells",                "",               1,      [0, 9],         "", "number of shells"],
               ["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"],
               ["sld_flat[n]",      "1e-6/Ang^2",     4.06,   [-inf, inf],    "", "sld function flat"],
               ["thick_flat[n]",    "Ang",            100.0,  [0, inf],       "", "flat layer_thickness"],
               ["func_inter[n]",    "",               0,      [0, 4],         "", "Erf:0, RPower:1, LPower:2, RExp:3, LExp:4"],
               ["thick_inter[n]",   "Ang",            50.0,   [0, inf],       "", "intern layer thickness"],
-              ["inter_nu[n]",      "",               2.5,    [-inf, inf],    "", "steepness parameter"],
+              ["nu_inter[n]",      "",               2.5,    [-inf, inf],    "", "steepness parameter"],
               ["npts_inter",       "",               35,     [0, 35],        "", "number of points in each sublayer Must be odd number"],
-              ["sld_solvent",      "1e-6/Ang^2",     1.0,    [-inf, inf],    "", "sld function solvent"],
               ]
 # pylint: enable=bad-whitespace, line-too-long
-#source = ["lib/librefl.c",  "lib/sph_j1c.c", "spherical_sld.c"]
-
-def Iq(q, *args, **kw):
-    return q
-
-def Iqxy(qx, *args, **kw):
-    return qx
+source = ["lib/librefl.c",  "lib/sph_j1c.c", "spherical_sld.c"]
+
+#def Iq(q, *args, **kw):
+#    return q
+
+#def Iqxy(qx, *args, **kw):
+#    return qx
+
+def profile(n_shells, radius_core,  sld_core,  sld_solvent, sld_flat,
+            thick_flat, func_inter, thick_inter, nu_inter, npts_inter):
+    """
+    Returns shape profile with x=radius, y=SLD.
+    """
+
+    #total_radius = 1.25*(sum(thickness[:n_shells]) + core_radius + 1)
+    #dr = total_radius/400  # 400 points for a smooth plot
+
+    z = []
+    beta = []
+    z0 = 0
+    # two sld points for core
+    z.append(0)
+    beta.append(sld_core)
+    z.append(radius_core)
+    beta.append(sld_core)
+    z0 += radius_core
+
+    for i in range(1, n_shells+2):
+        dz = thick_inter[i-1]/npts_inter
+        # j=0 for interface, j=1 for flat layer
+        for j in range(0, 2):
+            # interation for sub-layers
+            for n_s in range(0, npts_inter+1):
+                if j == 1:
+                    if i == n_shells+1:
+                        break
+                    # shift half sub thickness for the first point
+                    z0 -= dz#/2.0
+                        z.append(z0)
+                        #z0 -= dz/2.0
+                        z0 += thick_flat[i]
+                        sld_i = sld_flat[i]
+                        beta.append(sld_flat[i])
+                        dz = 0
+                    else:
+                        nu = nu_inter[i-1]
+                        # decide which sld is which, sld_r or sld_l
+                        if i == 1:
+                            sld_l = sld_core
+                        else:
+                            sld_l = sld_flat[i-1]
+                        if i == n_shells+1:
+                            sld_r = sld_solvent
+                        else:
+                            sld_r = sld_flat[i]
+                        # get function type
+                        func_idx = func_inter[i-1]
+                        # calculate the sld
+                        sld_i = intersldfunc(func_idx, npts_inter, n_s, nu,
+                                            sld_l, sld_r)
+                    # append to the list
+                    z.append(z0)
+                    beta.append(sld_i)
+                    z0 += dz
+                    if j == 1:
+                        break
+    # put sld of solvent
+    z.append(z0)
+    beta.append(sld_solvent)
+    z_ext = z0/5.0
+    z.append(z0+z_ext)
+    beta.append(sld_solvent)
+    # return sld profile (r, beta)
+    
+    return np.asarray(z), np.asarray(beta)*1e-6
+
+def ER(core_radius, n, thickness):
+    return np.sum(thickness[:n[0]], axis=0) + core_radius
+
+def VR(core_radius, n, thickness):
+    return 1.0, 1.0
 
 
