Index: example/oriented_usans.py
===================================================================
--- example/oriented_usans.py	(revision ea7504385b3e91f10b68265b8e6a5a5519f04260)
+++ example/oriented_usans.py	(revision ea7504385b3e91f10b68265b8e6a5a5519f04260)
@@ -0,0 +1,41 @@
+from bumps.names import *
+
+from sasmodels.data import load_data
+from sasmodels.core import load_model
+from sasmodels.bumps_model import Model, Experiment
+
+# Spherical particle data, not ellipsoids
+sans, usans = load_data('../../sasview/sasview/test/1d_data/latex_smeared.xml')
+usans.qmin, usans.qmax = np.min(usans.x), np.max(usans.x)
+usans.mask = (usans.x < 0.0)
+usans.oriented = True
+#print sans.dxl, usans.dxl
+#import pprint; pprint.pprint(sans.__dict__)
+
+kernel = load_model("ellipsoid")
+
+model = Model(
+    kernel,
+    scale=0.08, background=0,
+    sld=.291, sld_solvent=7.105,
+    r_polar=1800, r_polar_pd=0.222296, r_polar_pd_n=0,
+    r_equatorial=2600, r_equatorial_pd=0.28, r_equatorial_pd_n=0,
+    theta=60, theta_pd=0, theta_pd_n=0,
+    phi=60, phi_pd=0, phi_pd_n=0,
+    )
+
+# SET THE FITTING PARAMETERS
+model.r_polar.range(1000, 10000)
+model.r_equatorial.range(1000, 10000)
+model.theta.range(0, 360)
+model.phi.range(0, 360)
+model.background.range(0,1000)
+model.scale.range(0, 10)
+
+
+#cutoff = 0     # no cutoff on polydisperisity loops
+#cutoff = 1e-5  # default cutoff
+cutoff = 1e-3  # low precision cutoff
+M = Experiment(data=usans, model=model, cutoff=cutoff)
+problem = FitProblem(M)
+
Index: sasmodels/bumps_model.py
===================================================================
--- sasmodels/bumps_model.py	(revision 9217ef8174f5566df0d9259319dc1bde48372242)
+++ sasmodels/bumps_model.py	(revision ea7504385b3e91f10b68265b8e6a5a5519f04260)
@@ -223,5 +223,5 @@
         """
         data, theory, resid = self._data, self.theory(), self.residuals()
-        plot_theory(data, theory, resid, view)
+        plot_theory(data, theory, resid, view, Iq_calc = self.Iq_calc)
 
     def simulate_data(self, noise=None):
Index: sasmodels/compare.py
===================================================================
--- sasmodels/compare.py	(revision fa1582e600d4da130f72b038d8310888e0bb6894)
+++ sasmodels/compare.py	(revision ea7504385b3e91f10b68265b8e6a5a5519f04260)
@@ -473,5 +473,5 @@
         data = empty_data2D(np.linspace(-qmax, qmax, nq), resolution=res)
         data.accuracy = opts['accuracy']
-        set_beam_stop(data, 0.004)
+        set_beam_stop(data, 0.0004)
         index = ~data.mask
     else:
Index: sasmodels/data.py
===================================================================
--- sasmodels/data.py	(revision 1d61d077d41c07deb389a93edc5b1c71b823b3f0)
+++ sasmodels/data.py	(revision ea7504385b3e91f10b68265b8e6a5a5519f04260)
@@ -322,5 +322,5 @@
 
 def plot_theory(data, theory, resid=None, view='log',
-                use_data=True, limits=None):
+                use_data=True, limits=None, Iq_calc=None):
     """
     Plot theory calculation.
@@ -343,5 +343,6 @@
         _plot_result2D(data, theory, resid, view, use_data, limits=limits)
     else:
-        _plot_result1D(data, theory, resid, view, use_data, limits=limits)
+        _plot_result1D(data, theory, resid, view, use_data,
+                       limits=limits, Iq_calc=Iq_calc)
 
 
@@ -366,5 +367,6 @@
 
 @protect
