Index: example/model.py
===================================================================
--- example/model.py	(revision 15d2285669b798a65396200f7142e22014372cc8)
+++ example/model.py	(revision 15d2285669b798a65396200f7142e22014372cc8)
@@ -0,0 +1,40 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+import sys
+from bumps.names import *
+from sasmodels.core import load_model
+from sasmodels.bumps_model import Model, Experiment
+from sasmodels.data import load_data, set_beam_stop, set_top
+
+""" IMPORT THE DATA USED """
+radial_data = load_data('DEC07267.DAT')
+set_beam_stop(radial_data, 0.00669, outer=0.025)
+set_top(radial_data, -.0185)
+
+kernel = load_model("ellipsoid")
+
+model = Model(kernel,
+    scale=0.08,
+    rpolar=15, requatorial=800,
+    sld=.291, solvent_sld=7.105,
+    background=0,
+    theta=90, phi=0,
+    theta_pd=15, theta_pd_n=40, theta_pd_nsigma=3,
+    rpolar_pd=0.222296, rpolar_pd_n=1, rpolar_pd_nsigma=0,
+    requatorial_pd=.000128, requatorial_pd_n=1, requatorial_pd_nsigma=0,
+    phi_pd=0, phi_pd_n=20, phi_pd_nsigma=3,
+    )
+
+# SET THE FITTING PARAMETERS
+model.rpolar.range(15, 1000)
+model.requatorial.range(15, 1000)
+model.theta_pd.range(0, 360)
+model.background.range(0,1000)
+model.scale.range(0, 10)
+
+#cutoff = 0     # no cutoff on polydisperisity loops
+#cutoff = 1e-5  # default cutoff
+cutoff = 1e-3  # low precision cutoff
+M = Experiment(data=radial_data, model=model, cutoff=cutoff)
+problem = FitProblem(M)
Index: sasmodels/compare.py
===================================================================
--- sasmodels/compare.py	(revision a4a7308d68395cf588c20df822e1a84f1836d152)
+++ sasmodels/compare.py	(revision 6869cebe777d56ffcea7cfcffb97e63d951bbf84)
@@ -31,6 +31,4 @@
 import sys
 import math
-from os.path import basename, dirname, join as joinpath
-import glob
 import datetime
 import traceback
@@ -98,8 +96,5 @@
 kerneldll.ALLOW_SINGLE_PRECISION_DLLS = True
 
-# List of available models
-ROOT = dirname(__file__)
-MODELS = [basename(f)[:-3]
-          for f in sorted(glob.glob(joinpath(ROOT, "models", "[a-zA-Z]*.py")))]
+MODELS = core.list_models()
 
 # CRUFT python 2.6
