Index: sasmodels/models/_onion.py
===================================================================
--- sasmodels/models/_onion.py	(revision f22e50e7d1bc8ecbd1dc21b3650d30b579aca290)
+++ sasmodels/models/_onion.py	(revision f22e50e7d1bc8ecbd1dc21b3650d30b579aca290)
@@ -293,37 +293,35 @@
 
 #             ["name", "units", default, [lower, upper], "type","description"],
-parameters = [["core_sld", "1e-6/Ang^2", 1.0, [-inf, inf], "",
+parameters = [["sld_core", "1e-6/Ang^2", 1.0, [-inf, inf], "",
                "Core scattering length density"],
               ["core_radius", "Ang", 200., [0, inf], "volume",
                "Radius of the core"],
-              ["solvent_sld", "1e-6/Ang^2", 6.4, [-inf, inf], "",
+              ["sld_solvent", "1e-6/Ang^2", 6.4, [-inf, inf], "",
                "Solvent scattering length density"],
-              ["n", "", 1, [0, 10], "volume",
+              ["n_shells", "", 1, [0, 10], "volume",
                "number of shells"],
-              ["in_sld[n]", "1e-6/Ang^2", 1.7, [-inf, inf], "",
+              ["sld_in[n_shells]", "1e-6/Ang^2", 1.7, [-inf, inf], "",
                "scattering length density at the inner radius of shell k"],
-              ["out_sld[n]", "1e-6/Ang^2", 2.0, [-inf, inf], "",
+              ["sld_out[n_shells]", "1e-6/Ang^2", 2.0, [-inf, inf], "",
                "scattering length density at the outer radius of shell k"],
-              ["thickness[n]", "Ang", 40., [0, inf], "volume",
+              ["thickness[n_shells]", "Ang", 40., [0, inf], "volume",
                "Thickness of shell k"],
-              ["A[n]", "", 1.0, [-inf, inf], "",
+              ["A[n_shells]", "", 1.0, [-inf, inf], "",
                "Decay rate 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):
+source = ["lib/sph_j1c.c", "onion.c"]
+
+#def Iq(q, *args, **kw):
+#    return q
+
+profile_axes = ['Radius (A)', 'SLD (1e-6/A^2)']
+def profile(core_sld, core_radius, solvent_sld, n_shells,
+            in_sld, out_sld, thickness, A):
     """
-    SLD profile
+    Returns shape profile with x=radius, y=SLD.
     """
 
-    total_radius = 1.25*(sum(thickness[:n]) + core_radius + 1)
+    total_radius = 1.25*(sum(thickness[:n_shells]) + core_radius + 1)
     dr = total_radius/400  # 400 points for a smooth plot
 
@@ -338,5 +336,5 @@
 
     # add in the shells
-    for k in range(n):
+    for k in range(n_shells):
         # Left side of each shells
         r0 = r[-1]
@@ -365,5 +363,5 @@
     beta.append(solvent_sld)
 
-    return np.asarray(r), np.asarray(beta)
+    return np.asarray(r), np.asarray(beta)*1e-6
 
 def ER(core_radius, n, thickness):
@@ -374,10 +372,10 @@
 
 demo = {
-    "solvent_sld": 2.2,
-    "core_sld": 1.0,
+    "sld_solvent": 2.2,
+    "sld_core": 1.0,
     "core_radius": 100,
     "n": 4,
-    "in_sld": [0.5, 1.5, 0.9, 2.0],
-    "out_sld": [nan, 0.9, 1.2, 1.6],
+    "sld_in": [0.5, 1.5, 0.9, 2.0],
+    "sld_out": [nan, 0.9, 1.2, 1.6],
     "thickness": [50, 75, 150, 75],
     "A": [0, -1, 1e-4, 1],
Index: sasmodels/models/_spherical_sld.py
===================================================================
--- sasmodels/models/_spherical_sld.py	(revision f22e50e7d1bc8ecbd1dc21b3650d30b579aca290)
+++ sasmodels/models/_spherical_sld.py	(revision f22e50e7d1bc8ecbd1dc21b3650d30b579aca290)
@@ -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,12 +179,15 @@
 
 .. 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)
 
 """
 
+import numpy as np
 from numpy import inf
 
@@ -170,49 +202,121 @@
 # pylint: disable=bad-whitespace, line-too-long
 #            ["name", "units", default, [lower, upper], "type", "description"],
-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_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"],
-              ["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"],
+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],       "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],       "volume", "intern layer thickness"],
+              ["nu_inter[n_shells]",      "",               2.5,    [-inf, inf],    "", "steepness parameter"],
               ]
 # 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
-
-
-demo = dict(
-    n_shells=4,
-    scale=1.0,
-    solvent_sld=1.0,
-    background=0.0,
-    npts_inter=35.0,
-    )
+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):
+    """
+    Returns shape profile with x=radius, y=SLD.
+    """
+
+    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
+    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(n_shells, radius_core, thick_inter0, thick_inter, thick_flat):
+    total_thickness = thick_inter0
+    total_thickness += np.sum(thick_inter[:n_shells], axis=0)
+    total_thickness += np.sum(thick_flat[:n_shells], axis=0)
+    return total_thickness + radius_core
+
+
+demo = {
+    "n_shells": 4,
+    "npts_inter": 35.0,
+    "radius_core": 50.0,
+    "sld_core": 2.07,
+    "sld_solvent": 1.0,
+    "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],
+    }
 
 #TODO: Not working yet
 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],
 ]
Index: sasmodels/models/cylinder.c
===================================================================
--- sasmodels/models/cylinder.c	(revision 26141cba42ea53a3ca3b54e60b579f5525dd583d)
+++ sasmodels/models/cylinder.c	(revision 03cac087353e872f40154bf57ef8161eb5367fd4)
@@ -3,4 +3,6 @@
 double Iqxy(double qx, double qy, double sld, double solvent_sld,
     double radius, double length, double theta, double phi);
+
+#define INVALID(v) (v.radius<0 || v.length<0)
 
 double form_volume(double radius, double length)
@@ -15,5 +17,4 @@
     double length)
 {
-    // TODO: return NaN if radius<0 or length<0?
     // precompute qr and qh to save time in the loop
     const double qr = q*radius;
@@ -47,5 +48,4 @@
     double phi)
 {
-    // TODO: return NaN if radius<0 or length<0?
     double sn, cn; // slots to hold sincos function output
 
Index: sasmodels/models/cylinder.py
===================================================================
--- sasmodels/models/cylinder.py	(revision f2473145a6baa12aa3ab15421a80d2f5d004935c)
+++ sasmodels/models/cylinder.py	(revision 7ae2b7f85895f1eddd6ef18f37ff501467a34961)
@@ -82,6 +82,6 @@
 """
 
-import numpy as np
-from numpy import pi, inf
+import numpy as np  # type: ignore
+from numpy import pi, inf  # type: ignore
 
 name = "cylinder"
Index: sasmodels/models/flexible_cylinder.c
===================================================================
--- sasmodels/models/flexible_cylinder.c	(revision e6408d049d3598925e345d8adec1cdb0bb0c6adb)
+++ sasmodels/models/flexible_cylinder.c	(revision e6408d049d3598925e345d8adec1cdb0bb0c6adb)
@@ -1,59 +1,20 @@
-double form_volume(double length, double kuhn_length, double radius);
-double Iq(double q, double length, double kuhn_length, double radius,
-          double sld, double solvent_sld);
-double Iqxy(double qx, double qy, double length, double kuhn_length,
-            double radius, double sld, double solvent_sld);
-double flexible_cylinder_kernel(double q, double length, double kuhn_length,
-                                double radius, double sld, double solvent_sld);
-
-
-double form_volume(double length, double kuhn_length, double radius)
+static double
+form_volume(length, kuhn_length, radius)
 {
     return 1.0;
 }
 
-double flexible_cylinder_kernel(double q,
-          double length,
-          double kuhn_length,
-          double radius,
-          double sld,
-          double solvent_sld)
+static double
+Iq(double q,
+   double length,
+   double kuhn_length,
+   double radius,
+   double sld,
+   double solvent_sld)
 {
-
-    const double cont = sld-solvent_sld;
-    const double qr = q*radius;
-    //const double crossSect = (2.0*J1(qr)/qr)*(2.0*J1(qr)/qr);
-    const double crossSect = sas_J1c(qr);
-    double flex = Sk_WR(q,length,kuhn_length);
-    flex *= crossSect*crossSect;
-    flex *= M_PI*radius*radius*length;
-    flex *= cont*cont;
-    flex *= 1.0e-4;
-    return flex;
+    const double contrast = sld - solvent_sld;
+    const double cross_section = sas_J1c(q*radius);
+    const double volume = M_PI*radius*radius*length;
+    const double flex = Sk_WR(q, length, kuhn_length);
+    return 1.0e-4 * volume * square(contrast*cross_section) * flex;
 }
-
-double Iq(double q,
-          double length,
-          double kuhn_length,
-          double radius,
-          double sld,
-          double solvent_sld)
-{
-
-    double result = flexible_cylinder_kernel(q, length, kuhn_length, radius, sld, solvent_sld);
-    return result;
-}
-
-double Iqxy(double qx, double qy,
-            double length,
-            double kuhn_length,
-            double radius,
-            double sld,
-            double solvent_sld)
-{
-    double q;
-    q = sqrt(qx*qx+qy*qy);
-    double result = flexible_cylinder_kernel(q, length, kuhn_length, radius, sld, solvent_sld);
-
-    return result;
-}
Index: sasmodels/models/gel_fit.c
===================================================================
--- sasmodels/models/gel_fit.c	(revision 30b4ddfb88f4c0a299757def2ee628c2a3d50a42)
+++ sasmodels/models/gel_fit.c	(revision 03cac087353e872f40154bf57ef8161eb5367fd4)
@@ -1,19 +1,3 @@
-double form_volume(void);
-
-double Iq(double q,
-          double guinier_scale,
-          double lorentzian_scale,
-          double gyration_radius,
-          double fractal_exp,
-          double cor_length);
-
-double Iqxy(double qx, double qy,
-          double guinier_scale,
-          double lorentzian_scale,
-          double gyration_radius,
-          double fractal_exp,
-          double cor_length);
-
-static double _gel_fit_kernel(double q,
+static double Iq(double q,
           double guinier_scale,
           double lorentzian_scale,
@@ -24,5 +8,5 @@
     // Lorentzian Term
     ////////////////////////double a(x[i]*x[i]*zeta*zeta);
-    double lorentzian_term = q*q*cor_length*cor_length;
+    double lorentzian_term = square(q*cor_length);
     lorentzian_term = 1.0 + ((fractal_exp + 1.0)/3.0)*lorentzian_term;
     lorentzian_term = pow(lorentzian_term, (fractal_exp/2.0) );
@@ -30,5 +14,5 @@
     // Exponential Term
     ////////////////////////double d(x[i]*x[i]*rg*rg);
-    double exp_term = q*q*gyration_radius*gyration_radius;
+    double exp_term = square(q*gyration_radius);
     exp_term = exp(-1.0*(exp_term/3.0));
 
@@ -37,39 +21,2 @@
     return result;
 }
-double form_volume(void){
-    // Unused, so free to return garbage.
-    return NAN;
-}
-
-double Iq(double q,
-          double guinier_scale,
-          double lorentzian_scale,
-          double gyration_radius,
-          double fractal_exp,
-          double cor_length)
-{
-    return _gel_fit_kernel(q,
-                          guinier_scale,
-                          lorentzian_scale,
-                          gyration_radius,
-                          fractal_exp,
-                          cor_length);
-}
-
-// Iqxy is never called since no orientation or magnetic parameters.
-double Iqxy(double qx, double qy,
-          double guinier_scale,
-          double lorentzian_scale,
-          double gyration_radius,
-          double fractal_exp,
-          double cor_length)
-{
-    double q = sqrt(qx*qx + qy*qy);
-    return _gel_fit_kernel(q,
-                          guinier_scale,
-                          lorentzian_scale,
-                          gyration_radius,
-                          fractal_exp,
-                          cor_length);
-}
-
Index: sasmodels/models/hardsphere.py
===================================================================
--- sasmodels/models/hardsphere.py	(revision ec45c4f5ff6e17014d968da6eb80308c2c8791f0)
+++ sasmodels/models/hardsphere.py	(revision d2bb6049fb1ca2c1d60bcd3853b1d2ae889866b2)
@@ -149,9 +149,4 @@
    """
 
-Iqxy = """
-    // never called since no orientation or magnetic parameters.
-    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);
-    """
-
 # ER defaults to 0.0
 # VR defaults to 1.0
Index: sasmodels/models/hayter_msa.py
===================================================================
--- sasmodels/models/hayter_msa.py	(revision ec45c4f5ff6e17014d968da6eb80308c2c8791f0)
+++ sasmodels/models/hayter_msa.py	(revision d2bb6049fb1ca2c1d60bcd3853b1d2ae889866b2)
@@ -87,8 +87,4 @@
     return 1.0;
     """
-Iqxy = """
-    // never called since no orientation or magnetic parameters.
-    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);
-    """
 # ER defaults to 0.0
 # VR defaults to 1.0
Index: sasmodels/models/lamellar.py
===================================================================
--- sasmodels/models/lamellar.py	(revision 500128bae04dc9573c5f85c656f5ded54a1fa9a1)
+++ sasmodels/models/lamellar.py	(revision 500128bae04dc9573c5f85c656f5ded54a1fa9a1)
@@ -83,8 +83,4 @@
     """
 
-Iqxy = """
-    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);
-    """
-
 # ER defaults to 0.0
 # VR defaults to 1.0
Index: sasmodels/models/lamellar_hg.py
===================================================================
--- sasmodels/models/lamellar_hg.py	(revision ec45c4f5ff6e17014d968da6eb80308c2c8791f0)
+++ sasmodels/models/lamellar_hg.py	(revision d2bb6049fb1ca2c1d60bcd3853b1d2ae889866b2)
@@ -101,8 +101,4 @@
     """
 
-Iqxy = """
-    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);
-    """
-
 # ER defaults to 0.0
 # VR defaults to 1.0
Index: sasmodels/models/lamellar_hg_stack_caille.py
===================================================================
--- sasmodels/models/lamellar_hg_stack_caille.py	(revision ec45c4f5ff6e17014d968da6eb80308c2c8791f0)
+++ sasmodels/models/lamellar_hg_stack_caille.py	(revision d2bb6049fb1ca2c1d60bcd3853b1d2ae889866b2)
@@ -120,8 +120,4 @@
     """
 
-Iqxy = """
-    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);
-    """
-
 # ER defaults to 0.0
 # VR defaults to 1.0
Index: sasmodels/models/lamellar_stack_caille.py
===================================================================
--- sasmodels/models/lamellar_stack_caille.py	(revision ec45c4f5ff6e17014d968da6eb80308c2c8791f0)
+++ sasmodels/models/lamellar_stack_caille.py	(revision d2bb6049fb1ca2c1d60bcd3853b1d2ae889866b2)
@@ -104,8 +104,4 @@
     """
 
-Iqxy = """
-    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);
-    """
-
 # ER defaults to 0.0
 # VR defaults to 1.0
Index: sasmodels/models/lamellar_stack_paracrystal.py
===================================================================
--- sasmodels/models/lamellar_stack_paracrystal.py	(revision ec45c4f5ff6e17014d968da6eb80308c2c8791f0)
+++ sasmodels/models/lamellar_stack_paracrystal.py	(revision d2bb6049fb1ca2c1d60bcd3853b1d2ae889866b2)
@@ -132,8 +132,4 @@
     """
 
-Iqxy = """
-    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);
-    """
-
 # ER defaults to 0.0
 # VR defaults to 1.0
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/lib/sas_JN.c
===================================================================
--- sasmodels/models/lib/sas_JN.c	(revision e6408d049d3598925e345d8adec1cdb0bb0c6adb)
+++ sasmodels/models/lib/sas_JN.c	(revision e6408d049d3598925e345d8adec1cdb0bb0c6adb)
@@ -48,6 +48,7 @@
 */
 
-static double
-sas_JN( int n, double x ) {
+double sas_JN( int n, double x );
+
+double sas_JN( int n, double x ) {
 
     const double MACHEP = 1.11022302462515654042E-16;
Index: sasmodels/models/lib/sph_j1c.c
===================================================================
--- sasmodels/models/lib/sph_j1c.c	(revision e6f14104e82489f43d4acc1eeda0d50ad5408fb6)
+++ sasmodels/models/lib/sph_j1c.c	(revision ba32cddebeff04b66d828ca58b1843244d4425fb)
@@ -7,4 +7,5 @@
 * using double precision that are the source.
 */
+double sph_j1c(double q);
 
 // The choice of the number of terms in the series and the cutoff value for
@@ -43,5 +44,4 @@
 #endif
 
-double sph_j1c(double q);
 double sph_j1c(double q)
 {
Index: sasmodels/models/lib/sphere_form.c
===================================================================
--- sasmodels/models/lib/sphere_form.c	(revision ad90df9f3ed42019113506f1366fe4b851b5f389)
+++ sasmodels/models/lib/sphere_form.c	(revision ba32cddebeff04b66d828ca58b1843244d4425fb)
@@ -1,10 +1,11 @@
-inline double
-sphere_volume(double radius)
+double sphere_volume(double radius);
+double sphere_form(double q, double radius, double sld, double solvent_sld);
+
+double sphere_volume(double radius)
 {
     return M_4PI_3*cube(radius);
 }
 
-inline double
-sphere_form(double q, double radius, double sld, double solvent_sld)
+double sphere_form(double q, double radius, double sld, double solvent_sld)
 {
     const double fq = sphere_volume(radius) * sph_j1c(q*radius);
Index: sasmodels/models/lib/wrc_cyl.c
===================================================================
--- sasmodels/models/lib/wrc_cyl.c	(revision e7678b2f993854334aca4c45e1c7022b2fdbb509)
+++ sasmodels/models/lib/wrc_cyl.c	(revision ba32cddebeff04b66d828ca58b1843244d4425fb)
@@ -2,4 +2,7 @@
     Functions for WRC implementation of flexible cylinders
 */
+double Sk_WR(double q, double L, double b);
+
+
 static double
 AlphaSquare(double x)
@@ -363,5 +366,4 @@
 }
 
-double Sk_WR(double q, double L, double b);
 double Sk_WR(double q, double L, double b)
 {
Index: sasmodels/models/onion.c
===================================================================
--- sasmodels/models/onion.c	(revision abdd01cc0315960ae6e5951909d2866515e37bae)
+++ sasmodels/models/onion.c	(revision abdd01cc0315960ae6e5951909d2866515e37bae)
@@ -4,5 +4,5 @@
     double thickness, double A)
 {
-  const double vol = 4.0/3.0 * M_PI * r * r * r;
+  const double vol = M_4PI_3 * cube(r);
   const double qr = q * r;
   const double alpha = A * r/thickness;
@@ -19,5 +19,5 @@
     double sinqr, cosqr;
     SINCOS(qr, sinqr, cosqr);
-    fun = -3.0(
+    fun = -3.0*(
             ((alphasq - qrsq)*sinqr/qr - 2.0*alpha*cosqr) / sumsq
                 - (alpha*sinqr/qr - cosqr)
@@ -32,5 +32,5 @@
 f_linear(double q, double r, double sld, double slope)
 {
-  const double vol = 4.0/3.0 * M_PI * r * r * r;
+  const double vol = M_4PI_3 * cube(r);
   const double qr = q * r;
   const double bes = sph_j1c(qr);
@@ -52,5 +52,5 @@
 {
   const double bes = sph_j1c(q * r);
-  const double vol = 4.0/3.0 * M_PI * r * r * r;
+  const double vol = M_4PI_3 * cube(r);
   return sld * vol * bes;
 }
@@ -64,15 +64,15 @@
     r += thickness[i];
   }
-  return 4.0/3.0 * M_PI * r * r * r;
+  return M_4PI_3*cube(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[],
+Iq(double q, double sld_core, double core_radius, double sld_solvent,
+    double n, double sld_in[], double sld_out[], double thickness[],
     double A[])
 {
   int i;
-  r = core_radius;
-  f = f_constant(q, r, core_sld);
+  double r = core_radius;
+  double f = f_constant(q, r, sld_core);
   for (i=0; i<n; i++){
     const double r0 = r;
@@ -92,6 +92,6 @@
     }
   }
-  f -= f_constant(q, r, solvent_sld);
-  f2 = f * f * 1.0e-4;
+  f -= f_constant(q, r, sld_solvent);
+  const double f2 = f * f * 1.0e-4;
 
   return f2;
Index: sasmodels/models/rpa.c
===================================================================
--- sasmodels/models/rpa.c	(revision abdd01cc0315960ae6e5951909d2866515e37bae)
+++ sasmodels/models/rpa.c	(revision abdd01cc0315960ae6e5951909d2866515e37bae)
@@ -1,32 +1,10 @@
 double Iq(double q, double case_num,
-    double Na, double Phia, double va, double a_sld, double ba,
-    double Nb, double Phib, double vb, double b_sld, double bb,
-    double Nc, double Phic, double vc, double c_sld, double bc,
-    double Nd, double Phid, double vd, double d_sld, double bd,
+    double N[], double Phi[], double v[], double L[], double b[],
     double Kab, double Kac, double Kad,
     double Kbc, double Kbd, double Kcd
     );
 
-double Iqxy(double qx, double qy, double case_num,
-    double Na, double Phia, double va, double a_sld, double ba,
-    double Nb, double Phib, double vb, double b_sld, double bb,
-    double Nc, double Phic, double vc, double c_sld, double bc,
-    double Nd, double Phid, double vd, double d_sld, double bd,
-    double Kab, double Kac, double Kad,
-    double Kbc, double Kbd, double Kcd
-    );
-
-double form_volume(void);
-
-double form_volume(void)
-{
-    return 1.0;
-}
-
 double Iq(double q, double case_num,
-    double Na, double Phia, double va, double La, double ba,
-    double Nb, double Phib, double vb, double Lb, double bb,
-    double Nc, double Phic, double vc, double Lc, double bc,
-    double Nd, double Phid, double vd, double Ld, double bd,
+    double N[], double Phi[], double v[], double L[], double b[],
     double Kab, double Kac, double Kad,
     double Kbc, double Kbd, double Kcd
@@ -36,6 +14,6 @@
 #if 0  // Sasview defaults
   if (icase <= 1) {
-    Na=Nb=1000.0;
-    Phia=Phib=0.0000001;
+    N[0]=N[1]=1000.0;
+    Phi[0]=Phi[1]=0.0000001;
     Kab=Kac=Kad=Kbc=Kbd=-0.0004;
     La=Lb=1.0e-12;
@@ -43,5 +21,5 @@
     ba=bb=5.0;
   } else if (icase <= 4) {
-    Phia=0.0000001;
+    Phi[0]=0.0000001;
     Kab=Kac=Kad=-0.0004;
     La=1.0e-12;
@@ -51,24 +29,24 @@
 #else
   if (icase <= 1) {
-    Na=Nb=0.0;
-    Phia=Phib=0.0;
+    N[0]=N[1]=0.0;
+    Phi[0]=Phi[1]=0.0;
     Kab=Kac=Kad=Kbc=Kbd=0.0;
-    La=Lb=Ld;
-    va=vb=vd;
-    ba=bb=0.0;
+    L[0]=L[1]=L[3];
+    v[0]=v[1]=v[3];
+    b[0]=b[1]=0.0;
   } else if (icase <= 4) {
-    Na = 0.0;
-    Phia=0.0;
+    N[0] = 0.0;
+    Phi[0]=0.0;
     Kab=Kac=Kad=0.0;
-    La=Ld;
-    va=vd;
-    ba=0.0;
+    L[0]=L[3];
+    v[0]=v[3];
+    b[0]=0.0;
   }
 #endif
 
-  const double Xa = q*q*ba*ba*Na/6.0;
-  const double Xb = q*q*bb*bb*Nb/6.0;
-  const double Xc = q*q*bc*bc*Nc/6.0;
-  const double Xd = q*q*bd*bd*Nd/6.0;
+  const double Xa = q*q*b[0]*b[0]*N[0]/6.0;
+  const double Xb = q*q*b[1]*b[1]*N[1]/6.0;
+  const double Xc = q*q*b[2]*b[2]*N[2]/6.0;
+  const double Xd = q*q*b[3]*b[3]*N[3]/6.0;
 
   // limit as Xa goes to 0 is 1
@@ -98,33 +76,33 @@
 #if 0
   const double S0aa = icase<5
-                      ? 1.0 : Na*Phia*va*Paa;
+                      ? 1.0 : N[0]*Phi[0]*v[0]*Paa;
   const double S0bb = icase<2
-                      ? 1.0 : Nb*Phib*vb*Pbb;
-  const double S0cc = Nc*Phic*vc*Pcc;
-  const double S0dd = Nd*Phid*vd*Pdd;
+                      ? 1.0 : N[1]*Phi[1]*v[1]*Pbb;
+  const double S0cc = N[2]*Phi[2]*v[2]*Pcc;
+  const double S0dd = N[3]*Phi[3]*v[3]*Pdd;
   const double S0ab = icase<8
-                      ? 0.0 : sqrt(Na*va*Phia*Nb*vb*Phib)*Pa*Pb;
+                      ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[1]*v[1]*Phi[1])*Pa*Pb;
   const double S0ac = icase<9
-                      ? 0.0 : sqrt(Na*va*Phia*Nc*vc*Phic)*Pa*Pc*exp(-Xb);
+                      ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[2]*v[2]*Phi[2])*Pa*Pc*exp(-Xb);
   const double S0ad = icase<9
-                      ? 0.0 : sqrt(Na*va*Phia*Nd*vd*Phid)*Pa*Pd*exp(-Xb-Xc);
+                      ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[3]*v[3]*Phi[3])*Pa*Pd*exp(-Xb-Xc);
   const double S0bc = (icase!=4 && icase!=7 && icase!= 9)
-                      ? 0.0 : sqrt(Nb*vb*Phib*Nc*vc*Phic)*Pb*Pc;
+                      ? 0.0 : sqrt(N[1]*v[1]*Phi[1]*N[2]*v[2]*Phi[2])*Pb*Pc;
   const double S0bd = (icase!=4 && icase!=7 && icase!= 9)
-                      ? 0.0 : sqrt(Nb*vb*Phib*Nd*vd*Phid)*Pb*Pd*exp(-Xc);
+                      ? 0.0 : sqrt(N[1]*v[1]*Phi[1]*N[3]*v[3]*Phi[3])*Pb*Pd*exp(-Xc);
   const double S0cd = (icase==0 || icase==2 || icase==5)
-                      ? 0.0 : sqrt(Nc*vc*Phic*Nd*vd*Phid)*Pc*Pd;
+                      ? 0.0 : sqrt(N[2]*v[2]*Phi[2]*N[3]*v[3]*Phi[3])*Pc*Pd;
 #else  // sasview equivalent
-//printf("Xc=%g, S0cc=%g*%g*%g*%g\n",Xc,Nc,Phic,vc,Pcc);
-  double S0aa = Na*Phia*va*Paa;
-  double S0bb = Nb*Phib*vb*Pbb;
-  double S0cc = Nc*Phic*vc*Pcc;
-  double S0dd = Nd*Phid*vd*Pdd;
-  double S0ab = sqrt(Na*va*Phia*Nb*vb*Phib)*Pa*Pb;
-  double S0ac = sqrt(Na*va*Phia*Nc*vc*Phic)*Pa*Pc*exp(-Xb);
-  double S0ad = sqrt(Na*va*Phia*Nd*vd*Phid)*Pa*Pd*exp(-Xb-Xc);
-  double S0bc = sqrt(Nb*vb*Phib*Nc*vc*Phic)*Pb*Pc;
-  double S0bd = sqrt(Nb*vb*Phib*Nd*vd*Phid)*Pb*Pd*exp(-Xc);
-  double S0cd = sqrt(Nc*vc*Phic*Nd*vd*Phid)*Pc*Pd;
+//printf("Xc=%g, S0cc=%g*%g*%g*%g\n",Xc,N[2],Phi[2],v[2],Pcc);
+  double S0aa = N[0]*Phi[0]*v[0]*Paa;
+  double S0bb = N[1]*Phi[1]*v[1]*Pbb;
+  double S0cc = N[2]*Phi[2]*v[2]*Pcc;
+  double S0dd = N[3]*Phi[3]*v[3]*Pdd;
+  double S0ab = sqrt(N[0]*v[0]*Phi[0]*N[1]*v[1]*Phi[1])*Pa*Pb;
+  double S0ac = sqrt(N[0]*v[0]*Phi[0]*N[2]*v[2]*Phi[2])*Pa*Pc*exp(-Xb);
+  double S0ad = sqrt(N[0]*v[0]*Phi[0]*N[3]*v[3]*Phi[3])*Pa*Pd*exp(-Xb-Xc);
+  double S0bc = sqrt(N[1]*v[1]*Phi[1]*N[2]*v[2]*Phi[2])*Pb*Pc;
+  double S0bd = sqrt(N[1]*v[1]*Phi[1]*N[3]*v[3]*Phi[3])*Pb*Pd*exp(-Xc);
+  double S0cd = sqrt(N[2]*v[2]*Phi[2]*N[3]*v[3]*Phi[3])*Pc*Pd;
 switch(icase){
   case 0:
@@ -311,7 +289,7 @@
   // Note: 1e-13 to convert from fm to cm for scattering length
   const double sqrt_Nav=sqrt(6.022045e+23) * 1.0e-13;
-  const double Lad = icase<5 ? 0.0 : (La/va - Ld/vd)*sqrt_Nav;
-  const double Lbd = icase<2 ? 0.0 : (Lb/vb - Ld/vd)*sqrt_Nav;
-  const double Lcd = (Lc/vc - Ld/vd)*sqrt_Nav;
+  const double Lad = icase<5 ? 0.0 : (L[0]/v[0] - L[3]/v[3])*sqrt_Nav;
+  const double Lbd = icase<2 ? 0.0 : (L[1]/v[1] - L[3]/v[3])*sqrt_Nav;
+  const double Lcd = (L[2]/v[2] - L[3]/v[3])*sqrt_Nav;
 
   const double result=Lad*Lad*S11 + Lbd*Lbd*S22 + Lcd*Lcd*S33
@@ -321,23 +299,2 @@
 
 }
-
-double Iqxy(double qx, double qy,
-    double case_num,
-    double Na, double Phia, double va, double a_sld, double ba,
-    double Nb, double Phib, double vb, double b_sld, double bb,
-    double Nc, double Phic, double vc, double c_sld, double bc,
-    double Nd, double Phid, double vd, double d_sld, double bd,
-    double Kab, double Kac, double Kad,
-    double Kbc, double Kbd, double Kcd
-    )
-{
-    double q = sqrt(qx*qx + qy*qy);
-    return Iq(q,
-        case_num,
-        Na, Phia, va, a_sld, ba,
-        Nb, Phib, vb, b_sld, bb,
-        Nc, Phic, vc, c_sld, bc,
-        Nd, Phid, vd, d_sld, bd,
-        Kab, Kac, Kad,
-        Kbc, Kbd, Kcd);
-}
Index: sasmodels/models/rpa.py
===================================================================
--- sasmodels/models/rpa.py	(revision ec45c4f5ff6e17014d968da6eb80308c2c8791f0)
+++ sasmodels/models/rpa.py	(revision a5b8477c53695e7df03ecec7346f0c1bd770d3f7)
@@ -86,29 +86,11 @@
 #   ["name", "units", default, [lower, upper], "type","description"],
 parameters = [
-    ["case_num", CASES, 0, [0, 10], "", "Component organization"],
+    ["case_num", "", 1, [CASES], "", "Component organization"],
 
-    ["Na", "", 1000.0, [1, inf], "", "Degree of polymerization"],
-    ["Phia", "", 0.25, [0, 1], "", "volume fraction"],
-    ["va", "mL/mol", 100.0, [0, inf], "", "specific volume"],
-    ["La", "fm", 10.0, [-inf, inf], "", "scattering length"],
-    ["ba", "Ang", 5.0, [0, inf], "", "segment length"],
-
-    ["Nb", "", 1000.0, [1, inf], "", "Degree of polymerization"],
-    ["Phib", "", 0.25, [0, 1], "", "volume fraction"],
-    ["vb", "mL/mol", 100.0, [0, inf], "", "specific volume"],
-    ["Lb", "fm", 10.0, [-inf, inf], "", "scattering length"],
-    ["bb", "Ang", 5.0, [0, inf], "", "segment length"],
-
-    ["Nc", "", 1000.0, [1, inf], "", "Degree of polymerization"],
-    ["Phic", "", 0.25, [0, 1], "", "volume fraction"],
-    ["vc", "mL/mol", 100.0, [0, inf], "", "specific volume"],
-    ["Lc", "fm", 10.0, [-inf, inf], "", "scattering length"],
-    ["bc", "Ang", 5.0, [0, inf], "", "segment length"],
-
-    ["Nd", "", 1000.0, [1, inf], "", "Degree of polymerization"],
-    ["Phid", "", 0.25, [0, 1], "", "volume fraction"],
-    ["vd", "mL/mol", 100.0, [0, inf], "", "specific volume"],
-    ["Ld", "fm", 10.0, [-inf, inf], "", "scattering length"],
-    ["bd", "Ang", 5.0, [0, inf], "", "segment length"],
+    ["N[4]", "", 1000.0, [1, inf], "", "Degree of polymerization"],
+    ["Phi[4]", "", 0.25, [0, 1], "", "volume fraction"],
+    ["v[4]", "mL/mol", 100.0, [0, inf], "", "specific volume"],
+    ["L[4]", "fm", 10.0, [-inf, inf], "", "scattering length"],
+    ["b[4]", "Ang", 5.0, [0, inf], "", "segment length"],
 
     ["Kab", "", -0.0004, [-inf, inf], "", "Interaction parameter"],
Index: sasmodels/models/spherical_sld.c
===================================================================
--- sasmodels/models/spherical_sld.c	(revision dd71228a6caa11bbe84bc1a9c3e340b533f71059)
+++ sasmodels/models/spherical_sld.c	(revision dd71228a6caa11bbe84bc1a9c3e340b533f71059)
@@ -0,0 +1,218 @@
+static double form_volume(
+    int n_shells,
+    double radius_core,
+    double thick_inter0,
+    double thick_flat[],
+    double thick_inter[])
+{
+    int i;
+    double r = radius_core;
+    r += thick_inter0;
+    for (i=0; i < n_shells; i++) {
+        r += thick_inter[i];
+        r += thick_flat[i];
+    }
+    return M_4PI_3*cube(r);
+}
+
+
+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 onion 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;
+            }
+        }
+    }
+
+    f2 = f * f * 1.0e-4;
+    return (f2);
+}
+
+
+/**
+ * Function to evaluate 1D SphereSLD function
+ * @param q: q-value
+ * @return: function value
+ */
+static double Iq(double q,
+    int n_shells,
+    int npts_inter,
+    double radius_core,
+    double sld_core,
+    double sld_solvent,
+    double func_inter0,
+    double thick_inter0,
+    double nu_inter0,
+    double sld_flat[],
+    double thick_flat[],
+    double func_inter[],
+    double thick_inter[],
+    double nu_inter[] ) {
+
+    double 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);
+
+    return intensity;
+}
+
+/**
+ * Function to evaluate 2D SphereSLD function
+ * @param q_x: value of Q along x
+ * @param q_y: value of Q along y
+ * @return: function value
+ */
+
+/*static double Iqxy(double qx, double qy,
+    int n_shells,
+    int npts_inter,
+    double radius_core
+    double sld_core,
+    double sld_solvent,
+    double sld_flat[],
+    double thick_flat[],
+    double func_inter[],
+    double thick_inter[],
+    double nu_inter[],
+    ) {
+
+    double q = sqrt(qx*qx + qy*qy);
+    return Iq(q, n_shells, npts_inter, radius_core, sld_core, sld_solvent,
+    sld_flat[], thick_flat[], func_inter[], thick_inter[], nu_inter[])
+}*/
+
Index: sasmodels/models/squarewell.py
===================================================================
--- sasmodels/models/squarewell.py	(revision aa9e4e3d3288cd9e9d890f59aa6aa6fcdcffc2d6)
+++ sasmodels/models/squarewell.py	(revision aa9e4e3d3288cd9e9d890f59aa6aa6fcdcffc2d6)
@@ -123,8 +123,4 @@
 """
 
-Iqxy = """
-    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);
-    """
-
 # ER defaults to 0.0
 # VR defaults to 1.0
Index: sasmodels/models/stickyhardsphere.py
===================================================================
--- sasmodels/models/stickyhardsphere.py	(revision ec45c4f5ff6e17014d968da6eb80308c2c8791f0)
+++ sasmodels/models/stickyhardsphere.py	(revision d2bb6049fb1ca2c1d60bcd3853b1d2ae889866b2)
@@ -171,8 +171,4 @@
 """
 
-Iqxy = """
-    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);
-    """
-
 # ER defaults to 0.0
 # VR defaults to 1.0
