Index: sasmodels/resolution.py
===================================================================
--- sasmodels/resolution.py	(revision 33c8d731de1592fa12dc4ad97723f4bac228a31e)
+++ sasmodels/resolution.py	(revision 7954f5cf2eda1703234398d4956457761b379065)
@@ -4,13 +4,10 @@
 This defines classes for 1D and 2D resolution calculations.
 """
+from __future__ import division
 from scipy.special import erf
-from numpy import sqrt
+from numpy import sqrt, log, log10
 import numpy as np
 
-SLIT_SMEAR_POINTS = 500
 MINIMUM_RESOLUTION = 1e-8
-# TODO: Q_EXTEND_STEPS is much bigger than necessary
-Q_EXTEND_STEPS = 30   # number of extra q points above and below
-
 
 class Resolution1D(object):
@@ -34,4 +31,5 @@
         raise NotImplementedError("Subclass does not define the apply function")
 
+
 class Perfect1D(Resolution1D):
     """
@@ -46,4 +44,5 @@
         return theory
 
+
 class Pinhole1D(Resolution1D):
     r"""
@@ -55,8 +54,7 @@
 
     *q_calc* is the list of points to calculate, or None if this should
-    be autogenerated from the *q,q_width*.
+    be estimated from the *q* and *q_width*.
     """
     def __init__(self, q, q_width, q_calc=None):
-        #TODO: maybe add min_step=np.inf
         #*min_step* is the minimum point spacing to use when computing the
         #underlying model.  It should be on the order of
@@ -70,5 +68,6 @@
         # default to Perfect1D if the pinhole geometry is not defined.
         self.q, self.q_width = q, q_width
-        self.q_calc = pinhole_extend_q(q, q_width) if q_calc is None else q_calc
+        self.q_calc = pinhole_extend_q(q, q_width) \
+            if q_calc is None else np.sort(q_calc)
         self.weight_matrix = pinhole_resolution(self.q_calc,
                 self.q, np.maximum(q_width, MINIMUM_RESOLUTION))
@@ -77,4 +76,5 @@
         return apply_resolution_matrix(self.weight_matrix, theory)
 
+
 class Slit1D(Resolution1D):
     """
@@ -86,15 +86,30 @@
 
     *qy_height* slit height
-    """
-    def __init__(self, q, qx_width, qy_width):
-        if np.isscalar(qx_width):
-            qx_width = qx_width*np.ones(len(q))
-        if np.isscalar(qy_width):
-            qy_width = qy_width*np.ones(len(q))
-        self.q, self.qx_width, self.qy_width = [
-            v.flatten() for v in (q, qx_width, qy_width)]
-        self.q_calc = slit_extend_q(q, qx_width, qy_width)
+
+    *q_calc* is the list of points to calculate, or None if this should
+    be estimated from the *q* and *q_width*.
+
+    The *weight_matrix* is computed by :func:`slit1d_resolution`
+    """
+    def __init__(self, q, width, height, q_calc=None):
+        # TODO: maybe issue warnings rather than raising errors
+        if not np.isscalar(width):
+            if np.any(np.diff(width) > 0.0):
+                raise ValueError("Slit resolution requires fixed width slits")
+            width = width[0]
+        if not np.isscalar(height):
+            if np.any(np.diff(height) > 0.0):
+                raise ValueError("Slit resolution requires fixed height slits")
+            height = height[0]
+
+        # Remember what width/height was used even though we won't need them
+        # after the weight matrix is constructed
+        self.width, self.height = width, height
+
+        self.q = q.flatten()
+        self.q_calc = slit_extend_q(q, width, height) \
+            if q_calc is None else np.sort(q_calc)
         self.weight_matrix = \
-            slit_resolution(self.q_calc, self.q, self.qx_width, self.qy_width)
+            slit_resolution(self.q_calc, self.q, width, height)
 
     def apply(self, theory):
@@ -106,9 +121,10 @@
     Apply the resolution weight matrix to the computed theory function.
     """
-    #print "apply shapes",theory.shape, self.weight_matrix.shape
+    #print "apply shapes", theory.shape, weight_matrix.shape
     Iq = np.dot(theory[None,:], weight_matrix)
     #print "result shape",Iq.shape
     return Iq.flatten()
 
