Index: sasmodels/models/lib/librefl.c
===================================================================
--- sasmodels/models/lib/librefl.c	(revision 9244430cab1a706230c8d22de65e4c66b14526cd)
+++ sasmodels/models/lib/librefl.c	(revision 9244430cab1a706230c8d22de65e4c66b14526cd)
@@ -0,0 +1,335 @@
+//complex cassign(real, imag);
+//complex cplx_add(x,y);
+//complex rcmult(x,y);
+//complex cplx_sub(x,y);
+//complex cplx_mult(x,y);
+//complex cplx_div(x,y);
+//complex cplx_exp(b);
+
+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 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);
+
+complex cassign(real, imag)
+	double real, imag;
+{
+	complex x;
+	x.re = real;
+	x.im = imag;
+	return x;
+}
+
+
+complex cplx_add(x,y)
+	complex x,y;
+{
+	complex z;
+	z.re = x.re + y.re;
+	z.im = x.im + y.im;
+	return z;
+}
+
+complex rcmult(x,y)
+	double x;
+    complex y;
+{
+	complex z;
+	z.re = x*y.re;
+	z.im = x*y.im;
+	return z;
+}
+
+complex cplx_sub(x,y)
+	complex x,y;
+{
+	complex z;
+	z.re = x.re - y.re;
+	z.im = x.im - y.im;
+	return z;
+}
+
+
+complex cplx_mult(x,y)
+	complex x,y;
+{
+	complex z;
+	z.re = x.re*y.re - x.im*y.im;
+	z.im = x.re*y.im + x.im*y.re;
+	return z;
+}
+
+complex cplx_div(x,y)
+	complex x,y;
+{
+	complex z;
+	z.re = (x.re*y.re + x.im*y.im)/(y.re*y.re + y.im*y.im);
+	z.im = (x.im*y.re - x.re*y.im)/(y.re*y.re + y.im*y.im);
+	return z;
+}
+
+complex cplx_exp(b)
+	complex b;
+{
+	complex z;
+	double br,bi;
+	br=b.re;
+	bi=b.im;
+	z.re = exp(br)*cos(bi);
+	z.im = exp(br)*sin(bi);
+	return z;
+}
+
+
+complex cplx_sqrt(z)    //see Schaum`s Math Handbook p. 22, 6.6 and 6.10
+	complex z;
+{
+	complex c;
+	double zr,zi,x,y,r,w;
+
+	zr=z.re;
+	zi=z.im;
+
+	if (zr==0.0 && zi==0.0)
+	{
+    c.re=0.0;
+	c.im=0.0;
+    return c;
+	}
+	else
+	{
+		x=fabs(zr);
+		y=fabs(zi);
+		if (x>y)
+		{
+			r=y/x;
+			w=sqrt(x)*sqrt(0.5*(1.0+sqrt(1.0+r*r)));
+		}
+		else
+		{
+			r=x/y;
+			w=sqrt(y)*sqrt(0.5*(r+sqrt(1.0+r*r)));
+		}
+		if (zr >=0.0)
+		{
+			c.re=w;
+			c.im=zi/(2.0*w);
+		}
+		else
+		{
+			c.im=(zi >= 0) ? w : -w;
+			c.re=zi/(2.0*c.im);
+		}
+		return c;
+	}
+}
+
+complex cplx_cos(b)
+	complex b;
+{
+	complex zero,two,z,i,bi,negbi;
+	zero = cassign(0.0,0.0);
+	two = cassign(2.0,0.0);
+	i = cassign(0.0,1.0);
+	bi = cplx_mult(b,i);
+	negbi = cplx_sub(zero,bi);
+	z = cplx_div(cplx_add(cplx_exp(bi),cplx_exp(negbi)),two);
+	return z;
+}
+
+// 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 = 1e-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 = 1e-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 = 1e-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 = 1e-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 = 1e-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 = 1e-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 9244430cab1a706230c8d22de65e4c66b14526cd)
+++ sasmodels/models/spherical_sld.py	(revision 9244430cab1a706230c8d22de65e4c66b14526cd)
@@ -0,0 +1,205 @@
+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).
+Unlike the OnionExpShellModel (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
+in the GUI. The form factor is normalized by the total volume of the sphere.
+
+Definition
+----------
+
+The scattering intensity $I(q)$ in 1D is calculated as:
+
+.. math::
+
+    P(q) = \frac{f^2}{V_\text{particle}}
+    f = f_\text{core} + \sum_{\text{inter}_i=0}^N f_text{inter}_i +
+    \sum_{\text{flat}_i=0}^N f_text{flat}_i +f_\text{solvent}
+
+where, for a spherically symmetric particle with a particle density \rho(r)
+
+
+The scaling of the second power law region (coefficent C) is then automatically scaled to
+match the first by following formula
+
+.. math::
+    C = \frac{A}{qc^{-m1} qc^{-m2}}
+
+.. note::
+    Be sure to enter the power law exponents as positive values!
+
+For 2D data the scattering intensity is calculated in the same way as 1D,
+where the $q$ vector is defined as
+
+.. math::
+
+    q = \sqrt{q_x^2 + q_y^2}
+
+
+.. figure:: img/two_power_law_1d.jpg
+
+    1D plot using the default values (w/500 data point).
+
+References
+----------
+L A Feigin and D I Svergun, Structure Analysis by Small-Angle X-Ray and Neutron Scattering, Plenum Press, New York, (1987)
+
+"""
+
+from numpy import inf
+
+name = "spherical_sld"
+title = "Sperical SLD intensity calculation"
+description = """
+            I(q) =
+               background = Incoherent background [1/cm]
+        """
+category = "shere-based"
+
+# pylint: disable=bad-whitespace, line-too-long
+#            ["name", "units", default, [lower, upper], "type", "description"],
+parameters = [["n_shells",       "",         1, [0, 9], "", "number of shells"],
+              ["npts_imter",  "",         35, [0, inf], "", "number of points in each sublayer"],
+              ["sld_solv",     "1/Ang^2",    1E-6,[-inf, inf], "","sld function solvent"],
+              ["func_inter_0",       "",         0, [0, 4], "", "'Erf(|nu|*z)':0, 'RPower(z^|nu|)':1, 'LPower(z^|nu|)':2, 'RExp(-|nu|*z)':3, 'LExp(-|nu|*z)':4"],
+              ["nu_inter_0",       "",         2.5, [-inf, inf], "", "steepness parameter"],
+              ["thick_inter_0",       "Ang",         50, [-inf, inf], "", "intern layer thickness"],
+              ["sld_core_0",       "1/Ang^2",         2.07E-6, [-inf, inf],"", "sld function flat"],
+              ["thick_flat_0",       "Ang",         100, [-inf, inf], "", "flat layer_thickness"],
+              ["func_inter_1",       "",         0, [0, 4], "", "'Erf(|nu|*z)':0, 'RPower(z^|nu|)':1, 'LPower(z^|nu|)':2, 'RExp(-|nu|*z)':3, 'LExp(-|nu|*z)':4"],
+              ["nu_inter_1",       "",         2.5, [-inf, inf], "", "steepness parameter"],
+              ["thick_inter_1",       "Ang",         50.0, [-inf, inf], "", "intern layer thickness"],
+              ["sld_flat_1",       "1/Ang^2",         4.06E-6, [-inf, inf],"", "sld function flat"],
+              ["thick_flat_1",       "Ang",         100.0, [-inf, inf], "", "flat layer_thickness"],
+              ["func_inter_2",       "",         0, [0, 4], "", "'Erf(|nu|*z)':0, 'RPower(z^|nu|)':1, 'LPower(z^|nu|)':2, 'RExp(-|nu|*z)':3, 'LExp(-|nu|*z)':4"],
+              ["nu_inter_2",       "",         2.5, [-inf, inf], "", "steepness parameter"],
+              ["thick_inter_2",       "Ang",         50.0, [-inf, inf], "", "intern layer thickness"],
+              ["sld_flat_2",       "1/Ang^2",         3.56E-6, [-inf, inf],"", "sld function flat"],
+              ["thick_flat_2",       "Ang",         100.0, [-inf, inf], "", "flat layer_thickness"],
+              ["func_inter_3",       "",         0, [0, 4], "", "'Erf(|nu|*z)':0, 'RPower(z^|nu|)':1, 'LPower(z^|nu|)':2, 'RExp(-|nu|*z)':3, 'LExp(-|nu|*z)':4"],
+              ["nu_inter_3",       "",         2.5, [-inf, inf], "", "steepness parameter"],
+              ["thick_inter_3",       "Ang",         50.0, [-inf, inf], "", "intern layer thickness"],
+              ["sld_flat_3",       "1/Ang^2",         4.06E-6, [-inf, inf],"", "sld function flat"],
+              ["thick_flat_3",       "Ang",         100.0, [-inf, inf], "", "flat layer_thickness"],
+              ["func_inter_4",       "",         0, [0, 4], "", "'Erf(|nu|*z)':0, 'RPower(z^|nu|)':1, 'LPower(z^|nu|)':2, 'RExp(-|nu|*z)':3, 'LExp(-|nu|*z)':4"],
+              ["nu_inter_4",       "",         2.5, [-inf, inf], "", "steepness parameter"],
+              ["thick_inter_4",       "Ang",         50.0, [-inf, inf], "", "intern layer thickness"],
+              ["sld_flat_4",       "1/Ang^2",         3.5E-6, [-inf, inf],"", "sld function flat"],
+              ["thick_flat_4",       "Ang",         100.0, [-inf, inf], "", "flat layer_thickness"],
+              ["func_inter_5",       "",         0, [0, 4], "", "'Erf(|nu|*z)':0, 'RPower(z^|nu|)':1, 'LPower(z^|nu|)':2, 'RExp(-|nu|*z)':3, 'LExp(-|nu|*z)':4"],
+              ["nu_inter_5",       "",         2.5, [-inf, inf], "", "steepness parameter"],
+              ["thick_inter_5",       "Ang",         50.0, [-inf, inf], "", "intern layer thickness"],
+              ["sld_flat_5",       "1/Ang^2",         4.06E-6, [-inf, inf],"", "sld function flat"],
+              ["thick_flat_5",       "Ang",         100.0, [-inf, inf], "", "flat layer_thickness"],
+              ["func_inter_6",       "",         0, [0, 4], "", "'Erf(|nu|*z)':0, 'RPower(z^|nu|)':1, 'LPower(z^|nu|)':2, 'RExp(-|nu|*z)':3, 'LExp(-|nu|*z)':4"],
+              ["nu_inter_6",       "",         2.5, [-inf, inf], "", "steepness parameter"],
+              ["thick_inter_6",       "Ang",         50.0, [-inf, inf], "", "intern layer thickness"],
+              ["sld_flat_6",       "1/Ang^2",         3.5E-6, [-inf, inf],"", "sld function flat"],
+              ["thick_flat_6",       "Ang",         100.0, [-inf, inf], "", "flat layer_thickness"],
+              ["func_inter_7",       "",         0, [0, 4], "", "'Erf(|nu|*z)':0, 'RPower(z^|nu|)':1, 'LPower(z^|nu|)':2, 'RExp(-|nu|*z)':3, 'LExp(-|nu|*z)':4"],
+              ["nu_inter_7",       "",         2.5, [-inf, inf], "", "steepness parameter"],
+              ["thick_inter_7",       "Ang",         50.0, [-inf, inf], "", "intern layer thickness"],
+              ["sld_flat_7",       "1/Ang^2",         4.06E-6, [-inf, inf],"", "sld function flat"],
+              ["thick_flat_7",       "Ang",         100.0, [-inf, inf], "", "flat layer_thickness"],
+              ["func_inter_8",       "",         0, [0, 4], "", "'Erf(|nu|*z)':0, 'RPower(z^|nu|)':1, 'LPower(z^|nu|)':2, 'RExp(-|nu|*z)':3, 'LExp(-|nu|*z)':4"],
+              ["nu_inter_8",       "",         2.5, [-inf, inf], "", "steepness parameter"],
+              ["thick_inter_8",       "Ang",         50.0, [-inf, inf], "", "intern layer thickness"],
+              ["sld_flat_8",       "1/Ang^2",         3.5E-6, [-inf, inf],"", "sld function flat"],
+              ["thick_flat_8",       "Ang",         100.0, [-inf, inf], "", "flat layer_thickness"],
+              ["func_inter_9",       "",         0, [0, 4], "", "'Erf(|nu|*z)':0, 'RPower(z^|nu|)':1, 'LPower(z^|nu|)':2, 'RExp(-|nu|*z)':3, 'LExp(-|nu|*z)':4"],
+              ["nu_inter_9",       "",         2.5, [-inf, inf], "", "steepness parameter"],
+              ["thick_inter_9",       "Ang",         50.0, [-inf, inf], "", "intern layer thickness"],
+              ["sld_flat_9",       "1/Ang^2",         4.06E-6, [-inf, inf],"", "sld function flat"],
+              ["thick_flat_9",       "Ang",         100.0, [-inf, inf], "", "flat layer_thickness"],
+              ]
+# pylint: enable=bad-whitespace, line-too-long
+source = ["lib/librefl.c", "spherical_sld.c"]
+
+
+#TODO: Not clear if dispersion function is needed
+#def _set_dispersion(self):
+#        """
+#        Setting dispersion for all shells
+#        """
+#        ##set dispersion from model
+#        for name , value in self.model.dispersion.iteritems():
+#
+#            nshell = -1
+#            if name.split('_')[0] == 'thick':
+#                while nshell < 1:
+#                    nshell += 1
+#                    if name.split('_')[1] == 'inter%s' % str(nshell):
+#                        self.dispersion[name] = value
+#                    else:
+#                        continue
+#            else:
+#                self.dispersion[name] = value
+
+
+def ER(radius):
+        """
+        Calculate the effective radius for P(q)*S(q)
+        Tale different radius values from corresponding shells
+
+        :return: the value of the effective radius
+        """
+
+        return radius
+
+
+demo = dict(scale=1, background=0.1,
+        n_shells=1,
+        npts_inter=35,
+        sld_solv=1E-6,
+        func_inter_1=0.0,
+        nu_inter_1=2.5,
+        thick_inter_1=50,
+        sld_flat_1=4E-6,
+        thick_flat_1=100,
+        func_inter_0=0.0,
+        nu_inter_0=2.5,
+        rad_core_0=50.0,
+        std_core_0=2.07E-6,
+        thick_core_0=50
+    )
+
+oldname = "SphericalSLDModel"
+oldpars = dict(
+        npts_inter='npts_imter',
+        sld_solv='sld_solv',
+        func_inter_1='func_inter_1',
+        nu_inter_1='nu_inter',
+        thick_inter_1='thick_inter1',
+        sld_flat_1='sld_flat1',
+        thick_flat_1='thick_flat1',
+        func_inter_0='nu_inter0',
+        nu_inter_0='nu_inter0',
+        rad_core_0='rad_core0',
+        std_core_0='std_core0',
+        thick_core_0='thick_core0',
+        background='background')
+
+tests = [
+    # Accuracy tests based on content in test/utest_extra_models.py
+    [{'npts_iter':35,
+        'sld_solv':1E-6,
+        'func_inter_1':0.0,
+        'nu_inter':2.5,
+        'thick_inter_1':50,
+        'sld_flat_1':4E-6,
+        'thick_flat_1':100,
+        'func_inter_0':0.0,
+        'nu_inter_0':2.5,
+        'rad_core_0':50.0,
+        'std_core_0':2.07E-6,
+        'thick_core':50,
+        'background': 0.0,
+    }, 0.001, 1000],
+]
