Index: example/fit.py
===================================================================
--- example/fit.py	(revision 7cf2cfd8695ecfd2684525ab15c14bc00f23616a)
+++ example/fit.py	(revision 1182da5952c5f1fedd003e65e6b154df780d0afb)
@@ -31,11 +31,11 @@
     model = Model(kernel,
         scale=0.08,
-        rpolar=15, requatorial=800,
-        sld=.291, solvent_sld=7.105,
+        r_polar=15, r_equatorial=800,
+        sld=.291, sld_solvent=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,
+        r_polar_pd=0.222296, r_polar_pd_n=1, r_polar_pd_nsigma=0,
+        r_equatorial_pd=.000128, r_equatorial_pd_n=1, r_equatorial_pd_nsigma=0,
         phi_pd=0, phi_pd_n=20, phi_pd_nsigma=3,
         )
@@ -43,6 +43,6 @@
 
     # SET THE FITTING PARAMETERS
-    model.rpolar.range(15, 1000)
-    model.requatorial.range(15, 1000)
+    model.r_polar.range(15, 1000)
+    model.r_equatorial.range(15, 1000)
     model.theta_pd.range(0, 360)
     model.background.range(0,1000)
@@ -81,5 +81,5 @@
     pars = dict(
         scale=.01, background=35,
-        sld=.291, solvent_sld=5.77, 
+        sld=.291, sld_solvent=5.77,
         radius=250, length=178, 
         theta=90, phi=phi,
@@ -105,5 +105,5 @@
     model = Model(kernel,
         scale= .031, radius=19.5, thickness=30, length=22,
-        core_sld=7.105, shell_sld=.291, solvent_sld=7.105,
+        sld_core=7.105, sld_shell=.291, sdl_solvent=7.105,
         background=0, theta=0, phi=phi,
 
@@ -132,5 +132,5 @@
     model = Model(kernel,
         scale=.08, radius=20, cap_radius=40, length=400,
-        sld_capcyl=1, sld_solv=6.3,
+        sld=1, sld_solvent=6.3,
         background=0, theta=0, phi=phi,
         radius_pd=.1, radius_pd_n=5, radius_pd_nsigma=3,
@@ -147,5 +147,5 @@
     model = Model(kernel,
         scale=0.08, req_minor=15, req_major=20, rpolar=500,
-        sldEll=7.105, solvent_sld=.291,
+        sld=7.105, solvent_sld=.291,
         background=5, theta=0, phi=phi, psi=0,
         theta_pd=20, theta_pd_n=40, theta_pd_nsigma=3,
Index: example/model.py
===================================================================
--- example/model.py	(revision 15d2285669b798a65396200f7142e22014372cc8)
+++ example/model.py	(revision 1182da5952c5f1fedd003e65e6b154df780d0afb)
@@ -17,17 +17,17 @@
 model = Model(kernel,
     scale=0.08,
-    rpolar=15, requatorial=800,
-    sld=.291, solvent_sld=7.105,
+    r_polar=15, r_equatorial=800,
+    sld=.291, sld_solvent=7.105,
     background=0,
     theta=90, phi=0,
     theta_pd=15, theta_pd_n=40, theta_pd_nsigma=3,
-    rpolar_pd=0.222296, rpolar_pd_n=1, rpolar_pd_nsigma=0,
-    requatorial_pd=.000128, requatorial_pd_n=1, requatorial_pd_nsigma=0,
+    r_polar_pd=0.222296, r_polar_pd_n=1, r_polar_pd_nsigma=0,
+    r_equatorial_pd=.000128, r_equatorial_pd_n=1, r_equatorial_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.r_polar.range(15, 1000)
+model.r_equatorial.range(15, 1000)
 model.theta_pd.range(0, 360)
 model.background.range(0,1000)
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 d19962c0471a0eee538ce03020b5b226ba2752f6)
+++ sasmodels/bumps_model.py	(revision ea7504385b3e91f10b68265b8e6a5a5519f04260)
@@ -219,5 +219,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 d19962c0471a0eee538ce03020b5b226ba2752f6)
+++ sasmodels/compare.py	(revision ea7504385b3e91f10b68265b8e6a5a5519f04260)
@@ -468,5 +468,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:
@@ -784,4 +784,6 @@
     n1 = int(args[1]) if len(args) > 1 else 1
     n2 = int(args[2]) if len(args) > 2 else 1
+    use_sasview = any(engine=='sasview' and count>0
+                      for engine, count in zip(engines, [n1, n2]))
 
     # Get demo parameters from model definition, or use default parameters
@@ -811,5 +813,6 @@
     #import pprint; pprint.pprint(model_info)
     constrain_pars(model_info, pars)
-    constrain_new_to_old(model_info, pars)
+    if use_sasview:
+        constrain_new_to_old(model_info, pars)
     if opts['show_pars']:
         print(str(parlist(model_info, pars, opts['is2d'])))
Index: sasmodels/convert.py
===================================================================
--- sasmodels/convert.py	(revision ce896fde473387374a3fcfc3918ce7a13676b16c)
+++ sasmodels/convert.py	(revision 2a3e1f507e690f45c1529e7d23d9f8fccf7e5008)
@@ -138,5 +138,7 @@
     namelist = name.split('*') if '*' in name else [name]
     for name in namelist:
-        if name == 'pearl_necklace':
+        if name == 'stacked_disks':
+            _remove_pd(oldpars, 'n_stacking', name)
+        elif name == 'pearl_necklace':
             _remove_pd(oldpars, 'num_pearls', name)
             _remove_pd(oldpars, 'thick_string', name)
Index: sasmodels/data.py
===================================================================
--- sasmodels/data.py	(revision 092cb3ce952732db3425891b8c1619f81fd7a157)
+++ sasmodels/data.py	(revision d6f5da6a131a8de1b3e3d9a12af714a8d32aa4be)
@@ -41,5 +41,5 @@
     Load data using a sasview loader.
     """
-    from sas.dataloader.loader import Loader
+    from sas.sascalc.dataloader.loader import Loader
     loader = Loader()
     data = loader.load(filename)
@@ -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,9 +378,14 @@
     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
 
     if use_data or use_theory:
+        if num_plots > 1:
+            plt.subplot(1, num_plots, 1)
+
         #print(vmin, vmax)
         all_positive = True
@@ -399,5 +406,5 @@
             if view is 'log':
                 mtheory[mtheory <= 0] = masked
-            plt.plot(data.x/10, scale*mtheory, '-', hold=True)
+            plt.plot(data.x, scale*mtheory, '-', hold=True)
             all_positive = all_positive and (mtheory > 0).all()
             some_present = some_present or (mtheory.count() > 0)
@@ -406,6 +413,4 @@
             plt.ylim(*limits)
 
-        if num_plots > 1:
-            plt.subplot(1, num_plots, 1)
         plt.xscale('linear' if not some_present else view)
         plt.yscale('linear'
@@ -415,4 +420,15 @@
         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[qy>0], np.log10(Iqxy[qy>0,:]))
+        plt.xlabel("$q_x$/A$^{-1}$")
+        plt.xlabel("$q_y$/A$^{-1}$")
+        plt.xscale('log')
+        plt.yscale('log')
+        #plt.axis('equal')
+
     if use_resid:
         mresid = masked_array(resid, data.mask.copy())
@@ -421,5 +437,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}$")
@@ -561,13 +577,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 60eab2a7f2ffda1a9f4c3f3767712387c0e915cc)
+++ sasmodels/direct_model.py	(revision ea7504385b3e91f10b68265b8e6a5a5519f04260)
@@ -62,4 +62,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'
@@ -113,6 +115,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])
@@ -120,4 +122,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:
@@ -139,5 +158,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
@@ -152,13 +171,13 @@
         if self._kernel is None:
             self._kernel = self._model.make_kernel(self._kernel_inputs)
-            self._kernel_mono = (
-                self._model.make_kernel(self._kernel_mono_inputs)
-                if self._kernel_mono_inputs else None
-            )
+            self._kernel_mono = (self._model.make_kernel(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, 
@@ -166,4 +185,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/models/_spherepy.py
===================================================================
--- sasmodels/models/_spherepy.py	(revision 49da0794022419abc59514fd36d5fed434d4075b)
+++ sasmodels/models/_spherepy.py	(revision d7028dc22ec2cdba0543b07d743e678db43f44f1)
@@ -1,6 +1,5 @@
 r"""
-For information about polarised and magnetic scattering, click here_.
-
-.. _here: polar_mag_help.html
+For information about polarised and magnetic scattering, see 
+the :doc:`magnetic help <../sasgui/perspectives/fitting/mag_help>` documentation.
 
 Definition
Index: sasmodels/models/fuzzy_sphere.py
===================================================================
--- sasmodels/models/fuzzy_sphere.py	(revision 0cc31e163b64eca925abbec18fe8338cc1fc7450)
+++ sasmodels/models/fuzzy_sphere.py	(revision d7028dc22ec2cdba0543b07d743e678db43f44f1)
@@ -1,6 +1,5 @@
 r"""
-For information about polarised and magnetic scattering, click here_.
-
-.. _here: polar_mag_help.html
+For information about polarised and magnetic scattering, see 
+the :doc:`magnetic help <../sasgui/perspectives/fitting/mag_help>` documentation.
 
 Definition
Index: sasmodels/models/multilayer_vesicle.py
===================================================================
--- sasmodels/models/multilayer_vesicle.py	(revision c6ca41e077aa705fcd195f07f143a6daaf729840)
+++ sasmodels/models/multilayer_vesicle.py	(revision d7028dc22ec2cdba0543b07d743e678db43f44f1)
@@ -29,9 +29,6 @@
     is used as the effective radius for *S(Q)* when $P(Q) * S(Q)$ is applied.
 
-
-
-For information about polarised and magnetic scattering, click here_.
-
-.. _here: polar_mag_help.html
+For information about polarised and magnetic scattering, see 
+the :doc:`magnetic help <../sasgui/perspectives/fitting/mag_help>` documentation.
 
 This code is based on the form factor calculations implemented in the NIST
Index: sasmodels/models/sphere.py
===================================================================
--- sasmodels/models/sphere.py	(revision 49da0794022419abc59514fd36d5fed434d4075b)
+++ sasmodels/models/sphere.py	(revision d7028dc22ec2cdba0543b07d743e678db43f44f1)
@@ -1,6 +1,5 @@
 r"""
-For information about polarised and magnetic scattering, click here_.
-
-.. _here: polar_mag_help.html
+For information about polarised and magnetic scattering, see 
+the :doc:`magnetic help <../sasgui/perspectives/fitting/mag_help>` documentation.
 
 Definition
Index: sasmodels/models/stacked_disks.py
===================================================================
--- sasmodels/models/stacked_disks.py	(revision e664a1142e1d78d069b30f4e992cf69b4984ce9a)
+++ sasmodels/models/stacked_disks.py	(revision 53215cfe7c83e780da26e7549f7f4aa8555e8261)
@@ -100,4 +100,9 @@
 J S Higgins and H C Benoit, *Polymers and Neutron Scattering*, Clarendon, Oxford, 1994
 
+**Author:** NIST IGOR/DANSE **on:** pre 2010
+
+**Last Modified by:** Piotr Rozyczko **on:** February 18, 2016
+
+**Last Reviewed by:** Richard Heenan **on:** March 22, 2016
 """
 
@@ -157,5 +162,6 @@
 
 tests = [
-    # Accuracy tests based on content in test/utest_extra_models.py
+    # Accuracy tests based on content in test/utest_extra_models.py.
+    # Added 2 tests with n_stacked = 5 using SasView 3.1.2 - PDB
     [{'core_thick': 10.0,
       'layer_thick': 15.0,
@@ -171,4 +177,32 @@
       'background': 0.001,
      }, 0.001, 5075.12],
+
+    [{'core_thick': 10.0,
+      'layer_thick': 15.0,
+      'radius': 3000.0,
+      'n_stacking': 5.0,
+      'sigma_d': 0.0,
+      'sld_core': 4.0,
+      'sld_layer': -0.4,
+      'solvent_sd': 5.0,
+      'theta': 0.0,
+      'phi': 0.0,
+      'scale': 0.01,
+      'background': 0.001,
+     }, 0.001, 5065.12793824],
+
+    [{'core_thick': 10.0,
+      'layer_thick': 15.0,
+      'radius': 3000.0,
+      'n_stacking': 5.0,
+      'sigma_d': 0.0,
+      'sld_core': 4.0,
+      'sld_layer': -0.4,
+      'solvent_sd': 5.0,
+      'theta': 0.0,
+      'phi': 0.0,
+      'scale': 0.01,
+      'background': 0.001,
+     }, 0.164, 0.0127673597265],
 
     [{'core_thick': 10.0,
Index: sasmodels/resolution.py
===================================================================
--- sasmodels/resolution.py	(revision 48fbd50a25db17a9c6bf7770aeb061514b3e8879)
+++ 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:]
@@ -596,5 +596,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)
@@ -606,5 +606,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)
@@ -621,5 +621,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)
@@ -635,5 +635,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)
@@ -648,5 +648,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)
@@ -773,5 +773,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.
@@ -793,5 +793,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
@@ -811,5 +811,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))
@@ -1068,5 +1068,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 d6f5da6a131a8de1b3e3d9a12af714a8d32aa4be)
@@ -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':1000, 'high':500, 'med':200, 'low':50}
+N_SLIT_PERP_DOC = ", ".join("%s=%d"%(name,value) for value,name in
+                            sorted((2*v+1,k) for k,v in N_SLIT_PERP.items()))
 
 class Pinhole2D(Resolution):
@@ -167,2 +173,65 @@
         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_min, qy_max = np.log10(np.min(q)), np.log10(qy_width)
+        qy_calc = np.logspace(qy_min, qy_max, N_SLIT_PERP[accuracy])
+        qy_calc = np.hstack((-qy_calc[::-1], 0, qy_calc))
+        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.trapz(theory.reshape(self.ny, self.nx), axis=0, x=self.qy_calc)
+        if self.weights is not None:
+            Iq = resolution.apply_resolution_matrix(self.weights, Iq)
+        return Iq
+