+
 def pinhole_resolution(q_calc, q, q_width):
     """
@@ -119,17 +135,66 @@
     *W*, the resolution smearing can be computed using *dot(W,q)*.
 
-    *q_calc* must be increasing.
-    """
+    *q_calc* must be increasing.  *q_width* must be greater than zero.
+    """
+    # The current algorithm is a midpoint rectangle rule.  In the test case,
+    # neither trapezoid nor Simpson's rule improved the accuracy.
     edges = bin_edges(q_calc)
     edges[edges<0.0] = 0.0 # clip edges below zero
-    index = (q_width > 0.0) # treat perfect resolution differently
     G = erf( (edges[:,None] - q[None,:]) / (sqrt(2.0)*q_width)[None,:] )
     weights = G[1:] - G[:-1]
-    # w = np.sum(weights, axis=0); print "w",w.shape, w
     weights /= np.sum(weights, axis=0)[None,:]
     return weights
 
 
-def pinhole_extend_q(q, q_width):
+def slit_resolution(q_calc, q, width, height):
+    r"""
+    Build a weight matrix to compute *I_s(q)* from *I(q_calc)*, given
+    $q_v$ = *width* and $q_h$ = *height*.
+
+    *width* and *height* are scalars since current instruments use the
+    same slit settings for all measurement points.
+
+    If slit height is large relative to width, use:
+
+    .. math::
+
+        I_s(q_o) = \frac{1}{\Delta q_v}
+            \int_0^{\Delta q_v} I(\sqrt{q_o^2 + u^2} du
+
+    If slit width is large relative to height, use:
+
+    .. math::
+
+        I_s(q_o) = \frac{1}{2 \Delta q_v}
+            \int_{-\Delta q_v}^{\Delta q_v} I(u) du
+    """
+    if width == 0.0 and height == 0.0:
+        #print "condition zero"
+        return 1
+
+    q_edges = bin_edges(q_calc) # Note: requires q > 0
+    q_edges[q_edges<0.0] = 0.0 # clip edges below zero
+
+    #np.set_printoptions(linewidth=10000)
+    if width <= 100.0 * height or height == 0:
+        # The current algorithm is a midpoint rectangle rule.  In the test case,
+        # neither trapezoid nor Simpson's rule improved the accuracy.
+        #print "condition h", q_edges.shape, q.shape, q_calc.shape
+        weights = np.zeros((len(q), len(q_calc)), 'd')
+        for i, qi in enumerate(q):
+            weights[i, :] = np.diff(q_to_u(q_edges, qi))
+        weights /= width
+        weights = weights.T
+    else:
+        #print "condition w"
+        # Make q_calc into a row vector, and q into a column vector
+        q, q_calc = q[None,:], q_calc[:,None]
+        in_x = (q_calc >= q-width) * (q_calc <= q+width)
+        weights = np.diff(q_edges)[:,None] * in_x
+
+    return weights
+
+
+def pinhole_extend_q(q, q_width, nsigma=3):
     """
     Given *q* and *q_width*, find a set of sampling points *q_calc* so
@@ -137,66 +202,17 @@
     function.
     """
-    # If using min_step, you could do something like:
-    #     q_calc = np.arange(min, max, min_step)
-    #     index = np.searchsorted(q_calc, q)
-    #     q_calc[index] = q
-    # A refinement would be to assign q to the nearest neighbour.  A further
-    # refinement would be to guard multiple q points between q_calc points,
-    # either by removing duplicates or by only moving the values nearest the
-    # edges.  For really sparse samples, you will want to remove all remaining
-    # points that are not within three standard deviations of a measured q.
-    q_min, q_max = np.min(q), np.max(q)
-    extended_min, extended_max = np.min(q - 3*q_width), np.max(q + 3*q_width)
-    nbins_low, nbins_high = Q_EXTEND_STEPS, Q_EXTEND_STEPS
-    if extended_min < q_min:
-        q_low = np.linspace(extended_min, q_min, nbins_low+1)[:-1]
-        q_low = q_low[q_low>0.0]
-    else:
-        q_low = []
-    if extended_max > q_max:
-        q_high = np.linspace(q_max, extended_max, nbins_high+1)[:-1]
-    else:
-        q_high = []
-    return np.concatenate([q_low, q, q_high])
-
-
-def slit_resolution(q_calc, q, qx_width, qy_width):
-    edges = bin_edges(q_calc) # Note: requires q > 0
-    edges[edges<0.0] = 0.0 # clip edges below zero
-    qy_min, qy_max = 0.0, edges[-1]
-
-    # Make q_calc into a row vector, and q, qx_width, qy_width into columns
-    # Make weights match [ q_calc X q ]
-    weights = np.zeros((len(q),len(q_calc)), 'd')
-    q_calc = q_calc[None,:]
-    q, qx_width, qy_width, edges = [
-        v[:,None] for v in (q, qx_width, qy_width, edges)]
-
-    # Loop for width (height is analytical).
-    # Condition: height >>> width, otherwise, below is not accurate enough.
-    # Smear weight numerical iteration for width>0 when height>0.
-    # When width = 0, the numerical iteration will be skipped.
-    # The resolution calculation for the height is done by direct integration,
-    # assuming the I(q'=sqrt(q_j^2-(q+shift_w)^2)) is constant within
-    # a q' bin, [q_high, q_low].
-    # In general, this weight numerical iteration for width>0 might be a rough
-    # approximation, but it must be good enough when height >>> width.
-    E_sq = edges**2
-    y_points = SLIT_SMEAR_POINTS if np.any(qy_width>0) else 1
-    qy_step = 0 if y_points == 1 else qy_width/(y_points-1)
-    for k in range(-y_points+1,y_points):
-        qy = np.clip(q + qy_step*k, qy_min, qy_max)
-        qx_low = qy
-        qx_high = sqrt(qx_low**2 + qx_width**2)
-        in_x = (q_calc>=qx_low)*(q_calc<=qx_high)
-        qy_sq = qy**2
-        weights += (sqrt(E_sq[1:]-qy_sq) - sqrt(qy_sq - E_sq[:-1]))*in_x
-    #w = np.sum(weights, axis=1); print "w",w.shape, w
-    weights /= np.sum(weights, axis=1)[:,None]
-    return weights
-
-
-def slit_extend_q(q, qx_width, qy_width):
-    return q
+    q_min, q_max = np.min(q - nsigma*q_width), np.max(q + nsigma*q_width)
+    return geometric_extrapolation(q, q_min, q_max)
+
+
+def slit_extend_q(q, width, height):
+    """
+    Given *q*, *width* and *height*, find a set of sampling points *q_calc* so
+    that each point I(q) has sufficient support from the underlying
+    function.
+    """
+    height # keep lint happy
+    q_min, q_max = np.min(q), np.max(np.sqrt(q**2 + width**2))
+    return geometric_extrapolation(q, q_min, q_max)
 
 
@@ -218,4 +234,100 @@
     return edges
 
+
+def q_to_u(q, q0):
+    """
+    Convert *q* values to *u* values for the integral computed at *q0*.
+    """
+    # array([value])**2 - value**2 is not always zero
+    qpsq = q**2 - q0**2
+    qpsq[qpsq<0] = 0
+    return sqrt(qpsq)
+
+
+def interpolate(q, max_step):
+    """
+    Returns *q_calc* with points spaced at most max_step apart.
+    """
+    step = np.diff(q)
+    index = step>max_step
+    if np.any(index):
+        inserts = []
+        for q_i,step_i in zip(q[:-1][index],step[index]):
+            n = np.ceil(step_i/max_step)
+            inserts.extend(q_i + np.arange(1,n)*(step_i/n))
+        # Extend a couple of fringes beyond the end of the data
+        inserts.extend(q[-1] + np.arange(1,8)*max_step)
+        q_calc = np.sort(np.hstack((q,inserts)))
+    else:
+        q_calc = q
+    return q_calc
+
+
+def linear_extrapolation(q, q_min, q_max):
+    """
+    Extrapolate *q* out to [*q_min*, *q_max*] using the step size in *q* as
+    a guide.  Extrapolation below uses steps the same size as the first
+    interval.  Extrapolation above uses steps the same size as the final
+    interval.
+
+    if *q_min* is zero or less then *q[0]/10* is used instead.
+    """
+    q = np.sort(q)
+    if q_min < q[0]:
+        if q_min <= 0: q_min = q[0]/10
+        delta = q[1] - q[0]
+        q_low = np.arange(q_min, q[0], delta)
+    else:
+        q_low = []
+    if q_max > q[-1]:
+        delta = q[-1] - q[-2]
+        q_high = np.arange(q[-1]+delta, q_max, delta)
+    else:
+        q_high = []
+    return np.concatenate([q_low, q, q_high])
+
+
+def geometric_extrapolation(q, q_min, q_max):
+    r"""
+    Extrapolate *q* to [*q_min*, *q_max*] using geometric steps, with the
+    average geometric step size in *q* as the step size.
+
+    if *q_min* is zero or less then *q[0]/10* is used instead.
+
+    Starting at $q_1$ and stepping geometrically by $\Delta q$
+    to $q_n$ in $n$ points gives a geometric average of:
+
+    .. math::
+
+         \log \Delta q = (\log q_n - log q_1) / (n - 1)
+
+    From this we can compute the number of steps required to extend $q$
+    from $q_n$ to $q_\text{max}$ by $\Delta q$ as:
+
+    .. math::
+
+         n_\text{extend} = (\log q_\text{max} - \log q_n) / \log \Delta q
+
+    Substituting:
+
+        n_\text{extend} = (n-1) (\log q_\text{max} - \log q_n)
+            / (\log q_n - log q_1)
+    """
+    q = np.sort(q)
+    delta_q = (len(q)-1)/(log(q[-1]) - log(q[0]))
+    if q_min < q[0]:
+        if q_min < 0: q_min = q[0]/10
+        n_low = delta_q * (log(q[0])-log(q_min))
+        q_low  = np.logspace(log10(q_min), log10(q[0]), np.ceil(n_low)+1)[:-1]
+    else:
+        q_low = []
+    if q_max > q[-1]:
+        n_high = delta_q * (log(q_max)-log(q[-1]))
+        q_high = np.logspace(log10(q[-1]), log10(q_max), np.ceil(n_high)+1)[1:]
+    else:
+        q_high = []
+    return np.concatenate([q_low, q, q_high])
+
+
 ############################################################################
 # unit tests
@@ -223,4 +335,50 @@
 import unittest
 
+
+def eval_form(q, form, pars):
+    from sasmodels import core
+    kernel = core.make_kernel(form, [q])
+    theory = core.call_kernel(kernel, pars)
+    kernel.release()
+    return theory
+
+
+def gaussian(q, q0, dq):
+    from numpy import exp, pi
+    return exp(-0.5*((q-q0)/dq)**2)/(sqrt(2*pi)*dq)
+
+
+def romberg_slit_1d(q, delta_qv, form, pars):
+    """
+    Romberg integration for slit resolution.
+
+    This is an adaptive integration technique.  It is called with settings
+    that make it slow to evaluate but give it good accuracy.
+    """
+    from scipy.integrate import romberg
+    _fn = lambda u, q0: eval_form(sqrt(q0**2 + u**2), form, pars)
+    r = [romberg(_fn, 0, delta_qv[0], args=(qi,),
+                 divmax=100, vec_func=True, tol=0, rtol=1e-8)
+         for qi in q]
+    # r should be [float, ...], but it is [array([float]), array([float]),...]
+    return np.asarray(r).flatten()/delta_qv[0]
+
+
+def romberg_pinhole_1d(q, q_width, form, pars):
+    """
+    Romberg integration for pinhole resolution.
+
+    This is an adaptive integration technique.  It is called with settings
+    that make it slow to evaluate but give it good accuracy.
+    """
+    from scipy.integrate import romberg
+
+    _fn = lambda q, q0, dq: eval_form(q, form, pars)*gaussian(q, q0, dq)
+    r = [romberg(_fn, max(qi-5*dqi,0.01*q[0]), qi+5*dqi, args=(qi, dqi),
+                 divmax=100, vec_func=True, tol=0, rtol=1e-8)
+         for qi,dqi in zip(q,q_width)]
+    return np.asarray(r).flatten()
+
+
 class ResolutionTest(unittest.TestCase):
 
@@ -246,19 +404,55 @@
         Slit smearing with perfect resolution.
         """
-        resolution = Slit1D(self.x, 0., 0.)
+        resolution = Slit1D(self.x, width=0, height=0, q_calc=self.x)
         theory = self.Iq(resolution.q_calc)
         output = resolution.apply(theory)
         np.testing.assert_equal(output, self.y)
 
-    @unittest.skip("slit smearing is still broken")
-    def test_slit(self):
+    @unittest.skip("not yet supported")
+    def test_slit_high(self):
         """
         Slit smearing with height 0.005
         """
-        resolution = Slit1D(self.x, 0., 0.005)
+        resolution = Slit1D(self.x, width=0, height=0.005, q_calc=self.x)
         theory = self.Iq(resolution.q_calc)
         output = resolution.apply(theory)
-        answer = [ 9.0618, 8.6401, 8.1186, 7.1391, 6.1528,
-                   5.5555, 4.5584, 3.5606, 2.5623, 2. ]
+        answer = [ 9.0618, 8.6402, 8.1187, 7.1392, 6.1528,
+                   5.5555, 4.5584, 3.5606, 2.5623, 2.0000 ]
+        np.testing.assert_allclose(output, answer, atol=1e-4)
+
+    @unittest.skip("not yet supported")
+    def test_slit_both_high(self):
+        """
+        Slit smearing with width < 100*height.
+        """
+        q = np.logspace(-4, -1, 10)
+        resolution = Slit1D(q, width=0.2, height=np.inf)
+        theory = 1000*self.Iq(resolution.q_calc**4)
+        output = resolution.apply(theory)
+        answer = [ 8.85785, 8.43012, 7.92687, 6.94566, 6.03660,
+                   5.40363, 4.40655, 3.40880, 2.41058, 2.00000 ]
+        np.testing.assert_allclose(output, answer, atol=1e-4)
+
+    @unittest.skip("not yet supported")
+    def test_slit_wide(self):
+        """
+        Slit smearing with width 0.0002
+        """
+        resolution = Slit1D(self.x, width=0.0002, height=0, q_calc=self.x)
+        theory = self.Iq(resolution.q_calc)
+        output = resolution.apply(theory)
+        answer = [ 11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0 ]
+        np.testing.assert_allclose(output, answer, atol=1e-4)
+
+    @unittest.skip("not yet supported")
+    def test_slit_both_wide(self):
+        """
+        Slit smearing with width > 100*height.
+        """
+        resolution = Slit1D(self.x, width=0.0002, height=0.000001,
+                            q_calc=self.x)
+        theory = self.Iq(resolution.q_calc)
+        output = resolution.apply(theory)
+        answer = [ 11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0 ]
         np.testing.assert_allclose(output, answer, atol=1e-4)
 
@@ -285,8 +479,143 @@
         np.testing.assert_allclose(output, answer, atol=1e-8)
 
-# Q, dQ, I(Q) for Igor default sphere model.
-# combines CMSphere5.txt and CMSphere5smearsphere.txt from sasview/tests
-# TODO: move test data into its own file?
-SPHERE_RESOLUTION_TEST_DATA = """\
+
+class IgorComparisonTest(unittest.TestCase):
+
+    def setUp(self):
+        self.pars = TEST_PARS_PINHOLE_SPHERE
+        from sasmodels import core
+        from sasmodels.models import sphere
+        self.model = core.load_model(sphere, dtype='double')
+
+    def Iq_sphere(self, pars, resolution):
+        from sasmodels import core
+        kernel = core.make_kernel(self.model, [resolution.q_calc])
+        theory = core.call_kernel(kernel, pars)
+        result = resolution.apply(theory)
+        kernel.release()
+        return result
+
+    def compare(self, q, output, answer, tolerance):
+        err = (output - answer)/answer
+        idx = abs(err) >= tolerance
+        problem = zip(q[idx], output[idx], answer[idx], err[idx])
+        print "\n".join(str(v) for v in problem)
+        np.testing.assert_allclose(output, answer, rtol=tolerance)
+
+    def test_perfect(self):
+        """
+        Compare sphere model with NIST Igor SANS, no resolution smearing.
+        """
+        pars = TEST_PARS_SLIT_SPHERE
+        data_string = TEST_DATA_SLIT_SPHERE
+
+        data = np.loadtxt(data_string.split('\n')).T
+        q, width, answer, _ = data
+        resolution = Perfect1D(q)
+        output = self.Iq_sphere(pars, resolution)
+        self.compare(q, output, answer, 1e-6)
+
+    def test_pinhole(self):
+        """
+        Compare pinhole resolution smearing with NIST Igor SANS
+        """
+        pars = TEST_PARS_PINHOLE_SPHERE
+        data_string = TEST_DATA_PINHOLE_SPHERE
+
+        data = np.loadtxt(data_string.split('\n')).T
+        q, q_width, answer = data
+        resolution = Pinhole1D(q, q_width)
+        output = self.Iq_sphere(pars, resolution)
+        # TODO: relative error should be lower
+        self.compare(q, output, answer, 3e-4)
+
+    def test_pinhole_romberg(self):
+        """
+        Compare pinhole resolution smearing with romberg integration result.
+        """
+        pars = TEST_PARS_PINHOLE_SPHERE
+        data_string = TEST_DATA_PINHOLE_SPHERE
+        pars['radius'] *= 5
+        radius = pars['radius']
+
+        data = np.loadtxt(data_string.split('\n')).T
+        q, q_width, answer = data
+        answer = romberg_pinhole_1d(q, q_width, self.model, pars)
+        ## Getting 0.1% requires 5 sigma and 200 points per fringe
+        #q_calc = interpolate(pinhole_extend_q(q, q_width, nsigma=5),
+        #                     2*np.pi/radius/200)
+        #tol = 0.001
+        ## The default 3 sigma and no extra points gets 1%
+        q_calc, tol = None, 0.01
+        resolution = Pinhole1D(q, q_width, q_calc=q_calc)
+        output = self.Iq_sphere(pars, resolution)
+        if 0: # debug plot
+            import matplotlib.pyplot as plt
+            resolution = Perfect1D(q)
+            source = self.Iq_sphere(pars, resolution)
+            plt.loglog(q, source, '.')
+            plt.loglog(q, answer, '-', hold=True)
+            plt.loglog(q, output, '-', hold=True)
+            plt.show()
+        self.compare(q, output, answer, tol)
+
+    def test_slit(self):
+        """
+        Compare slit resolution smearing with NIST Igor SANS
+        """
+        pars = TEST_PARS_SLIT_SPHERE
+        data_string = TEST_DATA_SLIT_SPHERE
+
+        data = np.loadtxt(data_string.split('\n')).T
+        q, delta_qv, _, answer = data
+        resolution = Slit1D(q, width=delta_qv, height=0)
+        output = self.Iq_sphere(pars, resolution)
+        # TODO: eliminate Igor test since it is too inaccurate to be useful.
+        # This means we can eliminate the test data as well, and instead
+        # use a generated q vector.
+        self.compare(q, output, answer, 0.5)
+
+    def test_slit_romberg(self):
+        """
+        Compare slit resolution smearing with romberg integration result.
+        """
+        pars = TEST_PARS_SLIT_SPHERE
+        data_string = TEST_DATA_SLIT_SPHERE
+        radius = pars['radius']
+
+        data = np.loadtxt(data_string.split('\n')).T
+        q, delta_qv, _, answer = data
+        answer = romberg_slit_1d(q, delta_qv, self.model, pars)
+        q_calc = slit_extend_q(interpolate(q, 2*np.pi/radius/20),
+                               delta_qv[0], delta_qv[0])
+        resolution = Slit1D(q, width=delta_qv, height=0, q_calc=q_calc)
+        output = self.Iq_sphere(pars, resolution)
+        # TODO: relative error should be lower
+        self.compare(q, output, answer, 0.025)
+
+    #TODO: can sas q spacing be too sparse for the resolution calculation?
+    @unittest.skip("suppress sparse data test; not supported by current code")
+    def test_pinhole_sparse(self):
+        """
+        Compare pinhole resolution smearing with NIST Igor SANS on sparse data
+        """
+        pars = TEST_PARS_PINHOLE_SPHERE
+        data_string = TEST_DATA_PINHOLE_SPHERE
+
+        data = np.loadtxt(data_string.split('\n')).T
+        q, q_width, answer = data[:, ::20] # Take every nth point
+        resolution = Pinhole1D(q, q_width)
+        output = self.Iq_sphere(pars, resolution)
+        self.compare(q, output, answer, 1e-6)
+
+
+
+# pinhole sphere parameters
+TEST_PARS_PINHOLE_SPHERE = {
+    'scale': 1.0, 'background': 0.01,
+    'radius': 60.0, 'sld': 1, 'solvent_sld': 6.3,
+    }
+# Q, dQ, I(Q) calculated by NIST Igor SANS package
+TEST_DATA_PINHOLE_SPHERE = """\
 0.001278 0.0002847 2538.41176383
 0.001562 0.0002905 2536.91820405
@@ -391,41 +720,80 @@
 """
 
-class Pinhole1DTest(unittest.TestCase):
-
-    def setUp(self):
-        # sample q, q_width, I(q) calculated by NIST Igor SANS package
-        self.data = np.loadtxt(SPHERE_RESOLUTION_TEST_DATA.split('\n')).T
-        self.pars = dict(scale=1.0, background=0.01, radius=60.0,
-                         solvent_sld=6.3, sld=1)
-
-    def Iq(self, q, q_width):
-        from sasmodels import core
-        from sasmodels.models import sphere
-
-        model = core.load_model(sphere)
-        resolution = Pinhole1D(q, q_width)
-        kernel = core.make_kernel(model, [resolution.q_calc])
-        theory = core.call_kernel(kernel, self.pars)
-        result = resolution.apply(theory)
-        return result
-
-    def test_sphere(self):
-        """
-        Compare pinhole resolution smearing with NIST Igor SANS
-        """
-        q, q_width, answer = self.data
-        output = self.Iq(q, q_width)
-        np.testing.assert_allclose(output, answer, rtol=0.006)
-
-    #TODO: is all sas data sampled densely enough to support resolution calcs?
-    @unittest.skip("suppress sparse data test; not supported by current code")
-    def test_sphere_sparse(self):
-        """
-        Compare pinhole resolution smearing with NIST Igor SANS on sparse data
-        """
-        q, q_width, answer = self.data[:, ::20]  # Take every nth point
-        output = self.Iq(q, q_width)
-        np.testing.assert_allclose(output, answer, rtol=0.006)
-
+# Slit sphere parameters
+TEST_PARS_SLIT_SPHERE = {
+    'scale': 0.01, 'background': 0.01,
+    'radius': 60000, 'sld': 1, 'solvent_sld': 4,
+    }
+# Q dQ I(Q) I_smeared(Q)
+TEST_DATA_SLIT_SPHERE = """\
+2.26097e-05 0.117 5.5781372896e+09 1.4626077708e+06
+2.53847e-05 0.117 5.0363141458e+09 1.3117318023e+06
+2.81597e-05 0.117 4.4850108103e+09 1.1594863713e+06
+3.09347e-05 0.117 3.9364658459e+09 1.0094881630e+06
+3.37097e-05 0.117 3.4019975074e+09 8.6518941303e+05
+3.92597e-05 0.117 2.4139519814e+09 6.0232158311e+05
+4.48097e-05 0.117 1.5816877820e+09 3.8739994090e+05
+5.03597e-05 0.117 9.3715407224e+08 2.2745304775e+05
+5.59097e-05 0.117 4.8387917428e+08 1.2101295768e+05
+6.14597e-05 0.117 2.0193586928e+08 6.0055107771e+04
+6.70097e-05 0.117 5.5886110911e+07 3.2749521065e+04
+7.25597e-05 0.117 3.7782348010e+06 2.6350963616e+04
+7.81097e-05 0.117 5.3407817904e+06 2.9624963314e+04
+8.36597e-05 0.117 2.7975485523e+07 3.4403962254e+04
+8.92097e-05 0.117 4.9845448282e+07 3.6130017903e+04
+9.47597e-05 0.117 6.0092588905e+07 3.3495107849e+04
+1.00310e-04 0.117 5.6823430831e+07 2.7475726279e+04
+1.05860e-04 0.117 4.3857024036e+07 2.0144282226e+04
+1.11410e-04 0.117 2.7277144760e+07 1.3647403260e+04
+1.22510e-04 0.117 3.3119334113e+06 6.6519711526e+03
+1.33610e-04 0.117 1.4412859402e+06 6.9726212813e+03
+1.44710e-04 0.117 8.5620162463e+06 8.1441335775e+03
+1.55810e-04 0.117 9.6957429033e+06 6.4559996521e+03
+1.66910e-04 0.117 4.3818341914e+06 3.6252154156e+03
+1.78010e-04 0.117 2.7448997387e+05 2.4006505342e+03
+1.89110e-04 0.117 8.0472009936e+05 2.8187789089e+03
+2.00210e-04 0.117 2.8149552834e+06 3.0915662855e+03
+2.11310e-04 0.117 2.7510907861e+06 2.3722530293e+03
+2.22410e-04 0.117 1.0053133293e+06 1.4473468311e+03
+2.33510e-04 0.117 5.8428305052e+03 1.2048540556e+03
+2.44610e-04 0.117 5.1699305004e+05 1.4625670042e+03
+2.55710e-04 0.117 1.2120227268e+06 1.5010705973e+03
+2.66810e-04 0.117 9.7896842846e+05 1.1336343426e+03
+2.77910e-04 0.117 2.5507264791e+05 8.1848818080e+02
+3.05660e-04 0.117 5.2403101181e+05 7.4913374357e+02
+3.33410e-04 0.117 5.8699343809e+04 4.4669964560e+02
+3.61160e-04 0.117 3.0844327150e+05 4.6774007542e+02
+3.88910e-04 0.117 8.3360142970e+03 2.7169550220e+02
+4.16660e-04 0.117 1.8630080583e+05 3.0710983679e+02
+4.44410e-04 0.117 3.1616804732e-01 1.7959006831e+02
+4.72160e-04 0.117 1.1299016314e+05 2.0763952339e+02
+4.99910e-04 0.117 2.9952522747e+03 1.2536542765e+02
+5.27660e-04 0.117 6.7625695649e+04 1.4013969777e+02
+5.55410e-04 0.117 7.6927460089e+03 8.2145593180e+01
+6.10910e-04 0.117 1.1229057779e+04 8.4519745643e+01
+6.66410e-04 0.117 1.3035567943e+04 8.1554625609e+01
+7.21910e-04 0.117 1.3309931343e+04 7.4437319172e+01
+7.77410e-04 0.117 1.2462626212e+04 6.4697088261e+01
+8.32910e-04 0.117 1.0912927143e+04 5.3773301044e+01
+8.88410e-04 0.117 9.0172597469e+03 4.2843375753e+01
+9.43910e-04 0.117 7.0496495917e+03 3.2771032724e+01
+9.99410e-04 0.117 5.2030483682e+03 2.4113557144e+01
+1.05491e-03 0.117 3.5988976711e+03 1.7160773658e+01
+1.11041e-03 0.117 2.2996060652e+03 1.2016626459e+01
+1.22141e-03 0.117 6.4766590598e+02 6.0373017740e+00
+1.33241e-03 0.117 4.1963483264e+01 4.5215452974e+00
+1.44341e-03 0.117 6.3370708246e+01 5.1054681903e+00
+1.55441e-03 0.117 3.0736750577e+02 5.9176165298e+00
+1.66541e-03 0.117 5.0327682399e+02 5.9815000189e+00
+1.77641e-03 0.117 5.4084331454e+02 5.1634639625e+00
+1.88741e-03 0.117 4.3488671756e+02 3.8535158148e+00
+1.99841e-03 0.117 2.6322287860e+02 2.5824997753e+00
+2.10941e-03 0.117 1.0793633150e+02 1.7315517194e+00
+2.22041e-03 0.117 1.8474448850e+01 1.4077213604e+00
+2.33141e-03 0.117 1.5864062279e+00 1.4771560682e+00
+2.44241e-03 0.117 3.2267213848e+01 1.6916253448e+00
+2.55341e-03 0.117 7.4289116207e+01 1.8274751193e+00
+2.66441e-03 0.117 9.9000521929e+01 1.7706812289e+00
+"""
 
 def main():
@@ -494,4 +862,2 @@
     #demo()
     main()
-
-
