Index: compare.py
===================================================================
--- compare.py	(revision aa4946b68cdcd8686f2c8386cafe8a9e446f0ed2)
+++ compare.py	(revision 346bc886df88ad4f5c0e3441573497ef436dbf18)
@@ -9,5 +9,5 @@
 import numpy as np
 
-from sasmodels.bumps_model import BumpsModel, plot_data, tic
+from sasmodels.bumps_model import Model, Experiment, plot_theory, tic
 from sasmodels import core
 from sasmodels import kerneldll
@@ -97,11 +97,23 @@
 
 def eval_sasview(name, pars, data, Nevals=1):
+    from sas.models.qsmearing import smear_selection
     model = sasview_model(name, **pars)
+    smearer = smear_selection(data, model=model)
     toc = tic()
     for _ in range(Nevals):
         if hasattr(data, 'qx_data'):
-            value = model.evalDistribution([data.qx_data, data.qy_data])
+            q = np.sqrt(data.qx_data**2 + data.qy_data**2)
+            index = ((~data.mask) & (~np.isnan(data.data))
+                     & (q >= data.qmin) & (q <= data.qmax))
+            if smearer is not None:
+                smearer.model = model  # because smear_selection has a bug
+                smearer.set_index(index)
+                value = smearer.get_value()
+            else:
+                value = model.evalDistribution([data.qx_data[index], data.qy_data[index]])
         else:
             value = model.evalDistribution(data.x)
+            if smearer is not None:
+                value = smearer(value)
     average_time = toc()*1000./Nevals
     return value, average_time
@@ -114,5 +126,5 @@
         print "... trying again with single precision"
         model = core.load_model(model_definition, dtype='single', platform="ocl")
-    problem = BumpsModel(data, model, cutoff=cutoff, **pars)
+    problem = Experiment(data, Model(model, **pars), cutoff=cutoff)
     toc = tic()
     for _ in range(Nevals):
@@ -125,5 +137,5 @@
 def eval_ctypes(model_definition, pars, data, dtype='double', Nevals=1, cutoff=0):
     model = core.load_model(model_definition, dtype=dtype, platform="dll")
-    problem = BumpsModel(data, model, cutoff=cutoff, **pars)
+    problem = Experiment(data, Model(model, **pars), cutoff=cutoff)
     toc = tic()
     for _ in range(Nevals):
@@ -133,8 +145,8 @@
     return value, average_time
 
-def make_data(qmax, is2D, Nq=128, view='log'):
+def make_data(qmax, is2D, Nq=128, resolution=0.0, view='log'):
     if is2D:
         from sasmodels.bumps_model import empty_data2D, set_beam_stop
-        data = empty_data2D(np.linspace(-qmax, qmax, Nq))
+        data = empty_data2D(np.linspace(-qmax, qmax, Nq), resolution=resolution)
         set_beam_stop(data, 0.004)
         index = ~data.mask
@@ -146,5 +158,5 @@
         else:
             q = np.linspace(0.001*qmax, qmax, Nq)
-        data = empty_data1D(q)
+        data = empty_data1D(q, resolution=resolution)
         index = slice(None, None)
     return data, index
@@ -159,6 +171,7 @@
     qmax = 10.0 if '-exq' in opts else 1.0 if '-highq' in opts else 0.2 if '-midq' in opts else 0.05
     Nq = int(opt_values.get('-Nq', '128'))
+    res = float(opt_values.get('-res', '0'))
     is2D = not "-1d" in opts
-    data, index = make_data(qmax, is2D, Nq, view=view)
+    data, index = make_data(qmax, is2D, Nq, res, view=view)
 
 
@@ -185,5 +198,5 @@
         ocl, ocl_time = eval_opencl(model_definition, pars, data,
                                     dtype=dtype, cutoff=cutoff, Nevals=Nocl)
-        print "opencl t=%.1f ms, intensity=%.0f"%(ocl_time, sum(ocl[index]))
+        print "opencl t=%.1f ms, intensity=%.0f"%(ocl_time, sum(ocl))
         #print max(ocl), min(ocl)
 
@@ -193,9 +206,9 @@
                                     dtype=dtype, cutoff=cutoff, Nevals=Ncpu)
         comp = "ctypes"
-        print "ctypes t=%.1f ms, intensity=%.0f"%(cpu_time, sum(cpu[index]))
+        print "ctypes t=%.1f ms, intensity=%.0f"%(cpu_time, sum(cpu))
     elif Ncpu > 0:
         cpu, cpu_time = eval_sasview(model_definition, pars, data, Ncpu)
         comp = "sasview"
-        print "sasview t=%.1f ms, intensity=%.0f"%(cpu_time, sum(cpu[index]))
+        print "sasview t=%.1f ms, intensity=%.0f"%(cpu_time, sum(cpu))
 
     # Compare, but only if computing both forms
@@ -204,7 +217,6 @@
         #print "max |ocl/cpu|", max(abs(ocl/cpu)), "%.15g"%max(abs(ocl)), "%.15g"%max(abs(cpu))
         #cpu *= max(ocl/cpu)
-        resid, relerr = np.zeros_like(ocl), np.zeros_like(ocl)
-        resid[index] = (ocl - cpu)[index]
-        relerr[index] = resid[index]/cpu[index]
+        resid = (ocl - cpu)
+        relerr = resid/cpu
         #bad = (relerr>1e-4)
         #print relerr[bad],cpu[bad],ocl[bad],data.qx_data[bad],data.qy_data[bad]
@@ -221,6 +233,6 @@
                 ]
             print label,"  ".join(data)
-        stats("|ocl-%s|"%comp+(" "*(3+len(comp))), resid[index])
-        stats("|(ocl-%s)/%s|"%(comp,comp), relerr[index])
+        stats("|ocl-%s|"%comp+(" "*(3+len(comp))), resid)
+        stats("|(ocl-%s)/%s|"%(comp,comp), relerr)
 
     # Plot if requested
@@ -229,10 +241,10 @@
     if Ncpu > 0:
         if Nocl > 0: plt.subplot(131)
-        plot_data(data, cpu, view=view)
+        plot_theory(data, cpu, view=view)
         plt.title("%s t=%.1f ms"%(comp,cpu_time))
         cbar_title = "log I"
     if Nocl > 0:
         if Ncpu > 0: plt.subplot(132)
-        plot_data(data, ocl, view=view)
+        plot_theory(data, ocl, view=view)
         plt.title("opencl t=%.1f ms"%ocl_time)
         cbar_title = "log I"