-def _plot_result1D(data, theory, resid, view, use_data, limits=None):
+def _plot_result1D(data, theory, resid, view, use_data,
+                   limits=None, Iq_calc=None):
     """
     Plot the data and residuals for 1D data.
@@ -376,5 +378,7 @@
     use_theory = theory is not None
     use_resid = resid is not None
-    num_plots = (use_data or use_theory) + use_resid
+    use_calc = use_theory and Iq_calc is not None
+    num_plots = (use_data or use_theory) + use_calc + use_resid
+
 
     scale = data.x**4 if view == 'q4' else 1.0
@@ -416,4 +420,13 @@
         plt.ylabel('$I(q)$')
 
+    if use_calc:
+        # Only have use_calc if have use_theory
+        plt.subplot(1, num_plots, 2)
+        qx, qy, Iqxy = Iq_calc
+        plt.pcolormesh(qx, qy, np.log10(Iqxy))
+        plt.xlabel("$q_x$/A$^{-1}$")
+        plt.xlabel("$q_y$/A$^{-1}$")
+        #plt.axis('equal')
+
     if use_resid:
         mresid = masked_array(resid, data.mask.copy())
@@ -422,5 +435,5 @@
 
         if num_plots > 1:
-            plt.subplot(1, num_plots, (use_data or use_theory) + 1)
+            plt.subplot(1, num_plots, use_calc + 2)
         plt.plot(data.x, mresid, '-')
         plt.xlabel("$q$/A$^{-1}$")
@@ -562,13 +575,13 @@
     plottable = masked_array(image, ~valid | data.mask)
     # Divide range by 10 to convert from angstroms to nanometers
-    xmin, xmax = min(data.qx_data)/10, max(data.qx_data)/10
-    ymin, ymax = min(data.qy_data)/10, max(data.qy_data)/10
+    xmin, xmax = min(data.qx_data), max(data.qx_data)
+    ymin, ymax = min(data.qy_data), max(data.qy_data)
     if view == 'log':
         vmin, vmax = np.log10(vmin), np.log10(vmax)
     plt.imshow(plottable.reshape(len(data.x_bins), len(data.y_bins)),
-               interpolation='nearest', aspect=1, origin='upper',
+               interpolation='nearest', aspect=1, origin='lower',
                extent=[xmin, xmax, ymin, ymax], vmin=vmin, vmax=vmax)
-    plt.xlabel("$q_x$/nm$^{-1}$")
-    plt.ylabel("$q_y$/nm$^{-1}$")
+    plt.xlabel("$q_x$/A$^{-1}$")
+    plt.ylabel("$q_y$/A$^{-1}$")
     return vmin, vmax
 
Index: sasmodels/direct_model.py
===================================================================
--- sasmodels/direct_model.py	(revision 02e70ffbd2f37e9a91d5d027972b30069bd05c6b)
+++ sasmodels/direct_model.py	(revision ea7504385b3e91f10b68265b8e6a5a5519f04260)
@@ -63,4 +63,6 @@
         elif hasattr(data, 'qx_data'):
             self.data_type = 'Iqxy'
+        elif getattr(data, 'oriented', False):
+            self.data_type = 'Iq-oriented'
         else:
             self.data_type = 'Iq'
@@ -116,6 +118,6 @@
                   and getattr(data, 'dxw', None) is not None):
                 res = resolution.Slit1D(data.x[index],
-                                        width=data.dxh[index],
-                                        height=data.dxw[index])
+                                        qx_width=data.dxw[index],
+                                        qy_width=data.dxl[index])
             else:
                 res = resolution.Perfect1D(data.x[index])
@@ -123,4 +125,21 @@
             #self._theory = np.zeros_like(self.Iq)
             q_vectors = [res.q_calc]
+            q_mono = []
+        elif self.data_type == 'Iq-oriented':
+            index = (data.x >= data.qmin) & (data.x <= data.qmax)
+            if data.y is not None:
+                index &= ~np.isnan(data.y)
+                Iq = data.y[index]
+                dIq = data.dy[index]
+            else:
+                Iq, dIq = None, None
+            if (getattr(data, 'dxl', None) is None
+                or getattr(data, 'dxw', None) is None):
+                raise ValueError("oriented sample with 1D data needs slit resolution")
+
+            res = resolution2d.Slit2D(data.x[index],
+                                      qx_width=data.dxw[index],
+                                      qy_width=data.dxl[index])
+            q_vectors = res.q_calc
             q_mono = []
         else:
@@ -142,5 +161,5 @@
         y = Iq + np.random.randn(*dy.shape) * dy
         self.Iq = y
-        if self.data_type == 'Iq':
+        if self.data_type in ('Iq', 'Iq-oriented'):
             self._data.dy[self.index] = dy
             self._data.y[self.index] = y
@@ -155,9 +174,13 @@
         if self._kernel is None:
             self._kernel = make_kernel(self._model, self._kernel_inputs)  # pylint: disable=attribute-dedata_type
-            self._kernel_mono = make_kernel(self._model, self._kernel_mono_inputs) if self._kernel_mono_inputs else None
+            self._kernel_mono = (make_kernel(self._model, self._kernel_mono_inputs)
+                                 if self._kernel_mono_inputs else None)
 
         Iq_calc = call_kernel(self._kernel, pars, cutoff=cutoff)
-        Iq_mono = call_kernel(self._kernel_mono, pars, mono=True) if self._kernel_mono_inputs else None
+        # TODO: may want to plot the raw Iq for other than oriented usans
+        self.Iq_calc = None
         if self.data_type == 'sesans':
+            Iq_mono = (call_kernel(self._kernel_mono, pars, mono=True)
+                       if self._kernel_mono_inputs else None)
             result = sesans.transform(self._data,
                                    self._kernel_inputs[0], Iq_calc, 
@@ -165,4 +188,9 @@
         else:
             result = self.resolution.apply(Iq_calc)
+            if hasattr(self.resolution, 'nx'):
+                self.Iq_calc = (
+                    self.resolution.qx_calc, self.resolution.qy_calc,
+                    np.reshape(Iq_calc, (self.resolution.ny, self.resolution.nx))
+                )
         return result        
 
Index: sasmodels/resolution.py
===================================================================
--- sasmodels/resolution.py	(revision 8b935d12a0041bac8658cc0f3268fc22ca46a544)
+++ sasmodels/resolution.py	(revision ea7504385b3e91f10b68265b8e6a5a5519f04260)
@@ -82,6 +82,6 @@
         self.q_calc = (pinhole_extend_q(q, q_width, nsigma=nsigma)
                        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))
+        self.weight_matrix = pinhole_resolution(self.q_calc, self.q,
+                                np.maximum(q_width, MINIMUM_RESOLUTION))
 
     def apply(self, theory):
@@ -91,11 +91,11 @@
 class Slit1D(Resolution):
     """
-    Slit aperture with a complicated resolution function.
+    Slit aperture with resolution function.
 
     *q* points at which the data is measured.
 
-    *qx_width* slit width
-
-    *qy_height* slit height
+    *dqx* slit width in qx
+
+    *dqy* slit height in qy
 
     *q_calc* is the list of points to calculate, or None if this should
@@ -104,25 +104,25 @@
     The *weight_matrix* is computed by :func:`slit1d_resolution`
     """
-    def __init__(self, q, width, height=0., q_calc=None):
-        # Remember what width/height was used even though we won't need them
+    def __init__(self, q, qx_width, qy_width=0., q_calc=None):
+        # Remember what width/dqy was used even though we won't need them
         # after the weight matrix is constructed
-        self.width, self.height = width, height
+        self.qx_width, self.qy_width = qx_width, qy_width
 
         # Allow independent resolution on each point even though it is not
         # needed in practice.
-        if np.isscalar(width):
-            width = np.ones(len(q))*width
+        if np.isscalar(qx_width):
+            qx_width = np.ones(len(q))*qx_width
         else:
-            width = np.asarray(width)
-        if np.isscalar(height):
-            height = np.ones(len(q))*height
+            qx_width = np.asarray(qx_width)
+        if np.isscalar(qy_width):
+            qy_width = np.ones(len(q))*qy_width
         else:
-            height = np.asarray(height)
+            qy_width = np.asarray(qy_width)
 
         self.q = q.flatten()
-        self.q_calc = slit_extend_q(q, width, height) \
+        self.q_calc = slit_extend_q(q, qx_width, qy_width) \
             if q_calc is None else np.sort(q_calc)
         self.weight_matrix = \
-            slit_resolution(self.q_calc, self.q, width, height)
+            slit_resolution(self.q_calc, self.q, qx_width, qy_width)
 
     def apply(self, theory):
@@ -396,5 +396,5 @@
     """
     q = np.sort(q)
