Index: smodels/models/_core_multi_shell.py
===================================================================
--- sasmodels/models/_core_multi_shell.py	(revision f22e50e7d1bc8ecbd1dc21b3650d30b579aca290)
+++ 	(revision )
@@ -1,195 +1,0 @@
-r"""
-Definition
-----------
-
-This model is a trivial extension of the CoreShell function to a larger number
-of shells. The scattering length density profile for the default sld values 
-(w/ 4 shells).
-
-.. figure:: img/core_multi_shell_sld_default_profile.jpg
-
-    SLD profile of the core_multi_shell object from the center of sphere out
-    for the default SLDs.*
-
-The 2D scattering intensity is the same as *P(q)* above, regardless of the
-orientation of the *q* vector which is defined as
-
-.. math::
-
-    q = \sqrt{q_x^2 + q_y^2}
-
-.. note:: **Be careful!** The SLDs and scale can be highly correlated. Hold as
-         many of these parameters fixed as possible.
-
-.. note:: The outer most radius (= *radius* + *thickness*) is used as the
-          effective radius for *S(Q)* when *P(Q)* \* *S(Q)* is applied.
-
-For information about polarised and magnetic scattering, see 
-the :doc:`magnetic help <../sasgui/perspectives/fitting/mag_help>` documentation.
-
-Our model uses the form factor calculations implemented in a c-library provided
-by the NIST Center for Neutron Research (Kline, 2006).
-
-References
-----------
-See the :ref:`core_shell_sphere <core_shell_sphere>` model documentation.
-
-L A Feigin and D I Svergun,
-*Structure Analysis by Small-Angle X-Ray and Neutron Scattering*,
-Plenum Press, New York, 1987.
-
-**Author:** NIST IGOR/DANSE **on:** pre 2010
-
-**Last Modified by:** in progress **on:** March 20, 2016
-
-**Last Reviewed by:** in progress **on:** March 20, 2016
-"""
-
-
-
-from __future__ import division
-
-import numpy as np
-from numpy import inf, nan
-from math import fabs, exp, expm1
-
-name = "core_multi_shell"
-title = "This model provides the scattering from a spherical core with 1 to 4 \
- concentric shell structures. The SLDs of the core and each shell are \
- individually specified."
-
-description = """\
-Form factor for a core muti-shell (up to 4) sphere normalized by the volume.
-Each shell can have a unique thickness and sld.
-
-	background:background,
-	rad_core0: radius of sphere(core)
-	thick_shell#:the thickness of the shell#
-	sld_core0: the SLD of the sphere
-	sld_solv: the SLD of the solvent
-	sld_shell: the SLD of the shell#
-	A_shell#: the coefficient in the exponential function
-	
-	
-    scale: 1.0 if data is on absolute scale
-    volfraction: volume fraction of spheres
-    radius: the radius of the core
-    sld: the SLD of the core
-    thick_shelli: the thickness of the i'th shell from the core
-    sld_shelli: the SLD of the i'th shell from the core
-    sld_solvent: the SLD of the solvent
-    background: incoherent background
-
-"""
-
-category = "shape:sphere"
-
-
-#             ["name", "units", default, [lower, upper], "type","description"],
-parameters = [["volfraction", "", 0.05, [0,1],"",
-               "volume fraction of spheres"],
-              ["sld_core", "1e-6/Ang^2", 1.0, [-inf, inf], "",
-               "Core scattering length density"],
-              ["radius", "Ang", 200., [0, inf], "volume",
-               "Radius of the core"],
-              ["sld_solvent", "1e-6/Ang^2", 6.4, [-inf, inf], "",
-               "Solvent scattering length density"],
-              ["n", "", 1, [0, 10], "volume",
-               "number of shells"],
-              ["sld[n]", "1e-6/Ang^2", 1.7, [-inf, inf], "",
-               "scattering length density of shell k"],
-              ["thickness[n]", "Ang", 40., [0, inf], "volume",
-               "Thickness 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 profile(volfraction, sld_core, radius, sld_solvent, n, sld, thickness):
-    """
-    SLD profile
-    
-    :return: (r, beta) where r is a list of radius of the transition points\
-      and beta is a list of the corresponding SLD values.
-
-    """
-#        r = []
-#        beta = []
-#        # for core at r=0
-#        r.append(0)
-#        beta.append(self.params['sld_core0'])
-#        # for core at r=rad_core
-#        r.append(self.params['rad_core0'])
-#        beta.append(self.params['sld_core0'])
-#
-#        # for shells
-#        for n in range(1, self.n_shells+1):
-#            # Left side of each shells
-#            r0 = r[len(r)-1]
-#            r.append(r0)
-#            exec "beta.append(self.params['sld_shell%s'% str(n)])"
-#
-#            # Right side of each shells
-#            exec "r0 += self.params['thick_shell%s'% str(n)]"
-#            r.append(r0)
-#            exec "beta.append(self.params['sld_shell%s'% str(n)])"
-#
-#        # for solvent
-#        r0 = r[len(r)-1]            
-#        r.append(r0)
-#        beta.append(self.params['sld_solv'])
-#        r_solv = 5*r0/4
-#        r.append(r_solv)
-#        beta.append(self.params['sld_solv'])
-#
-#        return r, beta
-# above is directly from old code -- below is alotered from Woitek's first
-# cut an the onion.  Seems like we should be consistant?
-
-    total_radius = 1.25*(sum(thickness[:n]) + radius + 1)
-
-    r = []
-    beta = []
-
-    # add in the core
-    r.append(0)
-    beta.append(sld)
-    r.append(radius)
-    beta.append(sld)
-
-    # add in the shells
-    for k in range(n):
-        # Left side of each shells
-        r0 = r[-1]
-        r.append(r0)
-        beta.append(sld[k])
-        r.append(r0 + thickness[k])
-        beta.append(sld[k])
-   # add in the solvent
-    r.append(r[-1])
-    beta.append(sld_solvent)
-    r.append(total_radius)
-    beta.append(sld_solvent)
-
-    return np.asarray(r), np.asarray(beta)
-
-def ER(radius, n, thickness):
-    n = n[0]  # n cannot be polydisperse
-    return np.sum(thickness[:n], axis=0) + radius
-
-def VR(radius, n, thick_shell):
-    return 1.0, 1.0
-
-demo = dict(volfraction = 1.0,
-            sld = 6.4,
-            radius = 60,
-            sld_solvent = 6.4,
-            n = 1,
-            sld_shell = [2.0],
-            thick_shell = [10])
Index: smodels/models/_onion.py
===================================================================
--- sasmodels/models/_onion.py	(revision f22e50e7d1bc8ecbd1dc21b3650d30b579aca290)
+++ 	(revision )
@@ -1,385 +1,0 @@
-r"""
-This model provides the form factor, $P(q)$, for a multi-shell sphere where
-the scattering length density (SLD) of the each shell is described by an
-exponential, linear, or constant function. The form factor is normalized by
-the volume of the sphere where the SLD is not identical to the SLD of the
-solvent. We currently provide up to 9 shells with this model.
-
-NB: *radius* represents the core radius $r_0$ and
-*thickness[k]* represents the thickness of the shell, $r_{k+1} - r_k$.
-
-Definition
-----------
-
-The 1D scattering intensity is calculated in the following way
-
-.. math::
-
-    P(q) &= [f]^2 / V_\text{particle}
-
-where
-
-.. math::
-
-    f    &= f_\text{core}
-            + \left(\sum_{\text{shell}=1}^N f_\text{shell}\right)
-            + f_\text{solvent}
-
-
-The shells are spherically symmetric with particle density $\rho(r)$ and
-constant SLD within the core and solvent, so
-
-.. math::
-
-    f_\text{core}
-        &= 4\pi\int_0^{r_\text{core}} \rho_\text{core}
-            \frac{\sin(qr)}{qr}\, r^2\,\mathrm{d}r
-        &= 3\rho_\text{core} V(r_\text{core})
-            \frac{j_1(qr_\text{core})}{qr_\text{core}} \\
-    f_\text{shell}
-        &= 4\pi\int_{r_{\text{shell}-1}}^{r_\text{shell}}
-            \rho_\text{shell}(r)\frac{\sin(qr)}{qr}\,r^2\,\mathrm{d}r \\
-    f_\text{solvent}
-        &= 4\pi\int_{r_N}^\infty
-            \rho_\text{solvent}\frac{\sin(qr)}{qr}\,r^2\,\mathrm{d}r
-        &= -3\rho_\text{solvent}V(r_N)\frac{j_1(q r_N)}{q r_N}
-
-where the spherical bessel function $j_1$ is
-
-.. math::
-
-    j_1(x) = \frac{\sin(x)}{x^2} - \frac{\cos(x)}{x}
-
-and the volume is $V(r) = \frac{4\pi}{3}r^3$. The volume of the particle
-is determined by the radius of the outer shell, so $V_\text{particle} = V(r_N)$.
-
-Now lets consider the SLD of a shell defined by
-
-.. math::
-
-    \rho_\text{shell}(r) = \begin{cases}
-        B\exp\left(A(r-r_{\text{shell}-1})/\Delta t_\text{shell}\right)
-            + C & \mbox{for } A \neq 0 \\
-        \rho_\text{in} = \text{constant} & \mbox{for } A = 0
-    \end{cases}
-
-An example of a possible SLD profile is shown below where
-$\rho_\text{in}$ and $\Delta t_\text{shell}$ stand for the
-SLD of the inner side of the $k^\text{th}$ shell and the
-thickness of the $k^\text{th}$ shell in the equation above, respectively.
-
-For $A \gt 0$,
-
-.. math::
-
-    f_\text{shell} &= 4 \pi \int_{r_{\text{shell}-1}}^{r_\text{shell}}
-        \left[ B\exp
-            \left(A (r - r_{\text{shell}-1}) / \Delta t_\text{shell} \right) + C
-        \right] \frac{\sin(qr)}{qr}\,r^2\,\mathrm{d}r \\
-    &= 3BV(r_\text{shell}) e^A h(\alpha_\text{out},\beta_\text{out})
-        - 3BV(r_{\text{shell}-1}) h(\alpha_\text{in},\beta_\text{in})
-        + 3CV(r_{\text{shell}}) \frac{j_1(\beta_\text{out})}{\beta_\text{out}}
-        - 3CV(r_{\text{shell}-1}) \frac{j_1(\beta_\text{in})}{\beta_\text{in}}
-
-for
-
-.. math::
-    :nowrap:
-
-    \begin{align*}
-    B&=\frac{\rho_\text{out} - \rho_\text{in}}{e^A-1}
-         &C &= \frac{\rho_\text{in}e^A - \rho_\text{out}}{e^A-1} \\
-    \alpha_\text{in} &= A\frac{r_{\text{shell}-1}}{\Delta t_\text{shell}}
-         &\alpha_\text{out} &= A\frac{r_\text{shell}}{\Delta t_\text{shell}} \\
-    \beta_\text{in} &= qr_{\text{shell}-1}
-        &\beta_\text{out} &= qr_\text{shell}
-    \end{align*}
-
-where $h$ is
-
- .. math::
-
-    h(x,y) = \frac{x \sin(y) - y\cos(y)}{(x^2+y^2)y}
-               - \frac{(x^2-y^2)\sin(y) - 2xy\cos(y)}{(x^2+y^2)^2y}
-
-
-For $A \sim 0$, e.g., $A = -0.0001$, this function converges to that of the
-linear SLD profile with
-$\rho_\text{shell}(r) \approx A(r-r_{\text{shell}-1})/\Delta t_\text{shell})+B$,
-so this case is equivalent to
-
-.. math::
-
-
-    f_\text{shell}
-    &=
-      3 V(r_\text{shell}) \frac{\Delta\rho_\text{shell}}{\Delta t_\text{shell}}
-        \left[\frac{
-                2 \cos(qr_\text{out})
-                    + qr_\text{out} \sin(qr_\text{out})
-            }{
-                (qr_\text{out})^4
-            }\right] \\
-     &{}
-      -3 V(r_\text{shell}) \frac{\Delta\rho_\text{shell}}{\Delta t_\text{shell}}
-        \left[\frac{
-                    2\cos(qr_\text{in})
-                +qr_\text{in}\sin(qr_\text{in})
-            }{
-                (qr_\text{in})^4
-            }\right] \\
-    &{}
-      +3\rho_\text{out}V(r_\text{shell}) \frac{j_1(qr_\text{out})}{qr_\text{out}}
-      -3\rho_\text{in}V(r_{\text{shell}-1}) \frac{j_1(qr_\text{in})}{qr_\text{in}}
-
-For $A = 0$, the exponential function has no dependence on the radius (so that
-$\rho_\text{out}$ is ignored this case) and becomes flat. We set the constant
-to $\rho_\text{in}$ for convenience, and thus the form factor contributed by
-the shells is
-
-.. math::
-
-    f_\text{shell} =
-        3\rho_\text{in}V(r_\text{shell})
-           \frac{j_1(qr_\text{out})}{qr_\text{out}}
-        - 3\rho_\text{in}V(r_{\text{shell}-1})
-            \frac{j_1(qr_\text{in})}{qr_\text{in}}
-
-.. figure:: img/onion_geometry.png
-
-    Example of an onion model profile.
-
-The 2D scattering intensity is the same as $P(q)$ above, regardless of the
-orientation of the $q$ vector which is defined as
-
-.. math::
-
-    q = \sqrt{q_x^2 + q_y^2}
-
-NB: 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.
-"""
-
-#
-# Give a polynomial $\rho(r) = Ar^3 + Br^2 + Cr + D$ for density,
-#
-# .. math::
-#
-#    f = 4 \pi \int_a^b \rho(r) \sin(qr)/(qr) \mathrm{d}r  = h(b) - h(a)
-#
-# where
-#
-# .. math::
-#
-#    h(r) = \frac{4 \pi}{q^6}\left[
-#        (q^3(4Ar^3 + 3Br^2 + 2Cr + D) - q(24Ar + 6B)) \sin(qr)
-#      - (q^4(Ar^4 + Br^3 + Cr^2 + Dr) - q^2(12Ar^2 + 6Br + 2C) + 24A) \cos(qr)
-#    \right]
-#
-# Use the monotonic spline to get the polynomial coefficients for each shell.
-#
-# Order 0
-#
-# .. math::
-#
-#    h(r) = \frac{4 \pi}{q^3} \left[
-#       - \cos(qr) (Ar) q
-#       + \sin(qr) (A)
-#    \right]
-#
-# Order 1
-#
-# .. math::
-#
-#   h(r) = \frac{4 \pi}{q^4} \left[
-#       - \cos(qr) ( Ar^2 + Br) q^2
-#       + \sin(qr) ( Ar   + B ) q
-#       + \cos(qr) (2A        )
-#   \right]
-#
-# Order 2
-#
-# .. math::
-#  h(r) = \frac{4 \pi}{q^5} \left[
-#        - \cos(qr) ( Ar^3 +  Br^2 + Cr) q^3
-#        + \sin(qr) (3Ar^2 + 2Br   + C ) q^2
-#        + \cos(qr) (6Ar   + 2B        ) q
-#        - \sin(qr) (6A                )
-#
-# Order 3
-#
-#    h(r) = \frac{4 \pi}{q^6}\left[
-#      - \cos(qr) (  Ar^4 +  Br^3 +  Cr^2 + Dr) q^4
-#      + \sin(qr) ( 4Ar^3 + 3Br^2 + 2Cr   + D ) q^3
-#      + \cos(qr) (12Ar^2 + 6Br   + 2C        ) q^2
-#      - \sin(qr) (24Ar   + 6B                ) q
-#      - \cos(qr) (24A                        )
-#    \right]
-#
-# Order p
-#
-#    h(r) = \frac{4 \pi}{q^{2}}
-#      \sum_{k=0}^p -\frac{d^k\cos(qr)}{dr^k} \frac{d^k r\rho(r)}{dr^k} (qr)^{-k}
-#
-# Given the equation
-#
-#    f = sum_(k=0)^(n-1) h_k(r_(k+1)) - h_k(r_k)
-#
-# we can rearrange the terms so that
-#
-#    f = sum_0^(n-1) h_k(r_(k+1)) - sum_0^(n-1) h_k(r_k)
-#      = sum_1^n h_(k-1)(r_k) - sum_0^(n-1) h_k(r_k)
-#      = h_(n-1)(r_n) - h_0(r_0) + sum_1^(n-1) [h_(k-1)(r_k) - h_k(r_k)]
-#      = h_(n-1)(r_n) - h_0(r_0) - sum_1^(n-1) h_(Delta k)(r_k)
-#
-# where
-#
-#    h_(Delta k)(r) = h(Delta rho_k, r)
-#
-# for
-#
-#    Delta rho_k = (A_k-A_(k-1)) r^p + (B_k-B_(k-1)) r^(p-1) + ...
-#
-# Using l'H\^opital's Rule 6 times on the order 3 polynomial,
-#
-#   lim_(q->0) h(r) = (140D r^3 + 180C r^4 + 144B r^5 + 120A r^6)/720
-#
-
-
-from __future__ import division
-
-import numpy as np
-from numpy import inf, nan
-from math import fabs, exp, expm1
-
-name = "onion"
-title = "Onion shell model with constant, linear or exponential density"
-
-description = """\
-Form factor of mutishells normalized by the volume. Here each shell is
-described by an exponential function;
-
-	I) For A_shell != 0,
-		f(r) = B*exp(A_shell*(r-r_in)/thick_shell)+C
-	where
-		B=(sld_out-sld_in)/(exp(A_shell)-1)
-		C=sld_in-B.
-	Note that in the above case, the function becomes a linear function
-	as A_shell --> 0+ or 0-.
-
-	II) For the exact point of A_shell == 0,
-		f(r) = sld_in ,i.e., it crosses over flat function
-	Note that the 'sld_out' becaomes NULL in this case.
-
-	background:background,
-	rad_core0: radius of sphere(core)
-	thick_shell#:the thickness of the shell#
-	sld_core0: the SLD of the sphere
-	sld_solv: the SLD of the solvent
-	sld_shell: the SLD of the shell#
-	A_shell#: the coefficient in the exponential function
-"""
-
-category = "shape:sphere"
-
-# TODO: n is a volume parameter that is not polydisperse
-
-#             ["name", "units", default, [lower, upper], "type","description"],
-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"],
-              ["sld_solvent", "1e-6/Ang^2", 6.4, [-inf, inf], "",
-               "Solvent scattering length density"],
-              ["n_shells", "", 1, [0, 10], "volume",
-               "number of shells"],
-              ["sld_in[n_shells]", "1e-6/Ang^2", 1.7, [-inf, inf], "",
-               "scattering length density at the inner radius of shell k"],
-              ["sld_out[n_shells]", "1e-6/Ang^2", 2.0, [-inf, inf], "",
-               "scattering length density at the outer radius of shell k"],
-              ["thickness[n_shells]", "Ang", 40., [0, inf], "volume",
-               "Thickness of shell k"],
-              ["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
-
-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):
-    """
-    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
-
-    r = []
-    beta = []
-
-    # add in the core
-    r.append(0)
-    beta.append(core_sld)
-    r.append(core_radius)
-    beta.append(core_sld)
-
-    # add in the shells
-    for k in range(n_shells):
-        # Left side of each shells
-        r0 = r[-1]
-        r.append(r0)
-        beta.append(in_sld[k])
-
-        if fabs(A[k]) < 1.0e-16:
-            # flat shell
-            r.append(r0 + thickness[k])
-            beta.append(out_sld[k])
-        else:
-            # exponential shell
-            # num_steps must be at least 1, so use floor()+1 rather than ceil
-            # to protect against a thickness0.
-            num_steps = np.floor(thickness[k]/dr) + 1
-            slope = (out_sld[k] - in_sld[k])/expm1(A[k])
-            const = (in_sld[k] - slope)
-            for rk in np.linspace(0, thickness[k], num_steps+1):
-                r.append(r0+rk)
-                beta.append(slope*exp(A[k]*rk/thickness[k]) + const)
-
-    # add in the solvent
-    r.append(r[-1])
-    beta.append(solvent_sld)
-    r.append(total_radius)
-    beta.append(solvent_sld)
-
-    return np.asarray(r), 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
-
-demo = {
-    "sld_solvent": 2.2,
-    "sld_core": 1.0,
-    "core_radius": 100,
-    "n": 4,
-    "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],
-    # Could also specify them individually as
-    # "A[1]": 0, "A[2]": -1, "A[3]": 1e-4, "A[4]": 1,
-    }
-
Index: smodels/models/_spherical_sld.py
===================================================================
--- sasmodels/models/_spherical_sld.py	(revision f22e50e7d1bc8ecbd1dc21b3650d30b579aca290)
+++ 	(revision )
@@ -1,322 +1,0 @@
-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
-
-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
-----------
-
-The form factor $P(q)$ in 1D is calculated by:
-
-.. math::
-
-    P(q) = \frac{f^2}{V_\text{particle}} \text{ where }
-    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}
-
-For a spherically symmetric particle with a particle density $\rho_x(r)$
-the sld function can be defined as:
-
-.. math::
-
-    f_x = 4 \pi \int_{0}^{\infty} \rho_x(r)  \frac{\sin(qr)} {qr^2} r^2 dr
-
-
-so that individual terms can be calcualted as follows:
-
-.. math::
-    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 } ) ) }
-    {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 =
-    3 \rho_\text{solvent} V(r_N)
-    \Big[ \frac{\sin(qr_N) - qr_N \cos(qr_N)} {qr_N^3} \Big]
-
-
-Here we assumed that the SLDs of the core and solvent are constant against $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
-
-Exp:
-
-.. 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 \\
-    \end{cases}
-
-Power-Law
-
-.. 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 \\
-    \rho_{\text{flat}_{i+1}}  & \text{for} A = 0 \\
-    \end{cases}
-
-Erf:
-
-.. 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 \\
-    \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 =
-    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
-
-    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}) }
-    {\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}) }
-    {\beta_\text{in}^4 } \Big]
-
-    + 3 \rho_{ \text{inter}_i } ( r_{j+1} )  V ( r_j )
-    \Big[ \frac {\sin(\beta_\text{out}) - \cos(\beta_\text{out}) }
-    {\beta_\text{out}^4 } \Big]
-
-    - 3 \rho_{ \text{inter}_i } ( r_{j} )  V ( r_j )
-    \Big[ \frac {\sin(\beta_\text{in}) - \cos(\beta_\text{in}) }
-    {\beta_\text{in}^4 } \Big]
-    \Big]
-
-where
-
-.. math::
-    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}
-
-    \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$
-
-Finally form factor can be calculated by
-
-.. math::
-
-    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,
-where the $q$ vector is defined as
-
-.. math::
-
-    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.
-
-References
-----------
-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
-
-name = "spherical_sld"
-title = "Sperical SLD intensity calculation"
-description = """
-            I(q) =
-               background = Incoherent background [1/cm]
-        """
-category = "sphere-based"
-
-# pylint: disable=bad-whitespace, line-too-long
-#            ["name", "units", default, [lower, upper], "type", "description"],
-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"]
-
-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
-    [{"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/core_multi_shell.py
===================================================================
--- sasmodels/models/core_multi_shell.py	(revision 46ed76021e194fbd75cd7aeac6434e865025002f)
+++ sasmodels/models/core_multi_shell.py	(revision 46ed76021e194fbd75cd7aeac6434e865025002f)
@@ -0,0 +1,195 @@
+r"""
+Definition
+----------
+
+This model is a trivial extension of the CoreShell function to a larger number
+of shells. The scattering length density profile for the default sld values 
+(w/ 4 shells).
+
+.. figure:: img/core_multi_shell_sld_default_profile.jpg
+
+    SLD profile of the core_multi_shell object from the center of sphere out
+    for the default SLDs.*
+
+The 2D scattering intensity is the same as *P(q)* above, regardless of the
+orientation of the *q* vector which is defined as
+
+.. math::
+
+    q = \sqrt{q_x^2 + q_y^2}
+
+.. note:: **Be careful!** The SLDs and scale can be highly correlated. Hold as
+         many of these parameters fixed as possible.
+
+.. note:: The outer most radius (= *radius* + *thickness*) is used as the
+          effective radius for *S(Q)* when *P(Q)* \* *S(Q)* is applied.
+
+For information about polarised and magnetic scattering, see 
+the :doc:`magnetic help <../sasgui/perspectives/fitting/mag_help>` documentation.
+
+Our model uses the form factor calculations implemented in a c-library provided
+by the NIST Center for Neutron Research (Kline, 2006).
+
+References
+----------
+See the :ref:`core_shell_sphere <core_shell_sphere>` model documentation.
+
+L A Feigin and D I Svergun,
+*Structure Analysis by Small-Angle X-Ray and Neutron Scattering*,
+Plenum Press, New York, 1987.
+
+**Author:** NIST IGOR/DANSE **on:** pre 2010
+
+**Last Modified by:** in progress **on:** March 20, 2016
+
+**Last Reviewed by:** in progress **on:** March 20, 2016
+"""
+
+
+
+from __future__ import division
+
+import numpy as np
+from numpy import inf, nan
+from math import fabs, exp, expm1
+
+name = "core_multi_shell"
+title = "This model provides the scattering from a spherical core with 1 to 4 \
+ concentric shell structures. The SLDs of the core and each shell are \
+ individually specified."
+
+description = """\
+Form factor for a core muti-shell (up to 4) sphere normalized by the volume.
+Each shell can have a unique thickness and sld.
+
+	background:background,
+	rad_core0: radius of sphere(core)
+	thick_shell#:the thickness of the shell#
+	sld_core0: the SLD of the sphere
+	sld_solv: the SLD of the solvent
+	sld_shell: the SLD of the shell#
+	A_shell#: the coefficient in the exponential function
+	
+	
+    scale: 1.0 if data is on absolute scale
+    volfraction: volume fraction of spheres
+    radius: the radius of the core
+    sld: the SLD of the core
+    thick_shelli: the thickness of the i'th shell from the core
+    sld_shelli: the SLD of the i'th shell from the core
+    sld_solvent: the SLD of the solvent
+    background: incoherent background
+
+"""
+
+category = "shape:sphere"
+
+
+#             ["name", "units", default, [lower, upper], "type","description"],
+parameters = [["volfraction", "", 0.05, [0,1],"",
+               "volume fraction of spheres"],
+              ["sld_core", "1e-6/Ang^2", 1.0, [-inf, inf], "",
+               "Core scattering length density"],
+              ["radius", "Ang", 200., [0, inf], "volume",
+               "Radius of the core"],
+              ["sld_solvent", "1e-6/Ang^2", 6.4, [-inf, inf], "",
+               "Solvent scattering length density"],
+              ["n", "", 1, [0, 4], "volume",
+               "number of shells"],
+              ["sld[n]", "1e-6/Ang^2", 1.7, [-inf, inf], "",
+               "scattering length density of shell k"],
+              ["thickness[n]", "Ang", 40., [0, inf], "volume",
+               "Thickness 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 profile(volfraction, sld_core, radius, sld_solvent, n, sld, thickness):
+    """
+    SLD profile
+    
+    :return: (r, beta) where r is a list of radius of the transition points\
+      and beta is a list of the corresponding SLD values.
+
+    """
+#        r = []
+#        beta = []
+#        # for core at r=0
+#        r.append(0)
+#        beta.append(self.params['sld_core0'])
+#        # for core at r=rad_core
+#        r.append(self.params['rad_core0'])
+#        beta.append(self.params['sld_core0'])
+#
+#        # for shells
+#        for n in range(1, self.n_shells+1):
+#            # Left side of each shells
+#            r0 = r[len(r)-1]
+#            r.append(r0)
+#            exec "beta.append(self.params['sld_shell%s'% str(n)])"
+#
+#            # Right side of each shells
+#            exec "r0 += self.params['thick_shell%s'% str(n)]"
+#            r.append(r0)
+#            exec "beta.append(self.params['sld_shell%s'% str(n)])"
+#
+#        # for solvent
+#        r0 = r[len(r)-1]            
+#        r.append(r0)
+#        beta.append(self.params['sld_solv'])
+#        r_solv = 5*r0/4
+#        r.append(r_solv)
+#        beta.append(self.params['sld_solv'])
+#
+#        return r, beta
+# above is directly from old code -- below is alotered from Woitek's first
+# cut an the onion.  Seems like we should be consistant?
+
+    total_radius = 1.25*(sum(thickness[:n]) + radius + 1)
+
+    r = []
+    beta = []
+
+    # add in the core
+    r.append(0)
+    beta.append(sld)
+    r.append(radius)
+    beta.append(sld)
+
+    # add in the shells
+    for k in range(n):
+        # Left side of each shells
+        r0 = r[-1]
+        r.append(r0)
+        beta.append(sld[k])
+        r.append(r0 + thickness[k])
+        beta.append(sld[k])
+   # add in the solvent
+    r.append(r[-1])
+    beta.append(sld_solvent)
+    r.append(total_radius)
+    beta.append(sld_solvent)
+
+    return np.asarray(r), np.asarray(beta)
+
+def ER(radius, n, thickness):
+    n = n[0]  # n cannot be polydisperse
+    return np.sum(thickness[:n], axis=0) + radius
+
+def VR(radius, n, thickness):
+    return 1.0, 1.0
+
+demo = dict(volfraction = 1.0,
+            sld_core = 6.4,
+            radius = 60,
+            sld_solvent = 6.4,
+            n = 1,
+            sld = [2.0],
+            thickness = [10])
Index: sasmodels/models/onion.py
===================================================================
--- sasmodels/models/onion.py	(revision 46ed76021e194fbd75cd7aeac6434e865025002f)
+++ sasmodels/models/onion.py	(revision 46ed76021e194fbd75cd7aeac6434e865025002f)
@@ -0,0 +1,385 @@
+r"""
+This model provides the form factor, $P(q)$, for a multi-shell sphere where
+the scattering length density (SLD) of the each shell is described by an
+exponential, linear, or constant function. The form factor is normalized by
+the volume of the sphere where the SLD is not identical to the SLD of the
+solvent. We currently provide up to 9 shells with this model.
+
+NB: *radius* represents the core radius $r_0$ and
+*thickness[k]* represents the thickness of the shell, $r_{k+1} - r_k$.
+
+Definition
+----------
+
+The 1D scattering intensity is calculated in the following way
+
+.. math::
+
+    P(q) &= [f]^2 / V_\text{particle}
+
+where
+
+.. math::
+
+    f    &= f_\text{core}
+            + \left(\sum_{\text{shell}=1}^N f_\text{shell}\right)
+            + f_\text{solvent}
+
+
+The shells are spherically symmetric with particle density $\rho(r)$ and
+constant SLD within the core and solvent, so
+
+.. math::
+
+    f_\text{core}
+        &= 4\pi\int_0^{r_\text{core}} \rho_\text{core}
+            \frac{\sin(qr)}{qr}\, r^2\,\mathrm{d}r
+        &= 3\rho_\text{core} V(r_\text{core})
+            \frac{j_1(qr_\text{core})}{qr_\text{core}} \\
+    f_\text{shell}
+        &= 4\pi\int_{r_{\text{shell}-1}}^{r_\text{shell}}
+            \rho_\text{shell}(r)\frac{\sin(qr)}{qr}\,r^2\,\mathrm{d}r \\
+    f_\text{solvent}
+        &= 4\pi\int_{r_N}^\infty
+            \rho_\text{solvent}\frac{\sin(qr)}{qr}\,r^2\,\mathrm{d}r
+        &= -3\rho_\text{solvent}V(r_N)\frac{j_1(q r_N)}{q r_N}
+
+where the spherical bessel function $j_1$ is
+
+.. math::
+
+    j_1(x) = \frac{\sin(x)}{x^2} - \frac{\cos(x)}{x}
+
+and the volume is $V(r) = \frac{4\pi}{3}r^3$. The volume of the particle
+is determined by the radius of the outer shell, so $V_\text{particle} = V(r_N)$.
+
+Now lets consider the SLD of a shell defined by
+
+.. math::
+
+    \rho_\text{shell}(r) = \begin{cases}
+        B\exp\left(A(r-r_{\text{shell}-1})/\Delta t_\text{shell}\right)
+            + C & \mbox{for } A \neq 0 \\
+        \rho_\text{in} = \text{constant} & \mbox{for } A = 0
+    \end{cases}
+
+An example of a possible SLD profile is shown below where
+$\rho_\text{in}$ and $\Delta t_\text{shell}$ stand for the
+SLD of the inner side of the $k^\text{th}$ shell and the
+thickness of the $k^\text{th}$ shell in the equation above, respectively.
+
+For $A \gt 0$,
+
+.. math::
+
+    f_\text{shell} &= 4 \pi \int_{r_{\text{shell}-1}}^{r_\text{shell}}
+        \left[ B\exp
+            \left(A (r - r_{\text{shell}-1}) / \Delta t_\text{shell} \right) + C
+        \right] \frac{\sin(qr)}{qr}\,r^2\,\mathrm{d}r \\
+    &= 3BV(r_\text{shell}) e^A h(\alpha_\text{out},\beta_\text{out})
+        - 3BV(r_{\text{shell}-1}) h(\alpha_\text{in},\beta_\text{in})
+        + 3CV(r_{\text{shell}}) \frac{j_1(\beta_\text{out})}{\beta_\text{out}}
+        - 3CV(r_{\text{shell}-1}) \frac{j_1(\beta_\text{in})}{\beta_\text{in}}
+
+for
+
+.. math::
+    :nowrap:
+
+    \begin{align*}
+    B&=\frac{\rho_\text{out} - \rho_\text{in}}{e^A-1}
+         &C &= \frac{\rho_\text{in}e^A - \rho_\text{out}}{e^A-1} \\
+    \alpha_\text{in} &= A\frac{r_{\text{shell}-1}}{\Delta t_\text{shell}}
+         &\alpha_\text{out} &= A\frac{r_\text{shell}}{\Delta t_\text{shell}} \\
+    \beta_\text{in} &= qr_{\text{shell}-1}
+        &\beta_\text{out} &= qr_\text{shell}
+    \end{align*}
+
+where $h$ is
+
+ .. math::
+
+    h(x,y) = \frac{x \sin(y) - y\cos(y)}{(x^2+y^2)y}
+               - \frac{(x^2-y^2)\sin(y) - 2xy\cos(y)}{(x^2+y^2)^2y}
+
+
+For $A \sim 0$, e.g., $A = -0.0001$, this function converges to that of the
+linear SLD profile with
+$\rho_\text{shell}(r) \approx A(r-r_{\text{shell}-1})/\Delta t_\text{shell})+B$,
+so this case is equivalent to
+
+.. math::
+
+
+    f_\text{shell}
+    &=
+      3 V(r_\text{shell}) \frac{\Delta\rho_\text{shell}}{\Delta t_\text{shell}}
+        \left[\frac{
+                2 \cos(qr_\text{out})
+                    + qr_\text{out} \sin(qr_\text{out})
+            }{
+                (qr_\text{out})^4
+            }\right] \\
+     &{}
+      -3 V(r_\text{shell}) \frac{\Delta\rho_\text{shell}}{\Delta t_\text{shell}}
+        \left[\frac{
+                    2\cos(qr_\text{in})
+                +qr_\text{in}\sin(qr_\text{in})
+            }{
+                (qr_\text{in})^4
+            }\right] \\
+    &{}
+      +3\rho_\text{out}V(r_\text{shell}) \frac{j_1(qr_\text{out})}{qr_\text{out}}
+      -3\rho_\text{in}V(r_{\text{shell}-1}) \frac{j_1(qr_\text{in})}{qr_\text{in}}
+
+For $A = 0$, the exponential function has no dependence on the radius (so that
+$\rho_\text{out}$ is ignored this case) and becomes flat. We set the constant
+to $\rho_\text{in}$ for convenience, and thus the form factor contributed by
+the shells is
+
+.. math::
+
+    f_\text{shell} =
+        3\rho_\text{in}V(r_\text{shell})
+           \frac{j_1(qr_\text{out})}{qr_\text{out}}
+        - 3\rho_\text{in}V(r_{\text{shell}-1})
+            \frac{j_1(qr_\text{in})}{qr_\text{in}}
+
+.. figure:: img/onion_geometry.png
+
+    Example of an onion model profile.
+
+The 2D scattering intensity is the same as $P(q)$ above, regardless of the
+orientation of the $q$ vector which is defined as
+
+.. math::
+
+    q = \sqrt{q_x^2 + q_y^2}
+
+NB: 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.
+"""
+
+#
+# Give a polynomial $\rho(r) = Ar^3 + Br^2 + Cr + D$ for density,
+#
+# .. math::
+#
+#    f = 4 \pi \int_a^b \rho(r) \sin(qr)/(qr) \mathrm{d}r  = h(b) - h(a)
+#
+# where
+#
+# .. math::
+#
+#    h(r) = \frac{4 \pi}{q^6}\left[
+#        (q^3(4Ar^3 + 3Br^2 + 2Cr + D) - q(24Ar + 6B)) \sin(qr)
+#      - (q^4(Ar^4 + Br^3 + Cr^2 + Dr) - q^2(12Ar^2 + 6Br + 2C) + 24A) \cos(qr)
+#    \right]
+#
+# Use the monotonic spline to get the polynomial coefficients for each shell.
+#
+# Order 0
+#
+# .. math::
+#
+#    h(r) = \frac{4 \pi}{q^3} \left[
+#       - \cos(qr) (Ar) q
+#       + \sin(qr) (A)
+#    \right]
+#
+# Order 1
+#
+# .. math::
+#
+#   h(r) = \frac{4 \pi}{q^4} \left[
+#       - \cos(qr) ( Ar^2 + Br) q^2
+#       + \sin(qr) ( Ar   + B ) q
+#       + \cos(qr) (2A        )
+#   \right]
+#
+# Order 2
+#
+# .. math::
+#  h(r) = \frac{4 \pi}{q^5} \left[
+#        - \cos(qr) ( Ar^3 +  Br^2 + Cr) q^3
+#        + \sin(qr) (3Ar^2 + 2Br   + C ) q^2
+#        + \cos(qr) (6Ar   + 2B        ) q
+#        - \sin(qr) (6A                )
+#
+# Order 3
+#
+#    h(r) = \frac{4 \pi}{q^6}\left[
+#      - \cos(qr) (  Ar^4 +  Br^3 +  Cr^2 + Dr) q^4
+#      + \sin(qr) ( 4Ar^3 + 3Br^2 + 2Cr   + D ) q^3
+#      + \cos(qr) (12Ar^2 + 6Br   + 2C        ) q^2
+#      - \sin(qr) (24Ar   + 6B                ) q
+#      - \cos(qr) (24A                        )
+#    \right]
+#
+# Order p
+#
+#    h(r) = \frac{4 \pi}{q^{2}}
+#      \sum_{k=0}^p -\frac{d^k\cos(qr)}{dr^k} \frac{d^k r\rho(r)}{dr^k} (qr)^{-k}
+#
+# Given the equation
+#
+#    f = sum_(k=0)^(n-1) h_k(r_(k+1)) - h_k(r_k)
+#
+# we can rearrange the terms so that
+#
+#    f = sum_0^(n-1) h_k(r_(k+1)) - sum_0^(n-1) h_k(r_k)
+#      = sum_1^n h_(k-1)(r_k) - sum_0^(n-1) h_k(r_k)
+#      = h_(n-1)(r_n) - h_0(r_0) + sum_1^(n-1) [h_(k-1)(r_k) - h_k(r_k)]
+#      = h_(n-1)(r_n) - h_0(r_0) - sum_1^(n-1) h_(Delta k)(r_k)
+#
+# where
+#
+#    h_(Delta k)(r) = h(Delta rho_k, r)
+#
+# for
+#
+#    Delta rho_k = (A_k-A_(k-1)) r^p + (B_k-B_(k-1)) r^(p-1) + ...
+#
+# Using l'H\^opital's Rule 6 times on the order 3 polynomial,
+#
+#   lim_(q->0) h(r) = (140D r^3 + 180C r^4 + 144B r^5 + 120A r^6)/720
+#
+
+
+from __future__ import division
+
+import numpy as np
+from numpy import inf, nan
+from math import fabs, exp, expm1
+
+name = "onion"
+title = "Onion shell model with constant, linear or exponential density"
+
+description = """\
+Form factor of mutishells normalized by the volume. Here each shell is
+described by an exponential function;
+
+	I) For A_shell != 0,
+		f(r) = B*exp(A_shell*(r-r_in)/thick_shell)+C
+	where
+		B=(sld_out-sld_in)/(exp(A_shell)-1)
+		C=sld_in-B.
+	Note that in the above case, the function becomes a linear function
+	as A_shell --> 0+ or 0-.
+
+	II) For the exact point of A_shell == 0,
+		f(r) = sld_in ,i.e., it crosses over flat function
+	Note that the 'sld_out' becaomes NULL in this case.
+
+	background:background,
+	rad_core0: radius of sphere(core)
+	thick_shell#:the thickness of the shell#
+	sld_core0: the SLD of the sphere
+	sld_solv: the SLD of the solvent
+	sld_shell: the SLD of the shell#
+	A_shell#: the coefficient in the exponential function
+"""
+
+category = "shape:sphere"
+
+# TODO: n is a volume parameter that is not polydisperse
+
+#             ["name", "units", default, [lower, upper], "type","description"],
+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"],
+              ["sld_solvent", "1e-6/Ang^2", 6.4, [-inf, inf], "",
+               "Solvent scattering length density"],
+              ["n_shells", "", 1, [0, 10], "volume",
+               "number of shells"],
+              ["sld_in[n_shells]", "1e-6/Ang^2", 1.7, [-inf, inf], "",
+               "scattering length density at the inner radius of shell k"],
+              ["sld_out[n_shells]", "1e-6/Ang^2", 2.0, [-inf, inf], "",
+               "scattering length density at the outer radius of shell k"],
+              ["thickness[n_shells]", "Ang", 40., [0, inf], "volume",
+               "Thickness of shell k"],
+              ["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
+
+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):
+    """
+    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
+
+    r = []
+    beta = []
+
+    # add in the core
+    r.append(0)
+    beta.append(core_sld)
+    r.append(core_radius)
+    beta.append(core_sld)
+
+    # add in the shells
+    for k in range(n_shells):
+        # Left side of each shells
+        r0 = r[-1]
+        r.append(r0)
+        beta.append(in_sld[k])
+
+        if fabs(A[k]) < 1.0e-16:
+            # flat shell
+            r.append(r0 + thickness[k])
+            beta.append(out_sld[k])
+        else:
+            # exponential shell
+            # num_steps must be at least 1, so use floor()+1 rather than ceil
+            # to protect against a thickness0.
+            num_steps = np.floor(thickness[k]/dr) + 1
+            slope = (out_sld[k] - in_sld[k])/expm1(A[k])
+            const = (in_sld[k] - slope)
+            for rk in np.linspace(0, thickness[k], num_steps+1):
+                r.append(r0+rk)
+                beta.append(slope*exp(A[k]*rk/thickness[k]) + const)
+
+    # add in the solvent
+    r.append(r[-1])
+    beta.append(solvent_sld)
+    r.append(total_radius)
+    beta.append(solvent_sld)
+
+    return np.asarray(r), 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
+
+demo = {
+    "sld_solvent": 2.2,
+    "sld_core": 1.0,
+    "core_radius": 100,
+    "n": 4,
+    "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],
+    # Could also specify them individually as
+    # "A[1]": 0, "A[2]": -1, "A[3]": 1e-4, "A[4]": 1,
+    }
+
Index: sasmodels/models/rpa.py
===================================================================
--- sasmodels/models/rpa.py	(revision a5b8477c53695e7df03ecec7346f0c1bd770d3f7)
+++ sasmodels/models/rpa.py	(revision 46ed76021e194fbd75cd7aeac6434e865025002f)
@@ -106,4 +106,5 @@
 source = ["rpa.c"]
 
+control = "case_num"
 HIDE_NONE = set()
 HIDE_A = set("Na Phia va La Kab Kac Kad".split())
Index: sasmodels/models/spherical_sld.py
===================================================================
--- sasmodels/models/spherical_sld.py	(revision 46ed76021e194fbd75cd7aeac6434e865025002f)
+++ sasmodels/models/spherical_sld.py	(revision 46ed76021e194fbd75cd7aeac6434e865025002f)
@@ -0,0 +1,322 @@
+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
+
+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
+----------
+
+The form factor $P(q)$ in 1D is calculated by:
+
+.. math::
+
+    P(q) = \frac{f^2}{V_\text{particle}} \text{ where }
+    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}
+
+For a spherically symmetric particle with a particle density $\rho_x(r)$
+the sld function can be defined as:
+
+.. math::
+
+    f_x = 4 \pi \int_{0}^{\infty} \rho_x(r)  \frac{\sin(qr)} {qr^2} r^2 dr
+
+
+so that individual terms can be calcualted as follows:
+
+.. math::
+    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 } ) ) }
+    {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 =
+    3 \rho_\text{solvent} V(r_N)
+    \Big[ \frac{\sin(qr_N) - qr_N \cos(qr_N)} {qr_N^3} \Big]
+
+
+Here we assumed that the SLDs of the core and solvent are constant against $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
+
+Exp:
+
+.. 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 \\
+    \end{cases}
+
+Power-Law
+
+.. 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 \\
+    \rho_{\text{flat}_{i+1}}  & \text{for} A = 0 \\
+    \end{cases}
+
+Erf:
+
+.. 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 \\
+    \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 =
+    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
+
+    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}) }
+    {\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}) }
+    {\beta_\text{in}^4 } \Big]
+
+    + 3 \rho_{ \text{inter}_i } ( r_{j+1} )  V ( r_j )
+    \Big[ \frac {\sin(\beta_\text{out}) - \cos(\beta_\text{out}) }
+    {\beta_\text{out}^4 } \Big]
+
+    - 3 \rho_{ \text{inter}_i } ( r_{j} )  V ( r_j )
+    \Big[ \frac {\sin(\beta_\text{in}) - \cos(\beta_\text{in}) }
+    {\beta_\text{in}^4 } \Big]
+    \Big]
+
+where
+
+.. math::
+    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}
+
+    \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$
+
+Finally form factor can be calculated by
+
+.. math::
+
+    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,
+where the $q$ vector is defined as
+
+.. math::
+
+    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.
+
+References
+----------
+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
+
+name = "spherical_sld"
+title = "Sperical SLD intensity calculation"
+description = """
+            I(q) =
+               background = Incoherent background [1/cm]
+        """
+category = "sphere-based"
+
+# pylint: disable=bad-whitespace, line-too-long
+#            ["name", "units", default, [lower, upper], "type", "description"],
+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"]
+
+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
+    [{"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],
+]