@@ -244,6 +256,6 @@
             err,errstr,errview = abs(relerr), "rel err", "log"
         #err,errstr = ocl/cpu,"ratio"
-        plot_data(data, err, view=errview)
-        plt.title("max %s = %.3g"%(errstr, max(abs(err[index]))))
+        plot_theory(data, err, view=errview)
+        plt.title("max %s = %.3g"%(errstr, max(abs(err))))
         cbar_title = errstr if errview=="linear" else "log "+errstr
     if is2D:
@@ -253,5 +265,5 @@
     if Ncpu > 0 and Nocl > 0 and '-hist' in opts:
         plt.figure()
-        v = relerr[index]
+        v = relerr
         v[v==0] = 0.5*np.min(np.abs(v[v!=0]))
         plt.hist(np.log10(np.abs(v)), normed=1, bins=50);
@@ -289,4 +301,5 @@
     -linear/-log/-q4 intensity scaling
     -hist/-nohist* plot histogram of relative error
+    -res=0 sets the resolution width dQ/Q if calculating with resolution
 
 Key=value pairs allow you to set specific values to any of the model
@@ -313,5 +326,5 @@
 VALUE_OPTIONS = [
     # Note: random is both a name option and a value option
-    'cutoff', 'random', 'Nq',
+    'cutoff', 'random', 'Nq', 'res',
     ]
 
Index: example/fit.py
===================================================================
--- example/fit.py	(revision 5d3d7b49ed90178c7d8dbd9d26bdcd19db5cf66b)
+++ example/fit.py	(revision 346bc886df88ad4f5c0e3441573497ef436dbf18)
@@ -4,4 +4,5 @@
 import sys
 from bumps.names import *
+from sasmodels.core import load_model
 from sasmodels import bumps_model as sas
 
@@ -23,37 +24,39 @@
 data = radial_data if section != "tangential" else tan_data
 phi = 0 if section != "tangential" else 90
-kernel = sas.load_model(name, dtype="single")
+kernel = load_model(name, dtype="single")
 cutoff = 1e-3
 
 if name == "ellipsoid":
-    model = sas.BumpsModel(data, kernel,
-    scale=0.08,
-    rpolar=15, requatorial=800,
-    sld=.291, solvent_sld=7.105,
-    background=0,
-    theta=90, phi=phi,
-    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,
-    )
+    model = sas.Model(kernel,
+        scale=0.08,
+        rpolar=15, requatorial=800,
+        sld=.291, solvent_sld=7.105,
+        background=0,
+        theta=90, phi=phi,
+        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.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)
 
 
+
 elif name == "lamellar":
-    model = sas.BumpsModel(data, kernel,
-    scale=0.08,
-    thickness=19.2946,
-    sld=5.38,sld_sol=7.105,
-    background=0.003,
-    thickness_pd= 0.37765, thickness_pd_n=10, thickness_pd_nsigma=3,
-    )
+    model = sas.Model(kernel,
+        scale=0.08,
+        thickness=19.2946,
+        sld=5.38,sld_sol=7.105,
+        background=0.003,
+        thickness_pd= 0.37765, thickness_pd_n=10, thickness_pd_nsigma=3,
+        )
+
 
     # SET THE FITTING PARAMETERS
@@ -84,5 +87,5 @@
         theta_pd=10, theta_pd_n=50, theta_pd_nsigma=3,
         phi_pd=0, phi_pd_n=10, phi_pd_nsigma=3)
-    model = sas.BumpsModel(data, kernel, **pars)
+    model = sas.Model(kernel, **pars)
 
     # SET THE FITTING PARAMETERS
@@ -99,15 +102,15 @@
 
 elif name == "core_shell_cylinder":
-    model = sas.BumpsModel(data, kernel,
-    scale= .031, radius=19.5, thickness=30, length=22,
-    core_sld=7.105, shell_sld=.291, solvent_sld=7.105,
-    background=0, theta=0, phi=phi,
+    model = sas.Model(kernel,
+        scale= .031, radius=19.5, thickness=30, length=22,
+        core_sld=7.105, shell_sld=.291, solvent_sld=7.105,
+        background=0, theta=0, phi=phi,
 
-    radius_pd=0.26, radius_pd_n=10, radius_pd_nsigma=3,
-    length_pd=0.26, length_pd_n=10, length_pd_nsigma=3,
-    thickness_pd=1, thickness_pd_n=1, thickness_pd_nsigma=1,
-    theta_pd=1, theta_pd_n=10, theta_pd_nsigma=3,
-    phi_pd=0.1, phi_pd_n=1, phi_pd_nsigma=1,
-    )
+        radius_pd=0.26, radius_pd_n=10, radius_pd_nsigma=3,
+        length_pd=0.26, length_pd_n=10, length_pd_nsigma=3,
+        thickness_pd=1, thickness_pd_n=1, thickness_pd_nsigma=1,
+        theta_pd=1, theta_pd_n=10, theta_pd_nsigma=3,
+        phi_pd=0.1, phi_pd_n=1, phi_pd_nsigma=1,
+        )
 
     # SET THE FITTING PARAMETERS
@@ -126,14 +129,14 @@
 
 elif name == "capped_cylinder":
-    model = sas.BumpsModel(data, kernel,
-    scale=.08, radius=20, cap_radius=40, length=400,
-    sld_capcyl=1, sld_solv=6.3,
-    background=0, theta=0, phi=phi,
-    radius_pd=.1, radius_pd_n=5, radius_pd_nsigma=3,
-    cap_radius_pd=.1, cap_radius_pd_n=5, cap_radius_pd_nsigma=3,
-    length_pd=.1, length_pd_n=1, length_pd_nsigma=0,
-    theta_pd=.1, theta_pd_n=1, theta_pd_nsigma=0,
-    phi_pd=.1, phi_pd_n=1, phi_pd_nsigma=0,
-    dtype=dtype)
+    model = sas.Model(kernel,
+        scale=.08, radius=20, cap_radius=40, length=400,
+        sld_capcyl=1, sld_solv=6.3,
+        background=0, theta=0, phi=phi,
+        radius_pd=.1, radius_pd_n=5, radius_pd_nsigma=3,
+        cap_radius_pd=.1, cap_radius_pd_n=5, cap_radius_pd_nsigma=3,
+        length_pd=.1, length_pd_n=1, length_pd_nsigma=0,
+        theta_pd=.1, theta_pd_n=1, theta_pd_nsigma=0,
+        phi_pd=.1, phi_pd_n=1, phi_pd_nsigma=0,
+        )
 
     model.scale.range(0, 1)