-    if q_min < q[0]:
+    if q_min + 2*MINIMUM_RESOLUTION < q[0]:
         if q_min <= 0: q_min = q_min*MIN_Q_SCALE_FOR_NEGATIVE_Q_EXTRAPOLATION
         n_low = np.ceil((q[0]-q_min) / (q[1]-q[0])) if q[1] > q[0] else 15
@@ -402,5 +402,5 @@
     else:
         q_low = []
-    if q_max > q[-1]:
+    if q_max - 2*MINIMUM_RESOLUTION > q[-1]:
         n_high = np.ceil((q_max-q[-1]) / (q[-1]-q[-2])) if q[-1] > q[-2] else 15
         q_high = np.linspace(q[-1], q_max, n_high+1)[1:]
@@ -595,5 +595,5 @@
         Slit smearing with perfect resolution.
         """
-        resolution = Slit1D(self.x, width=0, height=0, q_calc=self.x)
+        resolution = Slit1D(self.x, qx_width=0, qy_width=0, q_calc=self.x)
         theory = self.Iq(resolution.q_calc)
         output = resolution.apply(theory)
@@ -605,5 +605,5 @@
         Slit smearing with height 0.005
         """
-        resolution = Slit1D(self.x, width=0, height=0.005, q_calc=self.x)
+        resolution = Slit1D(self.x, qx_width=0, qy_width=0.005, q_calc=self.x)
         theory = self.Iq(resolution.q_calc)
         output = resolution.apply(theory)
@@ -620,5 +620,5 @@
         """
         q = np.logspace(-4, -1, 10)
-        resolution = Slit1D(q, width=0.2, height=np.inf)
+        resolution = Slit1D(q, qx_width=0.2, qy_width=np.inf)
         theory = 1000*self.Iq(resolution.q_calc**4)
         output = resolution.apply(theory)
@@ -634,5 +634,5 @@
         Slit smearing with width 0.0002
         """
-        resolution = Slit1D(self.x, width=0.0002, height=0, q_calc=self.x)
+        resolution = Slit1D(self.x, qx_width=0.0002, qy_width=0, q_calc=self.x)
         theory = self.Iq(resolution.q_calc)
         output = resolution.apply(theory)