@@ -671,4 +666,5 @@
     Parse command line options.
     """
+    MODELS = core.list_models()
     flags = [arg for arg in sys.argv[1:]
              if arg.startswith('-')]
@@ -683,6 +679,11 @@
         print(columnize(MODELS, indent="  "))
         sys.exit(1)
-    if args[0] not in MODELS:
-        print("Model %r not available. Use one of:\n    %s"%(args[0], models))
+
+    name = args[0]
+    try:
+        model_definition = core.load_model_definition(name)
+    except ImportError, exc:
+        print(str(exc))
+        print("Use one of:\n    " + models)
         sys.exit(1)
     if len(args) > 3:
@@ -764,7 +765,4 @@
         del engines[2:]
 
-    name = args[0]
-    model_definition = core.load_model_definition(name)
-
     n1 = int(args[1]) if len(args) > 1 else 1
     n2 = int(args[2]) if len(args) > 2 else 1
Index: sasmodels/data.py
===================================================================
--- sasmodels/data.py	(revision d18582e2d66b12afb03e1c3a7c0af6e4449c6bfe)
+++ sasmodels/data.py	(revision 78242766ff138c39a869935c18f0396d8f9a3734)
@@ -548,4 +548,5 @@
     #plottable = Iq
     plottable = masked_array(image, ~valid | data.mask)
+    # Divide range by 10 to convert from angstroms to nanometers
     xmin, xmax = min(data.qx_data)/10, max(data.qx_data)/10
     ymin, ymax = min(data.qy_data)/10, max(data.qy_data)/10
Index: sasmodels/models/_onion.py
===================================================================
--- sasmodels/models/_onion.py	(revision fdb148728a4d4bdcb4cd711d87631f3ee41993fd)
+++ sasmodels/models/_onion.py	(revision fdb148728a4d4bdcb4cd711d87631f3ee41993fd)
@@ -0,0 +1,403 @@
+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_annotated_profile.gif
+
+    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.
+
+.. figure:: img/onion_1d.jpg
+
+    1D plot using the default values (w/400 point)
+
+.. figure:: img/onion_profile.jpg
+
+    SLD profile from the default values.
+
+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"
+
+
+#             ["name", "units", default, [lower, upper], "type","description"],
+parameters = [["core_sld", "1e-6/Ang^2", 1.0, [-inf, inf], "",
+               "Core scattering length density"],
+              ["core_radius", "Ang", 200., [0, inf], "",
+               "Radius of the core"],
+              ["solvent_sld", "1e-6/Ang^2", 6.4, [-inf, inf], "",
+               "Solvent scattering length density"],
+              ["n", "", 1, [0, 10], "volume",
+               "number of shells"],
+              ["in_sld[n]", "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], "",
+               "scattering length density at the outer radius of shell k"],
+              ["thickness[n]", "Ang", 40., [0, inf], "volume",
+               "Thickness of shell k"],
+              ["A[n]", "", 1.0, [-inf, inf], "",
+               "Decay rate of shell k"],
+              ]
+
+source = ["lib/sph_j1c.c", "onion.c"]
+
+
+def profile(core_sld, core_radius, solvent_sld, n, in_sld, out_sld, thickness, A):
+    """
+    Get SLD profile
+
+    Returns *(r, rho(r))* where *r* is the radius (Ang) and *rho(r)* is the
+    SLD (1/Ang^2).
+    """
+
+    total_radius = 1.25*(sum(thickness[:n]) + 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):
+        # 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)
+
+def ER(core_radius, n, thickness):
+    return np.sum(thickness[:n], axis=0) + core_radius
+
+def VR(core_radius, n, thickness):
+    return 1.0
+
+demo = {
+    "solvent_sld": 2.2,
+    "core_sld": 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],
+    "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,
+    }
+
+oldname = "OnionExpShellModel"
+oldpars = dict(
+    core_sld="sld_core0",
+    core_radius="rad_core0",
+    solvent_sld="sld_solv",
+    n="n_shells",
+    in_sld="sld_in_shell",
+    out_sld="sld_out_shell",
+    A="A_shell",
+    thickness="thick_shell",
+    # func_shell is always 2 in the user interface, so isn't included
+    )
Index: sasmodels/models/onion.c
===================================================================
--- sasmodels/models/onion.c	(revision fdb148728a4d4bdcb4cd711d87631f3ee41993fd)
+++ sasmodels/models/onion.c	(revision fdb148728a4d4bdcb4cd711d87631f3ee41993fd)
@@ -0,0 +1,98 @@
+
+static double
+f_exp(double q, double r, double sld_in, double sld_out,
+    double thickness, double A)
+{
+  const double vol = 4.0/3.0 * M_PI * r * r * r;
+  const double qr = q * r;
+  const double alpha = A * r/thickness;
+  const double bes = sph_j1c(qr);
+  const double B = (sld_out - sld_in)/expm1(A);
+  const double C = sld_in - B;
+  double fun;
+  if (qr == 0.0) {
+    fun = 1.0;
+  } else if (fabs(A) > 0.0) {
+    const double qrsq = qr*qr;
+    const double alphasq = alpha*alpha;
+    const double sumsq = alphasq + qrsq;
+    double sinqr, cosqr;
+    SINCOS(qr, sinqr, cosqr);
+    fun = -3.0(
+            ((alphasq - qrsq)*sinqr/qr - 2.0*alpha*cosqr) / sumsq
+                - (alpha*sinqr/qr - cosqr)
+        ) / sumsq;
+  } else {
+    fun = bes;
+  }
+  return vol * (B*fun + C*bes);
+}
+
+static double
+f_linear(double q, double r, double sld, double slope)
+{
+  const double vol = 4.0/3.0 * M_PI * r * r * r;
+  const double qr = q * r;
+  const double bes = sph_j1c(qr);
+  double fun = 0.0;
+  if (qr > 0.0) {
+    const double qrsq = qr*qr;
+    double sinqr, cosqr;
+    SINCOS(qr, sinqr, cosqr);
+    // Jae-He's code seems to simplify to this
+    //     fun = 3.0 * slope * r * (2.0*qr*sinqr - (qrsq-2.0)*cosqr)/(qrsq*qrsq);
+    // Rederiving the math, we get the following instead:
+    fun = 3.0 * slope * r * (2.0*cosqr + qr*sinqr)/(qrsq*qrsq);
+  }
+  return vol * (sld*bes + fun);
+}
+
+static double
+f_constant(double q, double r, double sld)
+{
+  const double bes = sph_j1c(q * r);
+  const double vol = 4.0/3.0 * M_PI * r * r * r;
+  return sld * vol * bes;
+}
+
+static double
+form_volume(double core_radius, double n, double thickness[])
+{
+  int i;
+  double r = core_radius;
+  for (i=0; i < n; i++) {
+    r += thickness[i];
+  }
+  return 4.0/3.0 * M_PI * r * r * 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[],
+    double A[])
+{
+  int i;
+  r = core_radius;
+  f = f_constant(q, r, core_sld);
+  for (i=0; i<n; i++){
+    const double r0 = r;
+    r += thickness[i];
+    if (r == r0) {
+      // no thickness, so nothing to add
+    } else if (fabs(A[i]) < 1e-16 || sld_out[i] == sld_in[i]) {
+      f -= f_constant(q, r0, sld_in[i]);
+      f += f_constant(q, r, sld_in[i]);
+    } else if (fabs(A[i]) < 1e-4) {
+      const double slope = (sld_out[i] - sld_in[i])/thickness[i];
+      f -= f_linear(q, r0, sld_in[i], slope);
+      f += f_linear(q, r, sld_out[i], slope);
+    } else {
+      f -= f_exp(q, r0, sld_in[i], sld_out[i], thickness[i], A[i]);
+      f += f_exp(q, r, sld_in[i], sld_out[i], thickness[i], A[i]);
+    }
+  }
+  f -= f_constant(q, r, solvent_sld);
+  f2 = f * f * 1.0e-4;
+
+  return f2;
+}
Index: sasmodels/models/rpa.py
===================================================================
--- sasmodels/models/rpa.py	(revision fa8011eb0e6b527e15b594945876ee39b6a46dbd)
+++ sasmodels/models/rpa.py	(revision 5cfda00e2127f7bc13806a82cd24ec5eacbae59c)
@@ -133,7 +133,7 @@
     case_num = pars.get("case_num", parameters[0][2])
     if case_num < 2:
+        return HIDE_AB
+    elif case_num < 5:
         return HIDE_A
-    elif case_num < 5:
-        return HIDE_AB
     else:
         return HIDE_NONE