@@ -141,14 +144,15 @@
 
 elif name == "triaxial_ellipsoid":
-    model = sas.BumpsModel(data, kernel,
-    scale=0.08, req_minor=15, req_major=20, rpolar=500,
-    sldEll=7.105, solvent_sld=.291,
-    background=5, theta=0, phi=phi, psi=0,
-    theta_pd=20, theta_pd_n=40, theta_pd_nsigma=3,
-    phi_pd=.1, phi_pd_n=1, phi_pd_nsigma=0,
-    psi_pd=30, psi_pd_n=1, psi_pd_nsigma=0,
-    req_minor_pd=.1, req_minor_pd_n=1, req_minor_pd_nsigma=0,
-    req_major_pd=.1, req_major_pd_n=1, req_major_pd_nsigma=0,
-    rpolar_pd=.1, rpolar_pd_n=1, rpolar_pd_nsigma=0, dtype=dtype)
+    model = sas.Model(kernel,
+        scale=0.08, req_minor=15, req_major=20, rpolar=500,
+        sldEll=7.105, solvent_sld=.291,
+        background=5, theta=0, phi=phi, psi=0,
+        theta_pd=20, theta_pd_n=40, theta_pd_nsigma=3,
+        phi_pd=.1, phi_pd_n=1, phi_pd_nsigma=0,
+        psi_pd=30, psi_pd_n=1, psi_pd_nsigma=0,
+        req_minor_pd=.1, req_minor_pd_n=1, req_minor_pd_nsigma=0,
+        req_major_pd=.1, req_major_pd_n=1, req_major_pd_nsigma=0,
+        rpolar_pd=.1, rpolar_pd_n=1, rpolar_pd_nsigma=0,
+        )
 
     # SET THE FITTING PARAMETERS
@@ -166,9 +170,11 @@
 
 model.cutoff = cutoff
+M = sas.Experiment(data=data, model=model)
 if section == "both":
-   tan_model = sas.BumpsModel(tan_data, model.model, model.parameters())
+   tan_model = sas.Model(model.kernel, **model.parameters())
    tan_model.phi = model.phi - 90
    tan_model.cutoff = cutoff
-   problem = FitProblem([model, tan_model])
+   tan_M = sas.Experiment(data=tan_data, model=tan_model)
+   problem = FitProblem([M, tan_M])
 else:
-   problem = FitProblem(model)
+   problem = FitProblem(M)
Index: example/sesansfit.py
===================================================================
--- example/sesansfit.py	(revision c210c9459fbe914e1d2c5138d127fd8187122cd9)
+++ example/sesansfit.py	(revision 346bc886df88ad4f5c0e3441573497ef436dbf18)
@@ -1,19 +1,14 @@
-import numpy as np
 from bumps.names import *
 
 from sasmodels import core, bumps_model
 
-kernel = core.load_model("sphere", dtype='single')
-
-
-
 if True: # fix when data loader exists
 #    from sas.dataloader.readers\
     from sas.dataloader.loader import Loader
-    loader=Loader()
+    loader = Loader()
     filename = 'testsasview1.ses'
-    data=loader.load(filename)
+    data = loader.load(filename)
     if data is None: raise IOError("Could not load file %r"%(filename,))
-    data.x /=10
+    data.x /= 10
 #    print data
 #    data = load_sesans('mydatfile.pz')
@@ -43,7 +38,9 @@
 ##  Sphere parameters
 
+kernel = core.load_model("sphere", dtype='single')
 phi = Parameter(0.1, name="phi")
-model = bumps_model.BumpsModel(data, kernel,
-    scale=phi*(1-phi), sld=7.0, solvent_sld=1.0, radius=radius)
+model = bumps_model.Model(kernel,
+    scale=phi*(1-phi), sld=7.0, solvent_sld=1.0, radius=radius,
+    )
 phi.range(0.001,0.5)
 #model.radius.pmp(40)
@@ -54,16 +51,11 @@
 #model.radius_pd_n=0
 
-if False: # have sans data
-    sansmodel = bumps_model.BumpsModel(sans_data, kernel, **model.parameters())
-    problem = FitProblem([model, sansmodel])
-else:
-    problem = FitProblem(model)
-
-
 ### Tri-Axial Ellipsoid
 #
+#kernel = core.load_model("triaxial_ellipsoid", dtype='single')
 #phi = Parameter(0.1, name='phi')
-#model = sas.BumpsModel(data, kernel,
-#    scale=phi*(1-phi), sld=7.0, solvent_sld=1.0, radius=radius)
+#model = bumps_model.Model(kernel,
+#    scale=phi*(1-phi), sld=7.0, solvent_sld=1.0, radius=radius,
+#    )
 #phi.range(0.001,0.90)
 ##model.radius.pmp(40)
@@ -71,10 +63,13 @@
 ##model.sld.pmp(5)
 ##model.background
-##model.radius_pd=0
-##model.radius_pd_n=0
-#
-#if False: # have sans data
-#    sansmodel = sas.BumpsModel(sans_data, kernel, **model.parameters())
-#    problem = FitProblem([model, sansmodel])
-#else:
-#    problem = FitProblem(model)
+##model.radius_pd = 0
+##model.radius_pd_n = 0
+
+if False: # have sans data
+    M_sesans = bumps_model.Experiment(data=data, model=model)
+    M_sans = bumps_model.Experiment(data=sans_data, model=model)
+    problem = FitProblem([M_sesans, M_sans])
+else:
+    M_sesans = bumps_model.Experiment(data=data, model=model)
+    problem = FitProblem(M_sesans)
+
Index: sasmodels/bumps_model.py
===================================================================
--- sasmodels/bumps_model.py	(revision 33d38d57ac52c5b89fc9d1954f5929f5f365d659)
+++ sasmodels/bumps_model.py	(revision 346bc886df88ad4f5c0e3441573497ef436dbf18)
@@ -1,11 +1,44 @@
 """
 Wrap sasmodels for direct use by bumps.
+
+:class:`Model` is a wrapper for the sasmodels kernel which defines a
+bumps *Parameter* box for each kernel parameter.  *Model* accepts keyword
+arguments to set the initial value for each parameter.
+
+:class:`Experiment` combines the *Model* function with a data file loaded by the
+sasview data loader.  *Experiment* takes a *cutoff* parameter controlling
+how far the polydispersity integral extends.
+
+A variety of helper functions are provided:
+
+    :func:`load_data` loads a sasview data file.
+
+    :func:`empty_data1D` creates an empty dataset, which is useful for plotting
+    a theory function before the data is measured.
+
+    :func:`empty_data2D` creates an empty 2D dataset.
+
+    :func:`set_beam_stop` masks the beam stop from the data.
+
+    :func:`set_half` selects the right or left half of the data, which can
+    be useful for shear measurements which have not been properly corrected
+    for path length and reflections.
+
+    :func:`set_top` cuts the top part off the data.
+
+    :func:`plot_data` plots the data file.
+
+    :func:`plot_theory` plots a calculated result from the model.
+
 """
 
 import datetime
+import warnings
 
 import numpy as np
 
 from . import sesans
+from .resolution import Perfect1D, Pinhole1D, Slit1D
+from .resolution2d import Pinhole2D
 
 # CRUFT python 2.6
@@ -20,4 +53,15 @@
 
 
+# CRUFT: old style bumps wrapper which doesn't separate data and model
+def BumpsModel(data, model, cutoff=1e-5, **kw):
+    warnings.warn("Use of BumpsModel is deprecated.  Use bumps_model.Experiment instead.")
+    model = Model(model, **kw)
+    experiment = Experiment(data=data, model=model, cutoff=cutoff)
+    for k in model._parameter_names:
+        setattr(experiment, k, getattr(model, k))
+    return experiment
+
+
+
 def tic():
     """
@@ -42,10 +86,28 @@
     return data
 
-
-def empty_data1D(q):
+def plot_data(data, view='log'):
+    """
+    Plot data loaded by the sasview loader.
+    """
+    if hasattr(data, 'qx_data'):
+        _plot_2d_signal(data, data.data, view=view)
+    else:
+        # Note: kind of weird using the _plot_result1D to plot just the
+        # data, but it handles the masking and graph markup already, so
+        # do not repeat.
+        _plot_result1D(data, None, None, view)
+
+def plot_theory(data, theory, view='log'):
+    if hasattr(data, 'qx_data'):
+        _plot_2d_signal(data, theory, view=view)
+    else:
+        _plot_result1D(data, theory, None, view, include_data=False)
+
+
+def empty_data1D(q, resolution=0.05):
     """
     Create empty 1D data using the given *q* as the x value.
 
-    Resolutions dq/q is 5%.
+    *resolution* dq/q defaults to 5%.
     """
 
@@ -54,11 +116,12 @@
     Iq = 100 * np.ones_like(q)
     dIq = np.sqrt(Iq)
-    data = Data1D(q, Iq, dx=0.05 * q, dy=dIq)
+    data = Data1D(q, Iq, dx=resolution * q, dy=dIq)
     data.filename = "fake data"
     data.qmin, data.qmax = q.min(), q.max()
+    data.mask = np.zeros(len(Iq), dtype='bool')
     return data
 
 
-def empty_data2D(qx, qy=None):
+def empty_data2D(qx, qy=None, resolution=0.05):
     """
     Create empty 2D data using the given mesh.
@@ -66,5 +129,5 @@
     If *qy* is missing, create a square mesh with *qy=qx*.
 
-    Resolution dq/q is 5%.
+    *resolution* dq/q defaults to 5%.
     """
     from sas.dataloader.data_info import Data2D, Detector
@@ -85,8 +148,11 @@
     data.err_data = dIq
     data.mask = mask
+    data.qmin = 1e-16
+    data.qmax = np.inf
 
     # 5% dQ/Q resolution
-    data.dqx_data = 0.05 * Qx
-    data.dqy_data = 0.05 * Qy
+    if resolution != 0:
+        data.dqx_data = resolution * Qx
+        data.dqy_data = resolution * Qy
 
     detector = Detector()
@@ -145,48 +211,5 @@
 
 
-def plot_data(data, Iq, vmin=None, vmax=None, view='log'):
-    """
-    Plot the target value for the data.  This could be the data itself,
-    the theory calculation, or the residuals.
-
-    *scale* can be 'log' for log scale data, or 'linear'.
-    """
-    from numpy.ma import masked_array
-    import matplotlib.pyplot as plt
-    if hasattr(data, 'qx_data'):
-        Iq = Iq + 0
-        valid = np.isfinite(Iq)
-        if view == 'log':
-            valid[valid] = (Iq[valid] > 0)
-            Iq[valid] = np.log10(Iq[valid])
-        elif view == 'q4':
-            Iq[valid] = Iq*(data.qx_data[valid]**2+data.qy_data[valid]**2)**2
-        Iq[~valid | data.mask] = 0
-        #plottable = Iq
-        plottable = masked_array(Iq, ~valid | data.mask)
-        xmin, xmax = min(data.qx_data), max(data.qx_data)
-        ymin, ymax = min(data.qy_data), max(data.qy_data)
-        try:
-            if vmin is None: vmin = Iq[valid & ~data.mask].min()
-            if vmax is None: vmax = Iq[valid & ~data.mask].max()
-        except:
-            vmin, vmax = 0, 1
-        plt.imshow(plottable.reshape(128, 128),
-                   interpolation='nearest', aspect=1, origin='upper',
-                   extent=[xmin, xmax, ymin, ymax], vmin=vmin, vmax=vmax)
-    else: # 1D data
-        if view == 'linear' or view == 'q4':
-            #idx = np.isfinite(Iq)
-            scale = data.x**4 if view == 'q4' else 1.0
-            plt.plot(data.x, scale*Iq) #, '.')
-        else:
-            # Find the values that are finite and positive
-            idx = np.isfinite(Iq)
-            idx[idx] = (Iq[idx] > 0)
-            Iq[~idx] = np.nan
-            plt.loglog(data.x, Iq)
-
-
-def _plot_result1D(data, theory, view):
+def _plot_result1D(data, theory, resid, view, include_data=True):
     """
     Plot the data and residuals for 1D data.
@@ -197,22 +220,32 @@
     #data.y[data.y<0.05] = 0.5
     mdata = masked_array(data.y, data.mask)
-    mdata[np.isnan(mdata)] = masked
+    mdata[~np.isfinite(mdata)] = masked
     if view is 'log':
         mdata[mdata <= 0] = masked
-    mtheory = masked_array(theory, mdata.mask)
-    mresid = masked_array((theory - data.y) / data.dy, mdata.mask)
 
     scale = data.x**4 if view == 'q4' else 1.0
-    plt.subplot(121)
-    plt.errorbar(data.x, scale*mdata, yerr=data.dy)
-    plt.plot(data.x, scale*mtheory, '-', hold=True)
+    if resid is not None:
+        plt.subplot(121)
+
+    if include_data:
+        plt.errorbar(data.x, scale*mdata, yerr=data.dy, fmt='.')
+    if theory is not None:
+        mtheory = masked_array(theory, mdata.mask)
+        plt.plot(data.x, scale*mtheory, '-', hold=True)
+    plt.xscale(view)
     plt.yscale('linear' if view == 'q4' else view)
-    plt.subplot(122)
-    plt.plot(data.x, mresid, 'x')
+    plt.xlabel('Q')
+    plt.ylabel('I(Q)')
+    if resid is not None:
+        mresid = masked_array(resid, mdata.mask)
+        plt.subplot(122)
+        plt.plot(data.x, mresid, 'x')
+        plt.ylabel('residuals')
+        plt.xlabel('Q')
+        plt.xscale(view)
 
 # pylint: disable=unused-argument
-def _plot_sesans(data, theory, view):
+def _plot_sesans(data, theory, resid, view):
     import matplotlib.pyplot as plt
-    resid = (theory - data.y) / data.dy
     plt.subplot(121)
     plt.errorbar(data.x, data.y, yerr=data.dy)
@@ -225,86 +258,75 @@
     plt.ylabel('residuals (P/P0)')
 
-def _plot_result2D(data, theory, view):
+def _plot_result2D(data, theory, resid, view):
     """
     Plot the data and residuals for 2D data.
     """
     import matplotlib.pyplot as plt
-    resid = (theory - data.data) / data.err_data
+    target = data.data[~data.mask]
+    if view == 'log':
+        vmin = min(target[target>0].min(), theory[theory>0].min())
+        vmax = max(target.max(), theory.max())
+    else:
+        vmin = min(target.min(), theory.min())
+        vmax = max(target.max(), theory.max())
+    #print vmin, vmax
     plt.subplot(131)
-    plot_data(data, data.data, view=view)
+    _plot_2d_signal(data, target, view=view, vmin=vmin, vmax=vmax)
+    plt.title('data')
     plt.colorbar()
     plt.subplot(132)
-    plot_data(data, theory, view=view)
+    _plot_2d_signal(data, theory, view=view, vmin=vmin, vmax=vmax)
+    plt.title('theory')
     plt.colorbar()
     plt.subplot(133)
-    plot_data(data, resid, view='linear')
+    _plot_2d_signal(data, resid, view='linear')
+    plt.title('residuals')
     plt.colorbar()
 
-class BumpsModel(object):
-    """
-    Return a bumps wrapper for a SAS model.
-
-    *data* is the data to be fitted.
-
-    *model* is the SAS model from :func:`core.load_model`.
-
-    *cutoff* is the integration cutoff, which avoids computing the
-    the SAS model where the polydispersity weight is low.
-
-    Model parameters can be initialized with additional keyword
-    arguments, or by assigning to model.parameter_name.value.
-
-    The resulting bumps model can be used directly in a FitProblem call.
-    """
-    def __init__(self, data, model, cutoff=1e-5, **kw):
+def _plot_2d_signal(data, signal, vmin=None, vmax=None, view='log'):
+    """
+    Plot the target value for the data.  This could be the data itself,
+    the theory calculation, or the residuals.
+
+    *scale* can be 'log' for log scale data, or 'linear'.
+    """
+    import matplotlib.pyplot as plt
+    from numpy.ma import masked_array
+
+    image = np.zeros_like(data.qx_data)
+    image[~data.mask] = signal
+    valid = np.isfinite(image)
+    if view == 'log':
+        valid[valid] = (image[valid] > 0)
+        image[valid] = np.log10(image[valid])
+    elif view == 'q4':
+        image[valid] *= (data.qx_data[valid]**2+data.qy_data[valid]**2)**2
+    image[~valid | data.mask] = 0
+    #plottable = Iq
+    plottable = masked_array(image, ~valid | data.mask)
+    xmin, xmax = min(data.qx_data), max(data.qx_data)
+    ymin, ymax = min(data.qy_data), max(data.qy_data)
+    # TODO: fix vmin, vmax so it is shared for theory/resid
+    vmin = vmax = None
+    try:
+        if vmin is None: vmin = image[valid & ~data.mask].min()
+        if vmax is None: vmax = image[valid & ~data.mask].max()
+    except:
+        vmin, vmax = 0, 1
+    #print vmin,vmax
+    plt.imshow(plottable.reshape(128, 128),
+               interpolation='nearest', aspect=1, origin='upper',
+               extent=[xmin, xmax, ymin, ymax], vmin=vmin, vmax=vmax)
+
+
+class Model(object):
+    def __init__(self, kernel, **kw):
         from bumps.names import Parameter
 
-        # remember inputs so we can inspect from outside
-        self.data = data
-        self.model = model
-        self.cutoff = cutoff
-        if hasattr(data, 'lam'):
-            self.data_type = 'sesans'
-        elif hasattr(data, 'qx_data'):
-            self.data_type = 'Iqxy'
-        else:
-            self.data_type = 'Iq'
-
-        partype = model.info['partype']
-
-        # interpret data
-        if self.data_type == 'sesans':
-            q = sesans.make_q(data.sample.zacceptance, data.Rmax)
-            self.index = slice(None, None)
-            self.Iq = data.y
-            self.dIq = data.dy
-            self._theory = np.zeros_like(q)
-            q_vectors = [q]
-        elif self.data_type == 'Iqxy':
-            self.index = (data.mask == 0) & (~np.isnan(data.data))
-            self.Iq = data.data[self.index]
-            self.dIq = data.err_data[self.index]
-            self._theory = np.zeros_like(data.data)
-            if not partype['orientation'] and not partype['magnetic']:
-                q_vectors = [np.sqrt(data.qx_data ** 2 + data.qy_data ** 2)]
-            else:
-                q_vectors = [data.qx_data, data.qy_data]
-        elif self.data_type == 'Iq':
-            self.index = (data.x >= data.qmin) & (data.x <= data.qmax) & ~np.isnan(data.y)
-            self.Iq = data.y[self.index]
-            self.dIq = data.dy[self.index]
-            self._theory = np.zeros_like(data.y)
-            q_vectors = [data.x]
-        else:
-            raise ValueError("Unknown data type") # never gets here
-
-        # Remember function inputs so we can delay loading the function and
-        # so we can save/restore state
-        self._fn_inputs = [v[self.index] for v in q_vectors]
-        self._fn = None
-
-        # define bumps parameters
+        self.kernel = kernel
+        partype = kernel.info['partype']
+
         pars = []
-        for p in model.info['parameters']:
+        for p in kernel.info['parameters']:
             name, default, limits = p[0], p[2], p[3]
             value = kw.pop(name, default)
@@ -313,8 +335,8 @@
         for name in partype['pd-2d']:
             for xpart, xdefault, xlimits in [
-                    ('_pd', 0, limits),
-                    ('_pd_n', 35, (0, 1000)),
-                    ('_pd_nsigma', 3, (0, 10)),
-                    ('_pd_type', 'gaussian', None),
+                ('_pd', 0, limits),
+                ('_pd_n', 35, (0, 1000)),
+                ('_pd_nsigma', 3, (0, 10)),
+                ('_pd_type', 'gaussian', None),
                 ]:
                 xname = name + xpart
@@ -328,4 +350,96 @@
             raise TypeError("unexpected parameters: %s"
                             % (", ".join(sorted(kw.keys()))))
+
+    def parameters(self):
+        """
+        Return a dictionary of parameters
+        """
+        return dict((k, getattr(self, k)) for k in self._parameter_names)
+
+class Experiment(object):
+    """
+    Return a bumps wrapper for a SAS model.
+
+    *data* is the data to be fitted.
+
+    *model* is the SAS model from :func:`core.load_model`.
+
+    *cutoff* is the integration cutoff, which avoids computing the
+    the SAS model where the polydispersity weight is low.
+
+    Model parameters can be initialized with additional keyword
+    arguments, or by assigning to model.parameter_name.value.
+
+    The resulting bumps model can be used directly in a FitProblem call.
+    """
+    def __init__(self, data, model, cutoff=1e-5):
+
+        # remember inputs so we can inspect from outside
+        self.data = data
+        self.model = model
+        self.cutoff = cutoff
+        if hasattr(data, 'lam'):
+            self.data_type = 'sesans'
+        elif hasattr(data, 'qx_data'):
+            self.data_type = 'Iqxy'
+        else:
+            self.data_type = 'Iq'
+
+        # interpret data
+        partype = model.kernel.info['partype']
+        if self.data_type == 'sesans':
+            q = sesans.make_q(data.sample.zacceptance, data.Rmax)
+            self.index = slice(None, None)
+            self.Iq = data.y
+            self.dIq = data.dy
+            #self._theory = np.zeros_like(q)
+            q_vectors = [q]
+        elif self.data_type == 'Iqxy':
+            q = np.sqrt(data.qx_data**2 + data.qy_data**2)
+            qmin = getattr(data, 'qmin', 1e-16)
+            qmax = getattr(data, 'qmax', np.inf)
+            self.index = (~data.mask) & (~np.isnan(data.data)) \
+                         & (q >= qmin) & (q <= qmax)
+            self.Iq = data.data[self.index]
+            self.dIq = data.err_data[self.index]
+            self.resolution = Pinhole2D(data=data, index=self.index,
+                                        nsigma=3.0, accuracy='Low')
+            #self._theory = np.zeros_like(self.Iq)
+            if not partype['orientation'] and not partype['magnetic']:
+                raise ValueError("not 2D without orientation or magnetic parameters")
+                #qx,qy = self.resolution.q_calc
+                #q_vectors = [np.sqrt(qx**2 + qy**2)]
+            else:
+                q_vectors = self.resolution.q_calc
+        elif self.data_type == 'Iq':
+            self.index = (data.x >= data.qmin) & (data.x <= data.qmax) & ~np.isnan(data.y)
+            self.Iq = data.y[self.index]
+            self.dIq = data.dy[self.index]
+            if getattr(data, 'dx', None) is not None:
+                q, dq = data.x[self.index], data.dx[self.index]
+                if (dq>0).any():
+                    self.resolution = Pinhole1D(q, dq)
+                else:
+                    self.resolution = Perfect1D(q)
+            elif (getattr(data, 'dxl', None) is not None and
+                  getattr(data, 'dxw', None) is not None):
+                q = data.x[self.index]
+                width = data.dxh[self.index]  # Note: dx
+                self.resolution = Slit1D(data.x[self.index],
+                                         width=data.dxh[self.index],
+                                         height=data.dxw[self.index])
+            else:
+                self.resolution = Perfect1D(data.x[self.index])
+
+            #self._theory = np.zeros_like(self.Iq)
+            q_vectors = [self.resolution.q_calc]
+        else:
+            raise ValueError("Unknown data type") # never gets here
+
+        # Remember function inputs so we can delay loading the function and
+        # so we can save/restore state
+        self._fn_inputs = [v for v in q_vectors]
+        self._fn = None
+
         self.update()
 
@@ -341,31 +455,32 @@
     def parameters(self):
         """
-            Return a dictionary of parameters
-        """
-        return dict((k, getattr(self, k)) for k in self._parameter_names)
+        Return a dictionary of parameters
+        """
+        return self.model.parameters()
 
     def theory(self):
         if 'theory' not in self._cache:
             if self._fn is None:
-                q_input = self.model.make_input(self._fn_inputs)
-                self._fn = self.model(q_input)
-
-            fixed_pars = [getattr(self, p).value for p in self._fn.fixed_pars]
+                q_input = self.model.kernel.make_input(self._fn_inputs)
+                self._fn = self.model.kernel(q_input)
+
+            fixed_pars = [getattr(self.model, p).value for p in self._fn.fixed_pars]
             pd_pars = [self._get_weights(p) for p in self._fn.pd_pars]
             #print fixed_pars,pd_pars
-            self._theory[self.index] = self._fn(fixed_pars, pd_pars, self.cutoff)
+            Iq_calc = self._fn(fixed_pars, pd_pars, self.cutoff)
             #self._theory[:] = self._fn.eval(pars, pd_pars)
             if self.data_type == 'sesans':
                 result = sesans.hankel(self.data.x, self.data.lam * 1e-9,
                                        self.data.sample.thickness / 10,
-                                       self._fn_inputs[0], self._theory)
+                                       self._fn_inputs[0], Iq_calc)
                 self._cache['theory'] = result
             else:
-                self._cache['theory'] = self._theory
+                Iq = self.resolution.apply(Iq_calc)
+                self._cache['theory'] = Iq
         return self._cache['theory']
 
     def residuals(self):
         #if np.any(self.err ==0): print "zeros in err"
-        return (self.theory()[self.index] - self.Iq) / self.dIq
+        return (self.theory() - self.Iq) / self.dIq
 
     def nllf(self):
@@ -381,24 +496,23 @@
         Plot the data and residuals.
         """
-        data, theory = self.data, self.theory()
+        data, theory, resid = self.data, self.theory(), self.residuals()
         if self.data_type == 'Iq':
-            _plot_result1D(data, theory, view)
+            _plot_result1D(data, theory, resid, view)
         elif self.data_type == 'Iqxy':
-            _plot_result2D(data, theory, view)
+            _plot_result2D(data, theory, resid, view)
         elif self.data_type == 'sesans':
-            _plot_sesans(data, theory, view)
+            _plot_sesans(data, theory, resid, view)
         else:
             raise ValueError("Unknown data type")
 
     def simulate_data(self, noise=None):
-        print "noise", noise
-        if noise is None:
-            noise = self.dIq[self.index]
-        else:
-            noise = 0.01 * noise
-            self.dIq[self.index] = noise
-        y = self.theory()
-        y += y * np.random.randn(*y.shape) * noise
+        theory = self.theory()
+        if noise is not None:
+            self.dIq = theory*noise*0.01
+        dy = self.dIq
+        y = theory + np.random.randn(*dy.shape) * dy
+        self.Iq = y
         if self.data_type == 'Iq':
+            self.data.dy[self.index] = dy
             self.data.y[self.index] = y
         elif self.data_type == 'Iqxy':
@@ -418,8 +532,8 @@
         from . import weights
 
-        relative = self.model.info['partype']['pd-rel']
-        limits = self.model.info['limits']
+        relative = self.model.kernel.info['partype']['pd-rel']
+        limits = self.model.kernel.info['limits']
         disperser, value, npts, width, nsigma = [
-            getattr(self, par + ext)
+            getattr(self.model, par + ext)
             for ext in ('_pd_type', '', '_pd_n', '_pd', '_pd_nsigma')]
         value, weight = weights.get_weights(
@@ -442,5 +556,5 @@
     data = load_data('DEC07086.DAT')
     set_beam_stop(data, 0.004)
-    plot_data(data, data.data)
+    plot_data(data)
     import matplotlib.pyplot as plt; plt.show()
 
Index: sasmodels/resolution.py
===================================================================
--- sasmodels/resolution.py	(revision f1b8c901d197c80b8cd672fa3a0346c6ec2daf9d)
+++ sasmodels/resolution.py	(revision 346bc886df88ad4f5c0e3441573497ef436dbf18)
@@ -11,5 +11,5 @@
 MINIMUM_RESOLUTION = 1e-8
 
-class Resolution1D(object):
+class Resolution(object):
     """
     Abstract base class defining a 1D resolution function.
@@ -32,5 +32,5 @@
 
 
-class Perfect1D(Resolution1D):
+class Perfect1D(Resolution):
     """
     Resolution function to use when there is no actual resolution smearing
@@ -45,5 +45,5 @@
 
 
-class Pinhole1D(Resolution1D):
+class Pinhole1D(Resolution):
     r"""
     Pinhole aperture with q-dependent gaussian resolution.
@@ -74,8 +74,10 @@
 
     def apply(self, theory):
+        print "q calc", self.q_calc
+        print "weights", self.weight_matrix.shape
         return apply_resolution_matrix(self.weight_matrix, theory)
 
 
-class Slit1D(Resolution1D):
+class Slit1D(Resolution):
     """
     Slit aperture with a complicated resolution function.
Index: sasmodels/resolution2d.py
===================================================================
--- sasmodels/resolution2d.py	(revision 346bc886df88ad4f5c0e3441573497ef436dbf18)
+++ sasmodels/resolution2d.py	(revision 346bc886df88ad4f5c0e3441573497ef436dbf18)
@@ -0,0 +1,236 @@
+"""
+#This software was developed by the University of Tennessee as part of the
+#Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
+#project funded by the US National Science Foundation. 
+#See the license text in license.txt
+"""
+from __future__ import division
+
+import numpy as np
+from numpy import pi, cos, sin, sqrt
+
+from .resolution import Resolution
+
+## Singular point
+SIGMA_ZERO = 1.0e-010
+## Limit of how many sigmas to be covered for the Gaussian smearing
+# default: 2.5 to cover 98.7% of Gaussian
+NSIGMA = 3.0
+## Defaults
+NR = {'Xhigh':10, 'High':5, 'Med':5, 'Low':3}
+NPHI ={'Xhigh':20, 'High':12, 'Med':6, 'Low':4}
+
+class Pinhole2D(Resolution):
+    """
+    Gaussian Q smearing class for SAS 2d data
+    """
+     
+    def __init__(self, data=None, index=None,
+                 nsigma=NSIGMA, accuracy='Low', coords='polar'):
+        """
+        Assumption: equally spaced bins in dq_r, dq_phi space.
+        
+        :param data: 2d data used to set the smearing parameters
+        :param index: 1d array with len(data) to define the range
+         of the calculation: elements are given as True or False
+        :param nr: number of bins in dq_r-axis
+        :param nphi: number of bins in dq_phi-axis
+        :param coord: coordinates [string], 'polar' or 'cartesian'
+        """
+        ## Accuracy: Higher stands for more sampling points in both directions
+        ## of r and phi.
+        ## number of bins in r axis for over-sampling
+        self.nr = NR[accuracy]
+        ## number of bins in phi axis for over-sampling 
+        self.nphi = NPHI[accuracy]
+        ## maximum nsigmas
+        self.nsigma= nsigma
+        self.coords = coords
+        self._init_data(data, index)
+
+    def _init_data(self, data, index):
+        """
+        Get qx_data, qy_data, dqx_data,dqy_data,
+        and calculate phi_data=arctan(qx_data/qy_data)
+        """
+        # TODO: maybe don't need to hold copy of qx,qy,dqx,dqy,data,index
+        # just need q_calc and weights
+        self.data = data
+        self.index = index
+
+        self.qx_data = data.qx_data[index]
+        self.qy_data = data.qy_data[index]
+        self.q_data = data.q_data[index]
+
+        dqx = getattr(data, 'dqx_data', None)
+        dqy = getattr(data, 'dqy_data', None)
+        if dqx is not None and dqy is not None:
+            # Here dqx and dqy mean dq_parr and dq_perp
+            self.dqx_data = dqx[index]
+            self.dqy_data = dqy[index]
+            ## Remove singular points if exists
+            self.dqx_data[self.dqx_data < SIGMA_ZERO] = SIGMA_ZERO
+            self.dqy_data[self.dqy_data < SIGMA_ZERO] = SIGMA_ZERO
+            qx_calc, qy_calc, weights = self._calc_res()
+            self.q_calc = [qx_calc, qy_calc]
+            self.q_calc_weights = weights
+        else:
+            # No resolution information
+            self.dqx_data = self.dqy_data = None
+            self.q_calc = [self.qx_data, self.qy_data]
+            self.q_calc_weights = None
+
+        #self.phi_data = np.arctan(self.qx_data / self.qy_data)
+
+    def _calc_res(self):
+        """
+        Over sampling of r_nbins times phi_nbins, calculate Gaussian weights, 
+        then find smeared intensity
+        """    
+        nr, nphi = self.nr, self.nphi
+        # Total number of bins = # of bins
+        nbins = nr * nphi
+        # Number of bins in the dqr direction (polar coordinate of dqx and dqy)
+        bin_size = self.nsigma // nr
+        # in dq_r-direction times # of bins in dq_phi-direction
+        # data length in the range of self.index
+        nq = len(self.qx_data)
+
+        # Mean values of dqr at each bins
+        # starting from the half of bin size
+        r = bin_size / 2.0 + np.arange(nr) * bin_size
+        # mean values of qphi at each bines
+        phi = np.arange(nphi)
+        dphi = phi * 2.0 * pi / nphi
+        dphi = dphi.repeat(nr)
+
+        ## Transform to polar coordinate,
+        #  and set dphi at each data points ; 1d array
+        dphi = dphi.repeat(nq)
+        q_phi = self.qy_data / self.qx_data
+
+        # Starting angle is different between polar
+        #  and cartesian coordinates.
+        #if self.coords != 'polar':
+        #    dphi += np.arctan( q_phi * self.dqx_data/ \
+        #                  self.dqy_data).repeat(nbins).reshape(nq,\
+        #                                nbins).transpose().flatten()
+
+        # The angle (phi) of the original q point
+        q_phi = np.arctan(q_phi).repeat(nbins)\
+            .reshape(nq, nbins).transpose().flatten()
+        ## Find Gaussian weight for each dq bins: The weight depends only
+        #  on r-direction (The integration may not need)
+        weight_res = (np.exp(-0.5 * (r - bin_size / 2.0)**2)  -
+                      np.exp(-0.5 * (r + bin_size / 2.0)**2))
+        # No needs of normalization here.
+        #weight_res /= np.sum(weight_res)
+        weight_res = weight_res.repeat(nphi).reshape(nr, nphi)
+        weight_res = weight_res.transpose().flatten()
+
+        ## Set dr for all dq bins for averaging
+        dr = r.repeat(nphi).reshape(nr, nphi).transpose().flatten()
+        ## Set dqr for all data points
+        dqx = np.outer(dr, self.dqx_data).flatten()
+        dqy = np.outer(dr, self.dqy_data).flatten()
+
+        qx = self.qx_data.repeat(nbins)\
+            .reshape(nq, nbins).transpose().flatten()
+        qy = self.qy_data.repeat(nbins)\
+            .reshape(nq, nbins).transpose().flatten()
+
+        # The polar needs rotation by -q_phi
+        if self.coords == 'polar':
+            q_r = sqrt(qx**2 + qy**2)
+            qx_res = ( (dqx*cos(dphi) + q_r) * cos(-q_phi) +
+                           dqy*sin(dphi) * sin(-q_phi))
+            qy_res = (-(dqx*cos(dphi) + q_r) * sin(-q_phi) +
+                           dqy*sin(dphi) * cos(-q_phi))
+        else:
+            qx_res = qx +  dqx*cos(dphi)
+            qy_res = qy +  dqy*sin(dphi)
+
+
+        return qx_res, qy_res, weight_res
+
+    def apply(self, theory):
+        if self.q_calc_weights is not None:
+            # TODO: interpolate rather than recomputing all the different qx,qy
+            # Resolution needs to be applied
+            nq, nbins = len(self.qx_data), self.nr * self.nphi
+            ## Reshape into 2d array to use np weighted averaging
+            theory = np.reshape(theory, (nbins, nq))
+            ## Averaging with Gaussian weighting: normalization included.
+            value =np.average(theory, axis=0, weights=self.q_calc_weights)
+            ## Return the smeared values in the range of self.index
+            return value
+        else:
+            return theory
+
+"""
+if __name__ == '__main__':
+    ## Test w/ 2D linear function
+    x = 0.001*np.arange(1, 11)
+    dx = np.ones(len(x))*0.0003
+    y = 0.001*np.arange(1, 11)
+    dy = np.ones(len(x))*0.001
+    z = np.ones(10)
+    dz = sqrt(z)
+
+    from sas.dataloader import Data2D
+    #for i in range(10): print i, 0.001 + i*0.008/9.0
+    #for i in range(100): print i, int(math.floor( (i/ (100/9.0)) ))
+    out = Data2D()
+    out.data = z
+    out.qx_data = x
+    out.qy_data = y
+    out.dqx_data = dx
+    out.dqy_data = dy
+    out.q_data = sqrt(dx * dx + dy * dy)
+    index = np.ones(len(x), dtype = bool)
+    out.mask = index
+    from sas.models.LineModel import LineModel
+    model = LineModel()
+    model.setParam("A", 0)
+
+    smear = Smearer2D(out, model, index)
+    #smear.set_accuracy('Xhigh')
+    value = smear.get_value()
+    ## All data are ones, so the smeared should also be ones.
+    print "Data length =", len(value)
+    print " 2D linear function, I = 0 + 1*qy"
+    text = " Gaussian weighted averaging on a 2D linear function will "
+    text += "provides the results same as without the averaging."
+    print text
+    print "qx_data", "qy_data", "I_nonsmear", "I_smeared"
+    for ind in range(len(value)):
+        print x[ind], y[ind], model.evalDistribution([x, y])[ind], value[ind]
+
+
+if __name__ == '__main__':
+    ## Another Test w/ constant function
+    x = 0.001*np.arange(1,11)
+    dx = np.ones(len(x))*0.001
+    y = 0.001*np.arange(1,11)
+    dy = np.ones(len(x))*0.001
+    z = np.ones(10)
+    dz = sqrt(z)
+
+    from DataLoader import Data2D
+    #for i in range(10): print i, 0.001 + i*0.008/9.0
+    #for i in range(100): print i, int(math.floor( (i/ (100/9.0)) ))
+    out = Data2D()
+    out.data = z
+    out.qx_data = x
+    out.qy_data = y
+    out.dqx_data = dx
+    out.dqy_data = dy
+    index = np.ones(len(x), dtype = bool)
+    out.mask = index
+    from sas.models.Constant import Constant
+    model = Constant()
+
+    value = Smearer2D(out,model,index).get_value()
+    ## All data are ones, so the smeared values should also be ones.
+    print "Data length =",len(value), ", Data=",value
+"""