@@ -647,5 +647,5 @@
         Slit smearing with width > 100*height.
         """
-        resolution = Slit1D(self.x, width=0.0002, height=0.000001,
+        resolution = Slit1D(self.x, qx_width=0.0002, qy_width=0.000001,
                             q_calc=self.x)
         theory = self.Iq(resolution.q_calc)
@@ -772,5 +772,5 @@
         data = np.loadtxt(data_string.split('\n')).T
         q, delta_qv, _, answer = data
-        resolution = Slit1D(q, width=delta_qv, height=0)
+        resolution = Slit1D(q, qx_width=delta_qv, qy_width=0)
         output = self._eval_sphere(pars, resolution)
         # TODO: eliminate Igor test since it is too inaccurate to be useful.
@@ -792,5 +792,5 @@
         q_calc = slit_extend_q(interpolate(q, 2*np.pi/radius/20),
                                delta_qv[0], 0.)
-        resolution = Slit1D(q, width=delta_qv, height=0, q_calc=q_calc)
+        resolution = Slit1D(q, qx_width=delta_qv, qy_width=0, q_calc=q_calc)
         output = self._eval_sphere(pars, resolution)
         # TODO: relative error should be lower
@@ -810,5 +810,5 @@
         q = np.logspace(log10(4e-5), log10(2.5e-2), 68)
         width, height = 0.117, 0.
-        resolution = Slit1D(q, width=width, height=height)
+        resolution = Slit1D(q, qx_width=width, qy_width=height)
         answer = romberg_slit_1d(q, width, height, form, pars)
         output = resolution.apply(eval_form(resolution.q_calc, form, pars))
@@ -1067,5 +1067,5 @@
 
     if isinstance(resolution, Slit1D):
-        width, height = resolution.width, resolution.height
+        width, height = resolution.dqx, resolution.dqy
         Iq_romb = romberg_slit_1d(resolution.q, width, height, model, pars)
     else:
Index: sasmodels/resolution2d.py
===================================================================
--- sasmodels/resolution2d.py	(revision 823e6204f35d6c2bf32f27a3a3c137a369a973cf)
+++ sasmodels/resolution2d.py	(revision ea7504385b3e91f10b68265b8e6a5a5519f04260)
@@ -10,4 +10,5 @@
 from numpy import pi, cos, sin, sqrt
 
+from . import resolution
 from .resolution import Resolution
 
@@ -20,4 +21,9 @@
 NR = {'xhigh':10, 'high':5, 'med':5, 'low':3}
 NPHI = {'xhigh':20, 'high':12, 'med':6, 'low':4}
+
+## Defaults
+N_SLIT_PERP = {'xhigh':2000, 'high':1000, 'med':500, 'low':250}
+N_SLIT_PERP_DOC = ", ".join("%s=%d"%(name,value) for value,name in
+                            sorted((v,k) for k,v in N_SLIT_PERP.items()))
 
 class Pinhole2D(Resolution):
@@ -167,2 +173,63 @@
         else:
             return theory
+
+
+class Slit2D(Resolution):
+    """
+    Slit aperture with resolution function on an oriented sample.
+
+    *q* points at which the data is measured.
+
+    *qx_width* slit width in qx
+
+    *qy_width* slit height in qy; current implementation requires a fixed
+    qy_width for all q points.
+
+    *q_calc* is the list of q points to calculate, or None if this
+    should be estimated from the *q* and *qx_width*.
+
+    *accuracy* determines the number of *qy* points to compute for each *q*.
+    The values are stored in sasmodels.resolution2d.N_SLIT_PERP.  The default
+    values are: %s
+    """
+    __doc__ = __doc__%N_SLIT_PERP_DOC
+    def __init__(self, q, qx_width, qy_width=0., q_calc=None, accuracy='low'):
+        # Remember what q and width was used even though we won't need them
+        # after the weight matrix is constructed
+        self.q, self.qx_width, self.qy_width = q, qx_width, qy_width
+
+        # Allow independent resolution on each qx point even though it is not
+        # needed in practice.  Set qy_width to the maximum qy width.
+        if np.isscalar(qx_width):
+            qx_width = np.ones(len(q))*qx_width
+        else:
+            qx_width = np.asarray(qx_width)
+        if not np.isscalar(qy_width):
+            qy_width = np.max(qy_width)
+
+        # Build grid of qx, qy points
+        if q_calc is not None:
+            qx_calc = np.sort(q_calc)
+        else:
+            qx_calc = resolution.pinhole_extend_q(q, qx_width, nsigma=3)
+        qy_calc = np.linspace(-qy_width, qy_width, N_SLIT_PERP[accuracy]*10)
+        self.q_calc = [v.flatten() for v in np.meshgrid(qx_calc, qy_calc)]
+        self.qx_calc, self.qy_calc = qx_calc, qy_calc
+        self.nx, self.ny = len(qx_calc), len(qy_calc)
+        self.dy = 2*qy_width/self.ny
+
+        # Build weight matrix for resolution integration
+        if np.any(qx_width > 0):
+            self.weights = resolution.pinhole_resolution(qx_calc, q,
+                    np.maximum(qx_width, resolution.MINIMUM_RESOLUTION))
+        elif len(qx_calc)==len(q) and np.all(qx_calc == q):
+            self.weights = None
+        else:
+            raise ValueError("Slit2D fails with q_calc != q")
+
+    def apply(self, theory):
+        Iq = np.sum(theory.reshape(self.ny, self.nx), axis=0)*self.dy
+        if self.weights is not None:
+            Iq = resolution.apply_resolution_matrix(self.weights, Iq)
+        return Iq
+
