-                      r33c8d73
+                      r7954f5c
 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):
 …
         raise NotImplementedError("Subclass does not define the apply function")
 class Perfect1D(Resolution1D):
     """
 …
         return theory
 class Pinhole1D(Resolution1D):
     r"""
 …
     *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
 …
         # 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))
 …
         return apply_resolution_matrix(self.weight_matrix, theory)
 class Slit1D(Resolution1D):
     """
 …
     *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):
 …
     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):
     """
 …
     *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
 …
     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)
 …
     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
 …
 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):
 …
         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,
+.5555, 4.5584, 3.5606, 2.5623, 2. ]
+        answer = [ 9.0618, 8.6402, 8.1187, 7.1392, 6.1528,
+.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,
+.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)
 …
         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 = """\
 .001278 0.0002847 2538.41176383
 .001562 0.0002905 2536.91820405
 …
 """
+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 = """\
+.26097e-05 0.117 5.5781372896e+09 1.4626077708e+06
+.53847e-05 0.117 5.0363141458e+09 1.3117318023e+06
+.81597e-05 0.117 4.4850108103e+09 1.1594863713e+06
+.09347e-05 0.117 3.9364658459e+09 1.0094881630e+06
+.37097e-05 0.117 3.4019975074e+09 8.6518941303e+05
+.92597e-05 0.117 2.4139519814e+09 6.0232158311e+05
+.48097e-05 0.117 1.5816877820e+09 3.8739994090e+05
+.03597e-05 0.117 9.3715407224e+08 2.2745304775e+05
+.59097e-05 0.117 4.8387917428e+08 1.2101295768e+05
+.14597e-05 0.117 2.0193586928e+08 6.0055107771e+04
+.70097e-05 0.117 5.5886110911e+07 3.2749521065e+04
+.25597e-05 0.117 3.7782348010e+06 2.6350963616e+04
+.81097e-05 0.117 5.3407817904e+06 2.9624963314e+04
+.36597e-05 0.117 2.7975485523e+07 3.4403962254e+04
+.92097e-05 0.117 4.9845448282e+07 3.6130017903e+04
+.47597e-05 0.117 6.0092588905e+07 3.3495107849e+04
+.00310e-04 0.117 5.6823430831e+07 2.7475726279e+04
+.05860e-04 0.117 4.3857024036e+07 2.0144282226e+04
+.11410e-04 0.117 2.7277144760e+07 1.3647403260e+04
+.22510e-04 0.117 3.3119334113e+06 6.6519711526e+03
+.33610e-04 0.117 1.4412859402e+06 6.9726212813e+03
+.44710e-04 0.117 8.5620162463e+06 8.1441335775e+03
+.55810e-04 0.117 9.6957429033e+06 6.4559996521e+03
+.66910e-04 0.117 4.3818341914e+06 3.6252154156e+03
+.78010e-04 0.117 2.7448997387e+05 2.4006505342e+03
+.89110e-04 0.117 8.0472009936e+05 2.8187789089e+03
+.00210e-04 0.117 2.8149552834e+06 3.0915662855e+03
+.11310e-04 0.117 2.7510907861e+06 2.3722530293e+03
+.22410e-04 0.117 1.0053133293e+06 1.4473468311e+03
+.33510e-04 0.117 5.8428305052e+03 1.2048540556e+03
+.44610e-04 0.117 5.1699305004e+05 1.4625670042e+03
+.55710e-04 0.117 1.2120227268e+06 1.5010705973e+03
+.66810e-04 0.117 9.7896842846e+05 1.1336343426e+03
+.77910e-04 0.117 2.5507264791e+05 8.1848818080e+02
+.05660e-04 0.117 5.2403101181e+05 7.4913374357e+02
+.33410e-04 0.117 5.8699343809e+04 4.4669964560e+02
+.61160e-04 0.117 3.0844327150e+05 4.6774007542e+02
+.88910e-04 0.117 8.3360142970e+03 2.7169550220e+02
+.16660e-04 0.117 1.8630080583e+05 3.0710983679e+02
+.44410e-04 0.117 3.1616804732e-01 1.7959006831e+02
+.72160e-04 0.117 1.1299016314e+05 2.0763952339e+02
+.99910e-04 0.117 2.9952522747e+03 1.2536542765e+02
+.27660e-04 0.117 6.7625695649e+04 1.4013969777e+02
+.55410e-04 0.117 7.6927460089e+03 8.2145593180e+01
+.10910e-04 0.117 1.1229057779e+04 8.4519745643e+01
+.66410e-04 0.117 1.3035567943e+04 8.1554625609e+01
+.21910e-04 0.117 1.3309931343e+04 7.4437319172e+01
+.77410e-04 0.117 1.2462626212e+04 6.4697088261e+01
+.32910e-04 0.117 1.0912927143e+04 5.3773301044e+01
+.88410e-04 0.117 9.0172597469e+03 4.2843375753e+01
+.43910e-04 0.117 7.0496495917e+03 3.2771032724e+01
+.99410e-04 0.117 5.2030483682e+03 2.4113557144e+01
+.05491e-03 0.117 3.5988976711e+03 1.7160773658e+01
+.11041e-03 0.117 2.2996060652e+03 1.2016626459e+01
+.22141e-03 0.117 6.4766590598e+02 6.0373017740e+00
+.33241e-03 0.117 4.1963483264e+01 4.5215452974e+00
+.44341e-03 0.117 6.3370708246e+01 5.1054681903e+00
+.55441e-03 0.117 3.0736750577e+02 5.9176165298e+00
+.66541e-03 0.117 5.0327682399e+02 5.9815000189e+00
+.77641e-03 0.117 5.4084331454e+02 5.1634639625e+00
+.88741e-03 0.117 4.3488671756e+02 3.8535158148e+00
+.99841e-03 0.117 2.6322287860e+02 2.5824997753e+00
+.10941e-03 0.117 1.0793633150e+02 1.7315517194e+00
+.22041e-03 0.117 1.8474448850e+01 1.4077213604e+00
+.33141e-03 0.117 1.5864062279e+00 1.4771560682e+00
+.44241e-03 0.117 3.2267213848e+01 1.6916253448e+00
+.55341e-03 0.117 7.4289116207e+01 1.8274751193e+00
+.66441e-03 0.117 9.9000521929e+01 1.7706812289e+00
+"""
 def main():
 …
     #demo()
     main()

Note: See TracChangeset for help on using the changeset viewer.

SasView

Changeset 7954f5c in sasmodels

Legend:

sasmodels/resolution.py

Download in other formats: