Index: extra/pylint.rc
===================================================================
--- extra/pylint.rc	(revision d4666ca59132213b621b9cd27fda96521894152e)
+++ extra/pylint.rc	(revision 823e6204f35d6c2bf32f27a3a3c137a369a973cf)
@@ -10,5 +10,5 @@
 
 # Profiled execution.
-profile=no
+#profile=no
 
 # Add files or directories to the blacklist. They should be base names, not
@@ -92,5 +92,5 @@
 # Add a comment according to your evaluation note. This is used by the global
 # evaluation report (RP0004).
-comment=no
+#comment=no
 
 # Template used to display messages. This is a python new-style format string
@@ -101,5 +101,5 @@
 
 # Required attributes for module, separated by a comma
-required-attributes=
+#required-attributes=
 
 # List of builtins function names that should not be used, separated by a comma
@@ -108,10 +108,9 @@
 # Good variable names which should always be accepted, separated by a comma
 good-names=_,
-    input,
-    h,i,j,k,n,w,x,y,z,
-    q,qx,qy,qz,
-    dt,dx,dy,dz,id,
-    Iq,dIq,Qx,Qy,Qz,Iqxy,
-    p,
+    input, id,
+    h, i, j, k, n, p, v, w, x, y, z,
+    q, qx, qy, qz, Q, Qx, Qy, Qz,
+    dt, dx, dy, dz, dq,
+    Iq, dIq, Iqxy, Iq_calc,
     ER, call_ER, VR, call_VR,
     Rmax, SElength,
@@ -281,5 +280,6 @@
 # (useful for modules/projects where namespaces are manipulated during runtime
 # and thus existing member attributes cannot be deduced by static analysis
-ignored-modules=numpy,np
+ignored-modules=numpy,np,numpy.random,
+    bumps,sas,
 
 # List of classes names for which member attributes should not be checked
@@ -289,5 +289,5 @@
 # When zope mode is activated, add a predefined set of Zope acquired attributes
 # to generated-members.
-zope=no
+#zope=no
 
 # List of members which are set dynamically and missed by pylint inference
@@ -319,5 +319,5 @@
 # List of interface methods to ignore, separated by a comma. This is used for
 # instance to not check methods defines in Zope's Interface base class.
-ignore-iface-methods=isImplementedBy,deferred,extends,names,namesAndDescriptions,queryDescriptionFor,getBases,getDescriptionFor,getDoc,getName,getTaggedValue,getTaggedValueTags,isEqualOrExtendedBy,setTaggedValue,isImplementedByInstancesOf,adaptWith,is_implemented_by
+#ignore-iface-methods=isImplementedBy,deferred,extends,names,namesAndDescriptions,queryDescriptionFor,getBases,getDescriptionFor,getDoc,getName,getTaggedValue,getTaggedValueTags,isEqualOrExtendedBy,setTaggedValue,isImplementedByInstancesOf,adaptWith,is_implemented_by
 
 # List of method names used to declare (i.e. assign) instance attributes.
Index: extra/pylint_numpy.py
===================================================================
--- extra/pylint_numpy.py	(revision 63b32bbc85e7dea7a022cf604ac020dee56d7a88)
+++ extra/pylint_numpy.py	(revision 823e6204f35d6c2bf32f27a3a3c137a369a973cf)
@@ -8,5 +8,8 @@
     #print("processing",module.name)
     if module.name.startswith('numpy'):
-        if module.name == 'numpy':  import numpy
-        elif module.name == 'numpy.random': import numpy.random
+        if module.name == 'numpy':
+            import numpy
+        elif module.name == 'numpy.random':
+            import numpy.random
+            from numpy.random import seed, get_state, set_state
 
Index: extra/pylint_pyopencl.py
===================================================================
--- extra/pylint_pyopencl.py	(revision 3c56da87f0995bdde77f4183b68e77b81584dde0)
+++ extra/pylint_pyopencl.py	(revision 823e6204f35d6c2bf32f27a3a3c137a369a973cf)
@@ -9,2 +9,3 @@
     if module.name == 'pyopencl':
         import pyopencl
+        import pyopencl as cl
Index: extra/run-pylint.py
===================================================================
--- extra/run-pylint.py	(revision f01e53d5a96b16346f85048d4f125bc41f154611)
+++ extra/run-pylint.py	(revision 823e6204f35d6c2bf32f27a3a3c137a369a973cf)
@@ -6,11 +6,21 @@
 
 def main():
+    extra = abspath(dirname(__file__))
+    root = abspath(joinpath(extra, '..'))
+
     envpath = os.environ.get('PYTHONPATH',None)
     path = [envpath] if envpath else []
-    path.append(abspath(dirname(__file__))) # so we can find the plugins
+    path.append(extra)
+
+    #bumps = abspath(joinpath(root, "..", "bumps"))
+    #periodictable = abspath(joinpath(root, "..", "periodictable"))
+    #sasview = abspath(joinpath(root, "..", "sasview", "src"))
+    #path.extend((bumps, periodictable, sasview))
+
     os.environ['PYTHONPATH'] = ':'.join(path)
-    root = abspath(joinpath(dirname(__file__), '..'))
+
+    # Run the lint command
+    cmd = "pylint --rcfile extra/pylint.rc -f parseable sasmodels"
     os.chdir(root)
-    cmd = "pylint --rcfile extra/pylint.rc -f parseable sasmodels"
     status = os.system(cmd)
     sys.exit(status)
Index: sasmodels/compare.py
===================================================================
--- sasmodels/compare.py	(revision 190fc2b167d3d89ca9c3b7d3128bdf8bde9f67d0)
+++ sasmodels/compare.py	(revision 190fc2b167d3d89ca9c3b7d3128bdf8bde9f67d0)
@@ -604,5 +604,7 @@
 def columnize(L, indent="", width=79):
     """
-    Format a list of strings into columns for printing.
+    Format a list of strings into columns.
+
+    Returns a string with carriage returns ready for printing.
     """
     column_width = max(len(w) for w in L) + 1
Index: sasmodels/core.py
===================================================================
--- sasmodels/core.py	(revision 190fc2b167d3d89ca9c3b7d3128bdf8bde9f67d0)
+++ sasmodels/core.py	(revision eafc9fa8b200a3b754e7dafb960c8891361eacef)
@@ -38,4 +38,7 @@
     """
     Load a model definition given the model name.
+
+    This returns a handle to the module defining the model.  This can be
+    used with functions in generate to build the docs or extract model info.
     """
     __import__('sasmodels.models.'+model_name)
@@ -115,6 +118,5 @@
     Return a computation kernel from the model definition and the q input.
     """
-    model_input = model.make_input(q_vectors)
-    return model(model_input)
+    return model(q_vectors)
 
 def get_weights(info, pars, name):
Index: sasmodels/data.py
===================================================================
--- sasmodels/data.py	(revision 69ec80f641366fc1c80fed9bcfff0d283905e6ef)
+++ sasmodels/data.py	(revision 5c962dff13d0a0a2f6d779be8df775f540958dbc)
@@ -87,4 +87,30 @@
 
 class Data1D(object):
+    """
+    1D data object.
+
+    Note that this definition matches the attributes from sasview, with
+    some generic 1D data vectors and some SAS specific definitions.  Some
+    refactoring to allow consistent naming conventions between 1D, 2D and
+    SESANS data would be helpful.
+
+    **Attributes**
+
+    *x*, *dx*: $q$ vector and gaussian resolution
+
+    *y*, *dy*: $I(q)$ vector and measurement uncertainty
+
+    *mask*: values to include in plotting/analysis
+
+    *dxl*: slit widths for slit smeared data, with *dx* ignored
+
+    *qmin*, *qmax*: range of $q$ values in *x*
+
+    *filename*: label for the data line
+
+    *_xaxis*, *_xunit*: label and units for the *x* axis
+
+    *_yaxis*, *_yunit*: label and units for the *y* axis
+    """
     def __init__(self, x=None, y=None, dx=None, dy=None):
         self.x, self.y, self.dx, self.dy = x, y, dx, dy
@@ -93,5 +119,8 @@
         self.qmin = x.min() if x is not None else np.NaN
         self.qmax = x.max() if x is not None else np.NaN
-        self.mask = np.isnan(y) if y is not None else None
+        # TODO: why is 1D mask False and 2D mask True?
+        self.mask = (np.isnan(y) if y is not None
+                     else np.zeros_like(x, 'b') if x is not None
+                     else None)
         self._xaxis, self._xunit = "x", ""
         self._yaxis, self._yunit = "y", ""
@@ -114,9 +143,42 @@
 
 class Data2D(object):
+    """
+    2D data object.
+
+    Note that this definition matches the attributes from sasview. Some
+    refactoring to allow consistent naming conventions between 1D, 2D and
+    SESANS data would be helpful.
+
+    **Attributes**
+
+    *qx_data*, *dqx_data*: $q_x$ matrix and gaussian resolution
+
+    *qy_data*, *dqy_data*: $q_y$ matrix and gaussian resolution
+
+    *data*, *err_data*: $I(q)$ matrix and measurement uncertainty
+
+    *mask*: values to exclude from plotting/analysis
+
+    *qmin*, *qmax*: range of $q$ values in *x*
+
+    *filename*: label for the data line
+
+    *_xaxis*, *_xunit*: label and units for the *x* axis
+
+    *_yaxis*, *_yunit*: label and units for the *y* axis
+
+    *_zaxis*, *_zunit*: label and units for the *y* axis
+
+    *Q_unit*, *I_unit*: units for Q and intensity
+
+    *x_bins*, *y_bins*: grid steps in *x* and *y* directions
+    """
     def __init__(self, x=None, y=None, z=None, dx=None, dy=None, dz=None):
         self.qx_data, self.dqx_data = x, dx
         self.qy_data, self.dqy_data = y, dy
         self.data, self.err_data = z, dz
-        self.mask = ~np.isnan(z) if z is not None else None
+        self.mask = (~np.isnan(z) if z is not None
+                     else np.ones_like(x) if x is not None
+                     else None)
         self.q_data = np.sqrt(x**2 + y**2)
         self.qmin = 1e-16
@@ -126,7 +188,7 @@
         self.Q_unit = "1/A"
         self.I_unit = "1/cm"
-        self.xaxis("Q_x", "A^{-1}")
-        self.yaxis("Q_y", "A^{-1}")
-        self.zaxis("Intensity", r"\text{cm}^{-1}")
+        self.xaxis("Q_x", "1/A")
+        self.yaxis("Q_y", "1/A")
+        self.zaxis("Intensity", "1/cm")
         self._xaxis, self._xunit = "x", ""
         self._yaxis, self._yunit = "y", ""
@@ -157,4 +219,7 @@
 
 class Vector(object):
+    """
+    3-space vector of *x*, *y*, *z*
+    """
     def __init__(self, x=None, y=None, z=None):
         self.x, self.y, self.z = x, y, z
@@ -235,4 +300,11 @@
     """
     Plot data loaded by the sasview loader.
+
+    *data* is a sasview data object, either 1D, 2D or SESANS.
+
+    *view* is log or linear.
+
+    *limits* sets the intensity limits on the plot; if None then the limits
+    are inferred from the data.
     """
     # Note: kind of weird using the plot result functions to plot just the
@@ -249,4 +321,19 @@
 def plot_theory(data, theory, resid=None, view='log',
                 use_data=True, limits=None):
+    """
+    Plot theory calculation.
+
+    *data* is needed to define the graph properties such as labels and
+    units, and to define the data mask.
+
+    *theory* is a matrix of the same shape as the data.
+
+    *view* is log or linear
+
+    *use_data* is True if the data should be plotted as well as the theory.
+
+    *limits* sets the intensity limits on the plot; if None then the limits
+    are inferred from the data.
+    """
     if hasattr(data, 'lam'):
         _plot_result_sesans(data, theory, resid, use_data=True, limits=limits)
@@ -258,7 +345,16 @@
 
 def protect(fn):
+    """
+    Decorator to wrap calls in an exception trapper which prints the
+    exception and continues.  Keyboard interrupts are ignored.
+    """
     def wrapper(*args, **kw):
+        """
+        Trap and print errors from function.
+        """
         try:
             return fn(*args, **kw)
+        except KeyboardInterrupt:
+            raise
         except:
             traceback.print_exc()
@@ -332,4 +428,7 @@
 @protect
 def _plot_result_sesans(data, theory, resid, use_data, limits=None):
+    """
+    Plot SESANS results.
+    """
     import matplotlib.pyplot as plt
     use_data = use_data and data.y is not None
@@ -459,4 +558,7 @@
 
 def demo():
+    """
+    Load and plot a SAS dataset.
+    """
     data = load_data('DEC07086.DAT')
     set_beam_stop(data, 0.004)
Index: sasmodels/direct_model.py
===================================================================
--- sasmodels/direct_model.py	(revision 9404dd3f599fcbe72790aae8a6d176e64f292ea6)
+++ sasmodels/direct_model.py	(revision eafc9fa8b200a3b754e7dafb960c8891361eacef)
@@ -1,3 +1,25 @@
-import warnings
+"""
+Class interface to the model calculator.
+
+Calling a model is somewhat non-trivial since the functions called depend
+on the data type.  For 1D data the *Iq* kernel needs to be called, for
+2D data the *Iqxy* kernel needs to be called, and for SESANS data the
+*Iq* kernel needs to be called followed by a Hankel transform.  Before
+the kernel is called an appropriate *q* calculation vector needs to be
+constructed.  This is not the simple *q* vector where you have measured
+the data since the resolution calculation will require values beyond the
+range of the measured data.  After the calculation the resolution calculator
+must be called to return the predicted value for each measured data point.
+
+:class:`DirectModel` is a callable object that takes *parameter=value*
+keyword arguments and returns the appropriate theory values for the data.
+
+:class:`DataMixin` does the real work of interpreting the data and calling
+the model calculator.  This is used by :class:`DirectModel`, which uses
+direct parameter values and by :class:`bumps_model.Experiment` which wraps
+the parameter values in boxes so that the user can set fitting ranges, etc.
+on the individual parameters and send the model to the Bumps optimizers.
+"""
+from __future__ import print_function
 
 import numpy as np
@@ -19,6 +41,18 @@
     currently used by :class:`sasview_model.SasviewModel` since this will
     require a number of changes to SasView before we can do it.
+
+    :meth:`_interpret_data` initializes the data structures necessary
+    to manage the calculations.  This sets attributes in the child class
+    such as *data_type* and *resolution*.
+
+    :meth:`_calc_theory` evaluates the model at the given control values.
+
+    :meth:`_set_data` sets the intensity data in the data object,
+    possibly with random noise added.  This is useful for simulating a
+    dataset with the results from :meth:`_calc_theory`.
     """
     def _interpret_data(self, data, model):
+        # pylint: disable=attribute-defined-outside-init
+
         self._data = data
         self._model = model
@@ -36,8 +70,10 @@
         if self.data_type == 'sesans':
             q = sesans.make_q(data.sample.zacceptance, data.Rmax)
-            self.index = slice(None, None)
+            index = slice(None, None)
+            res = None
             if data.y is not None:
-                self.Iq = data.y
-                self.dIq = data.dy
+                Iq, dIq = data.y, data.dy
+            else:
+                Iq, dIq = None, None
             #self._theory = np.zeros_like(q)
             q_vectors = [q]
@@ -49,35 +85,39 @@
             qmax = getattr(data, 'qmax', np.inf)
             accuracy = getattr(data, 'accuracy', 'Low')
-            self.index = ~data.mask & (q >= qmin) & (q <= qmax)
+            index = ~data.mask & (q >= qmin) & (q <= qmax)
             if data.data is not None:
-                self.index &= ~np.isnan(data.data)
-                self.Iq = data.data[self.index]
-                self.dIq = data.err_data[self.index]
-            self.resolution = resolution2d.Pinhole2D(data=data, index=self.index,
-                                                     nsigma=3.0, accuracy=accuracy)
+                index &= ~np.isnan(data.data)
+                Iq = data.data[index]
+                dIq = data.err_data[index]
+            else:
+                Iq, dIq = None, None
+            res = resolution2d.Pinhole2D(data=data, index=index,
+                                         nsigma=3.0, accuracy=accuracy)
             #self._theory = np.zeros_like(self.Iq)
-            q_vectors = self.resolution.q_calc
+            q_vectors = res.q_calc
         elif self.data_type == 'Iq':
-            self.index = (data.x >= data.qmin) & (data.x <= data.qmax)
+            index = (data.x >= data.qmin) & (data.x <= data.qmax)
             if data.y is not None:
-                self.index &= ~np.isnan(data.y)
-                self.Iq = data.y[self.index]
-                self.dIq = data.dy[self.index]
+                index &= ~np.isnan(data.y)
+                Iq = data.y[index]
+                dIq = data.dy[index]
+            else:
+                Iq, dIq = None, None
             if getattr(data, 'dx', None) is not None:
-                q, dq = data.x[self.index], data.dx[self.index]
-                if (dq>0).any():
-                    self.resolution = resolution.Pinhole1D(q, dq)
+                q, dq = data.x[index], data.dx[index]
+                if (dq > 0).any():
+                    res = resolution.Pinhole1D(q, dq)
                 else:
-                    self.resolution = resolution.Perfect1D(q)
-            elif (getattr(data, 'dxl', None) is not None and
-                          getattr(data, 'dxw', None) is not None):
-                self.resolution = resolution.Slit1D(data.x[self.index],
-                                                    width=data.dxh[self.index],
-                                                    height=data.dxw[self.index])
-            else:
-                self.resolution = resolution.Perfect1D(data.x[self.index])
+                    res = resolution.Perfect1D(q)
+            elif (getattr(data, 'dxl', None) is not None
+                  and getattr(data, 'dxw', None) is not None):
+                res = resolution.Slit1D(data.x[index],
+                                        width=data.dxh[index],
+                                        height=data.dxw[index])
+            else:
+                res = resolution.Perfect1D(data.x[index])
 
             #self._theory = np.zeros_like(self.Iq)
-            q_vectors = [self.resolution.q_calc]
+            q_vectors = [res.q_calc]
         else:
             raise ValueError("Unknown data type") # never gets here
@@ -87,6 +127,9 @@
         self._kernel_inputs = [v for v in q_vectors]
         self._kernel = None
+        self.Iq, self.dIq, self.index = Iq, dIq, index
+        self.resolution = res
 
     def _set_data(self, Iq, noise=None):
+        # pylint: disable=attribute-defined-outside-init
         if noise is not None:
             self.dIq = Iq*noise*0.01
@@ -106,6 +149,5 @@
     def _calc_theory(self, pars, cutoff=0.0):
         if self._kernel is None:
-            q_input = self._model.make_input(self._kernel_inputs)
-            self._kernel = self._model(q_input)
+            self._kernel = make_kernel(self._model, self._kernel_inputs)  # pylint: disable=attribute-defined-outside-init
 
         Iq_calc = call_kernel(self._kernel, pars, cutoff=cutoff)
@@ -120,20 +162,50 @@
 
 class DirectModel(DataMixin):
+    """
+    Create a calculator object for a model.
+
+    *data* is 1D SAS, 2D SAS or SESANS data
+
+    *model* is a model calculator return from :func:`generate.load_model`
+
+    *cutoff* is the polydispersity weight cutoff.
+    """
     def __init__(self, data, model, cutoff=1e-5):
         self.model = model
         self.cutoff = cutoff
+        # Note: _interpret_data defines the model attributes
         self._interpret_data(data, model)
-        self.kernel = make_kernel(self.model, self._kernel_inputs)
+
     def __call__(self, **pars):
         return self._calc_theory(pars, cutoff=self.cutoff)
+
     def ER(self, **pars):
+        """
+        Compute the equivalent radius for the model.
+
+        Return 0. if not defined.
+        """
         return call_ER(self.model.info, pars)
+
     def VR(self, **pars):
+        """
+        Compute the equivalent volume for the model, including polydispersity
+        effects.
+
+        Return 1. if not defined.
+        """
         return call_VR(self.model.info, pars)
+
     def simulate_data(self, noise=None, **pars):
+        """
+        Generate simulated data for the model.
+        """
         Iq = self.__call__(**pars)
         self._set_data(Iq, noise=noise)
 
-def demo():
+def main():
+    """
+    Program to evaluate a particular model at a set of q values.
+    """
     import sys
     from .data import empty_data1D, empty_data2D
@@ -144,5 +216,5 @@
     model_name = sys.argv[1]
     call = sys.argv[2].upper()
-    if call not in ("ER","VR"):
+    if call not in ("ER", "VR"):
         try:
             values = [float(v) for v in call.split(',')]
@@ -153,6 +225,6 @@
             data = empty_data1D([q])
         elif len(values) == 2:
-            qx,qy = values
-            data = empty_data2D([qx],[qy])
+            qx, qy = values
+            data = empty_data2D([qx], [qy])
         else:
             print("use q or qx,qy or ER or VR")
@@ -164,7 +236,7 @@
     model = load_model(model_definition, dtype='single')
     calculator = DirectModel(data, model)
-    pars = dict((k,float(v))
+    pars = dict((k, float(v))
                 for pair in sys.argv[3:]
-                for k,v in [pair.split('=')])
+                for k, v in [pair.split('=')])
     if call == "ER":
         print(calculator.ER(**pars))
@@ -176,3 +248,3 @@
 
 if __name__ == "__main__":
-    demo()
+    main()
Index: sasmodels/exception.py
===================================================================
--- sasmodels/exception.py	(revision eaca9ebe6d0022b8ac658c30e21ef334fb4f3861)
+++ sasmodels/exception.py	(revision 823e6204f35d6c2bf32f27a3a3c137a369a973cf)
@@ -8,4 +8,7 @@
 except NameError:
     class WindowsError(Exception):
+        """
+        Fake WindowsException when not on Windows.
+        """
         pass
 
@@ -31,5 +34,5 @@
     # TODO: try to incorporate current stack trace in the raised exception
     if isinstance(exc, WindowsError):
-        raise OSError(str(exc)+" "+msg)
+        raise OSError(str(exc) + " " + msg)
 
     args = exc.args
@@ -38,6 +41,6 @@
     else:
         try:
-            arg0 = " ".join((args[0],msg))
+            arg0 = " ".join((args[0], msg))
             exc.args = tuple([arg0] + list(args[1:]))
         except:
-            exc.args = (" ".join((str(exc),msg)),)
+            exc.args = (" ".join((str(exc), msg)),)
Index: sasmodels/generate.py
===================================================================
--- sasmodels/generate.py	(revision 190fc2b167d3d89ca9c3b7d3128bdf8bde9f67d0)
+++ sasmodels/generate.py	(revision eafc9fa8b200a3b754e7dafb960c8891361eacef)
@@ -7,5 +7,5 @@
     particular dimensions averaged over all orientations.
 
-    *Iqxy(qx, qy, p1, p2, ...)* returns the scattering at qx,qy for a form
+    *Iqxy(qx, qy, p1, p2, ...)* returns the scattering at qx, qy for a form
     with particular dimensions for a single orientation.
 
@@ -47,5 +47,5 @@
 the *sin* and *cos* values are needed for a particular argument.  Since
 this function does not exist in C99, all use of *sincos* should be
-replaced by the macro *SINCOS(value,sn,cn)* where *sn* and *cn* are
+replaced by the macro *SINCOS(value, sn, cn)* where *sn* and *cn* are
 previously declared *double* variables.  When compiled for systems without
 OpenCL, *SINCOS* will be replaced by *sin* and *cos* calls.   If *value* is
@@ -194,4 +194,5 @@
 returns a list of files required by the model.
 """
+from __future__ import print_function
 
 # TODO: identify model files which have changed since loading and reload them.
@@ -312,5 +313,5 @@
     raise ValueError("%r not found in %s" % (filename, search_path))
 
-def sources(info):
+def model_sources(info):
     """
     Return a list of the sources file paths for the module.
@@ -372,9 +373,9 @@
 
 
-def kernel_name(info, is_2D):
+def kernel_name(info, is_2d):
     """
     Name of the exported kernel symbol.
     """
-    return info['name'] + "_" + ("Iqxy" if is_2D else "Iq")
+    return info['name'] + "_" + ("Iqxy" if is_2d else "Iq")
 
 
@@ -419,15 +420,15 @@
 
 
-def build_polydispersity_loops(pd_pars):
-    """
-    Build polydispersity loops
-
-    Returns loop opening and loop closing
-    """
-    LOOP_OPEN = """\
+LOOP_OPEN = """\
 for (int %(name)s_i=0; %(name)s_i < N%(name)s; %(name)s_i++) {
   const double %(name)s = loops[2*(%(name)s_i%(offset)s)];
   const double %(name)s_w = loops[2*(%(name)s_i%(offset)s)+1];\
 """
+def build_polydispersity_loops(pd_pars):
+    """
+    Build polydispersity loops
+
+    Returns loop opening and loop closing
+    """
     depth = 4
     offset = ""
@@ -464,5 +465,5 @@
 
     # Load additional sources
-    source = [open(f).read() for f in sources(info)]
+    source = [open(f).read() for f in model_sources(info)]
 
     # Prepare defines
@@ -572,9 +573,9 @@
 
     #for d in defines: print(d)
-    DEFINES = '\n'.join('#define %s %s' % (k, v) for k, v in defines)
-    SOURCES = '\n\n'.join(source)
+    defines = '\n'.join('#define %s %s' % (k, v) for k, v in defines)
+    sources = '\n\n'.join(source)
     return C_KERNEL_TEMPLATE % {
-        'DEFINES':DEFINES,
-        'SOURCES':SOURCES,
+        'DEFINES': defines,
+        'SOURCES': sources,
         }
 
@@ -590,5 +591,5 @@
     demo_parameters = getattr(kernel_module, 'demo', None)
     if demo_parameters is None:
-        demo_parameters = dict((p[0],p[2]) for p in parameters)
+        demo_parameters = dict((p[0], p[2]) for p in parameters)
     filename = abspath(kernel_module.__file__)
     kernel_id = splitext(basename(filename))[0]
@@ -597,5 +598,5 @@
         name = " ".join(w.capitalize() for w in kernel_id.split('_'))
     info = dict(
-        id = kernel_id,  # string used to load the kernel
+        id=kernel_id,  # string used to load the kernel
         filename=abspath(kernel_module.__file__),
         name=name,
@@ -643,5 +644,5 @@
         elif section_marker.match(line):
             if len(line) >= len(prior):
-                yield "".join( ("**",prior,"**") )
+                yield "".join(("**", prior, "**"))
                 prior = None
             else:
Index: sasmodels/kernelcl.py
===================================================================
--- sasmodels/kernelcl.py	(revision cde11f0f857d1861f37ded681777699e3aa91738)
+++ sasmodels/kernelcl.py	(revision eafc9fa8b200a3b754e7dafb960c8891361eacef)
@@ -1,4 +1,4 @@
 """
-GPU support through OpenCL
+GPU driver for C kernels
 
 There should be a single GPU environment running on the system.  This
@@ -152,10 +152,4 @@
 
 
-def make_result(self, size):
-    self.res = np.empty(size, dtype=self.dtype)
-    self.res_b = cl.Buffer(self.program.context, mf.READ_WRITE, self.res.nbytes)
-    return self.res, self.res_b
-
-
 # for now, this returns one device in the context
 # TODO: create a context that contains all devices on all platforms
@@ -183,8 +177,16 @@
 
     def has_type(self, dtype):
+        """
+        Return True if all devices support a given type.
+        """
         dtype = generate.F32 if dtype == 'fast' else np.dtype(dtype)
         return all(has_type(d, dtype) for d in self.context.devices)
 
     def _create_some_context(self):
+        """
+        Protected call to cl.create_some_context without interactivity.  Use
+        this if PYOPENCL_CTX is set in the environment.  Sets the *context*
+        attribute.
+        """
         try:
             self.context = cl.create_some_context(interactive=False)
@@ -195,4 +197,7 @@
 
     def compile_program(self, name, source, dtype, fast=False):
+        """
+        Compile the program for the device in the given context.
+        """
         key = "%s-%s-%s"%(name, dtype, fast)
         if key not in self.compiled:
@@ -204,4 +209,7 @@
 
     def release_program(self, name):
+        """
+        Free memory associated with the program on the device.
+        """
         if name in self.compiled:
             self.compiled[name].release()
@@ -209,4 +217,7 @@
 
 def _get_default_context():
+    """
+    Get an OpenCL context, preferring GPU over CPU.
+    """
     default = None
     for platform in cl.get_platforms():
@@ -241,42 +252,35 @@
         self.info = info
         self.source = source
-        self.dtype = generate.F32 if dtype=='fast' else np.dtype(dtype)
+        self.dtype = generate.F32 if dtype == 'fast' else np.dtype(dtype)
         self.fast = (dtype == 'fast')
         self.program = None # delay program creation
 
     def __getstate__(self):
-        state = self.__dict__.copy()
-        state['program'] = None
-        return state
+        return self.info, self.source, self.dtype, self.fast
 
     def __setstate__(self, state):
-        self.__dict__ = state.copy()
-
-    def __call__(self, q_input):
-        if self.dtype != q_input.dtype:
-            raise TypeError("data is %s kernel is %s"
-                            % (q_input.dtype, self.dtype))
+        self.info, self.source, self.dtype, self.fast = state
+        self.program = None
+
+    def __call__(self, q_vectors):
         if self.program is None:
             compiler = environment().compile_program
             self.program = compiler(self.info['name'], self.source, self.dtype,
                                     self.fast)
-        kernel_name = generate.kernel_name(self.info, q_input.is_2D)
+        is_2d = len(q_vectors) == 2
+        kernel_name = generate.kernel_name(self.info, is_2d)
         kernel = getattr(self.program, kernel_name)
-        return GpuKernel(kernel, self.info, q_input)
+        return GpuKernel(kernel, self.info, q_vectors, self.dtype)
 
     def release(self):
+        """
+        Free the resources associated with the model.
+        """
         if self.program is not None:
             environment().release_program(self.info['name'])
             self.program = None
 
-    def make_input(self, q_vectors):
-        """
-        Make q input vectors available to the model.
-
-        Note that each model needs its own q vector even if the case of
-        mixture models because some models may be OpenCL, some may be
-        ctypes and some may be pure python.
-        """
-        return GpuInput(q_vectors, dtype=self.dtype)
+    def __del__(self):
+        self.release()
 
 # TODO: check that we don't need a destructor for buffers which go out of scope
@@ -304,5 +308,5 @@
         self.nq = q_vectors[0].size
         self.dtype = np.dtype(dtype)
-        self.is_2D = (len(q_vectors) == 2)
+        self.is_2d = (len(q_vectors) == 2)
         # TODO: stretch input based on get_warp()
         # not doing it now since warp depends on kernel, which is not known
@@ -317,18 +321,25 @@
 
     def release(self):
+        """
+        Free the memory.
+        """
         for b in self.q_buffers:
             b.release()
         self.q_buffers = []
 
+    def __del__(self):
+        self.release()
+
 class GpuKernel(object):
     """
     Callable SAS kernel.
 
-    *kernel* is the GpuKernel object to call.
+    *kernel* is the GpuKernel object to call
 
     *info* is the module information
 
-    *q_input* is the DllInput q vectors at which the kernel should be
-    evaluated.
+    *q_vectors* is the q vectors at which the kernel should be evaluated
+
+    *dtype* is the kernel precision
 
     The resulting call method takes the *pars*, a list of values for
@@ -340,10 +351,10 @@
     Call :meth:`release` when done with the kernel instance.
     """
-    def __init__(self, kernel, info, q_input):
-        self.q_input = q_input
+    def __init__(self, kernel, info, q_vectors, dtype):
+        q_input = GpuInput(q_vectors, dtype)
         self.kernel = kernel
         self.info = info
         self.res = np.empty(q_input.nq, q_input.dtype)
-        dim = '2d' if q_input.is_2D else '1d'
+        dim = '2d' if q_input.is_2d else '1d'
         self.fixed_pars = info['partype']['fixed-' + dim]
         self.pd_pars = info['partype']['pd-' + dim]
@@ -358,5 +369,5 @@
                                 q_input.global_size[0] * q_input.dtype.itemsize)
                       for _ in env.queues]
-
+        self.q_input = q_input
 
     def __call__(self, fixed_pars, pd_pars, cutoff=1e-5):
@@ -399,4 +410,7 @@
 
     def release(self):
+        """
+        Release resources associated with the kernel.
+        """
         for b in self.loops_b:
             b.release()
@@ -405,4 +419,5 @@
             b.release()
         self.res_b = []
+        self.q_input.release()
 
     def __del__(self):
Index: sasmodels/kerneldll.py
===================================================================
--- sasmodels/kerneldll.py	(revision 74667d322dd33e8f42805e77bc688d1a9a99b4fa)
+++ sasmodels/kerneldll.py	(revision eafc9fa8b200a3b754e7dafb960c8891361eacef)
@@ -1,4 +1,4 @@
 r"""
-C types wrapper for sasview models.
+DLL driver for C kernels
 
 The global attribute *ALLOW_SINGLE_PRECISION_DLLS* should be set to *True* if
@@ -44,4 +44,5 @@
 directory for this application, then OpenMP should be supported.
 """
+from __future__ import print_function
 
 import sys
@@ -122,5 +123,5 @@
     models are allowed as DLLs.
     """
-    if callable(info.get('Iq',None)):
+    if callable(info.get('Iq', None)):
         return PyModel(info)
 
@@ -139,11 +140,11 @@
 
     source = generate.convert_type(source, dtype)
-    source_files = generate.sources(info) + [info['filename']]
-    dll= dll_path(info, dtype)
+    source_files = generate.model_sources(info) + [info['filename']]
+    dll = dll_path(info, dtype)
     newest = max(os.path.getmtime(f) for f in source_files)
-    if not os.path.exists(dll) or os.path.getmtime(dll)<newest:
+    if not os.path.exists(dll) or os.path.getmtime(dll) < newest:
         # Replace with a proper temp file
-        fid, filename = tempfile.mkstemp(suffix=".c",prefix=tempfile_prefix)
-        os.fdopen(fid,"w").write(source)
+        fid, filename = tempfile.mkstemp(suffix=".c", prefix=tempfile_prefix)
+        os.fdopen(fid, "w").write(source)
         command = COMPILE%{"source":filename, "output":dll}
         print("Compile command: "+command)
@@ -160,5 +161,5 @@
 def load_dll(source, info, dtype="double"):
     """
-    Create and load a dll corresponding to the source,info pair returned
+    Create and load a dll corresponding to the source, info pair returned
     from :func:`sasmodels.generate.make` compiled for the target precision.
 
@@ -199,5 +200,5 @@
         Npd2d = len(self.info['partype']['pd-2d'])
 
-        #print("dll",self.dllpath)
+        #print("dll", self.dllpath)
         try:
             self.dll = ct.CDLL(self.dllpath)
@@ -210,5 +211,5 @@
               else c_longdouble)
         pd_args_1d = [c_void_p, fp] + [c_int]*Npd1d if Npd1d else []
-        pd_args_2d= [c_void_p, fp] + [c_int]*Npd2d if Npd2d else []
+        pd_args_2d = [c_void_p, fp] + [c_int]*Npd2d if Npd2d else []
         self.Iq = self.dll[generate.kernel_name(self.info, False)]
         self.Iq.argtypes = IQ_ARGS + pd_args_1d + [fp]*Nfixed1d
@@ -218,28 +219,20 @@
 
     def __getstate__(self):
-        return {'info': self.info, 'dllpath': self.dllpath, 'dll': None}
+        return self.info, self.dllpath
 
     def __setstate__(self, state):
-        self.__dict__ = state
-
-    def __call__(self, q_input):
-        if self.dtype != q_input.dtype:
-            raise TypeError("data is %s kernel is %s" % (q_input.dtype, self.dtype))
+        self.info, self.dllpath = state
+        self.dll = None
+
+    def __call__(self, q_vectors):
+        q_input = PyInput(q_vectors, self.dtype)
         if self.dll is None: self._load_dll()
-        kernel = self.Iqxy if q_input.is_2D else self.Iq
+        kernel = self.Iqxy if q_input.is_2d else self.Iq
         return DllKernel(kernel, self.info, q_input)
 
-    # pylint: disable=no-self-use
-    def make_input(self, q_vectors):
-        """
-        Make q input vectors available to the model.
-
-        Note that each model needs its own q vector even if the case of
-        mixture models because some models may be OpenCL, some may be
-        ctypes and some may be pure python.
-        """
-        return PyInput(q_vectors, dtype=self.dtype)
-
     def release(self):
+        """
+        Release any resources associated with the model.
+        """
         pass # TODO: should release the dll
 
@@ -257,5 +250,5 @@
 
     The resulting call method takes the *pars*, a list of values for
-    the fixed parameters to the kernel, and *pd_pars*, a list of (value,weight)
+    the fixed parameters to the kernel, and *pd_pars*, a list of (value, weight)
     vectors for the polydisperse parameters.  *cutoff* determines the
     integration limits: any points with combined weight less than *cutoff*
@@ -269,5 +262,5 @@
         self.kernel = kernel
         self.res = np.empty(q_input.nq, q_input.dtype)
-        dim = '2d' if q_input.is_2D else '1d'
+        dim = '2d' if q_input.is_2d else '1d'
         self.fixed_pars = info['partype']['fixed-'+dim]
         self.pd_pars = info['partype']['pd-'+dim]
@@ -299,3 +292,6 @@
 
     def release(self):
+        """
+        Release any resources associated with the kernel.
+        """
         pass
Index: sasmodels/kernelpy.py
===================================================================
--- sasmodels/kernelpy.py	(revision 4c2c535bfe4a2a460e9aaf4ef0ed8b22d96965fd)
+++ sasmodels/kernelpy.py	(revision eafc9fa8b200a3b754e7dafb960c8891361eacef)
@@ -1,2 +1,10 @@
+"""
+Python driver for python kernels
+
+Calls the kernel with a vector of $q$ values for a single parameter set.
+Polydispersity is supported by looping over different parameter sets and
+summing the results.  The interface to :class:`PyModel` matches those for
+:class:`kernelcl.GpuModel` and :class:`kerneldll.DllModel`.
+"""
 import numpy as np
 from numpy import pi, cos
@@ -5,16 +13,19 @@
 
 class PyModel(object):
+    """
+    Wrapper for pure python models.
+    """
     def __init__(self, info):
         self.info = info
 
-    def __call__(self, q_input):
-        kernel = self.info['Iqxy'] if q_input.is_2D else self.info['Iq']
+    def __call__(self, q_vectors):
+        q_input = PyInput(q_vectors, dtype=F64)
+        kernel = self.info['Iqxy'] if q_input.is_2d else self.info['Iq']
         return PyKernel(kernel, self.info, q_input)
 
-    # pylint: disable=no-self-use
-    def make_input(self, q_vectors):
-        return PyInput(q_vectors, dtype=F64)
-
     def release(self):
+        """
+        Free resources associated with the model.
+        """
         pass
 
@@ -41,9 +52,12 @@
         self.nq = q_vectors[0].size
         self.dtype = dtype
-        self.is_2D = (len(q_vectors) == 2)
+        self.is_2d = (len(q_vectors) == 2)
         self.q_vectors = [np.ascontiguousarray(q, self.dtype) for q in q_vectors]
         self.q_pointers = [q.ctypes.data for q in self.q_vectors]
 
     def release(self):
+        """
+        Free resources associated with the model inputs.
+        """
         self.q_vectors = []
 
@@ -71,5 +85,5 @@
         self.q_input = q_input
         self.res = np.empty(q_input.nq, q_input.dtype)
-        dim = '2d' if q_input.is_2D else '1d'
+        dim = '2d' if q_input.is_2d else '1d'
         # Loop over q unless user promises that the kernel is vectorized by
         # taggining it with vectorized=True
@@ -77,8 +91,14 @@
             if dim == '2d':
                 def vector_kernel(qx, qy, *args):
+                    """
+                    Vectorized 2D kernel.
+                    """
                     return np.array([kernel(qxi, qyi, *args)
                                      for qxi, qyi in zip(qx, qy)])
             else:
                 def vector_kernel(q, *args):
+                    """
+                    Vectorized 1D kernel.
+                    """
                     return np.array([kernel(qi, *args)
                                      for qi in q])
@@ -122,4 +142,7 @@
 
     def release(self):
+        """
+        Free resources associated with the kernel.
+        """
         self.q_input = None
 
Index: sasmodels/list_pars.py
===================================================================
--- sasmodels/list_pars.py	(revision 5c962dff13d0a0a2f6d779be8df775f540958dbc)
+++ sasmodels/list_pars.py	(revision 5c962dff13d0a0a2f6d779be8df775f540958dbc)
@@ -0,0 +1,63 @@
+"""
+List all parameters used along with the models which use them.
+
+Usage:
+
+    python -m sasmodels.list_pars [-v]
+
+If '-v' is given, then list the models containing the parameter in
+addition to just the parameter name.
+"""
+from __future__ import print_function
+
+import sys
+
+from .core import load_model_definition
+from .generate import make_info
+from .compare import MODELS, columnize
+
+def find_pars():
+    """
+    Find all parameters in all models.
+
+    Returns the reference table *{parameter: [model, model, ...]}*
+    """
+    partable = {}
+    for name in sorted(MODELS):
+        definition = load_model_definition(name)
+        info = make_info(definition)
+        for p in info['parameters']:
+            pname = p[0]
+            partable.setdefault(pname, [])
+            partable[pname].append(name)
+    return partable
+
+def list_pars(names_only=True):
+    """
+    Print all parameters in all models.
+
+    If *names_only* then only print the parameter name, not the models it
+    occurs in.
+    """
+    partable = find_pars()
+    if names_only:
+        print(columnize(list(sorted(partable.keys()))))
+    else:
+        for k, v in sorted(partable.items()):
+            print("%s: %s"%(k, ", ".join(v)))
+
+def main():
+    """
+    Program to list the parameters used across all models.
+    """
+    if len(sys.argv) == 2 and sys.argv[1] == '-v':
+        verbose = True
+    elif len(sys.argv) == 1:
+        verbose = False
+    else:
+        print(__doc__)
+        sys.exit(1)
+    list_pars(names_only=not verbose)
+
+if __name__ == "__main__":
+    main()
Index: sasmodels/model_test.py
===================================================================
--- sasmodels/model_test.py	(revision 9404dd3f599fcbe72790aae8a6d176e64f292ea6)
+++ sasmodels/model_test.py	(revision 5c962dff13d0a0a2f6d779be8df775f540958dbc)
@@ -9,5 +9,5 @@
     if model1 is 'all', then all except the remaining models will be tested
 
-Each model is tested using the default parameters at q=0.1, (qx,qy)=(0.1,0.1),
+Each model is tested using the default parameters at q=0.1, (qx, qy)=(0.1, 0.1),
 and the ER and VR are computed.  The return values at these points are not
 considered.  The test is only to verify that the models run to completion,
@@ -30,5 +30,5 @@
         [ {parameters}, (qx, qy), I(qx, Iqy)],
         [ {parameters}, [(qx1, qy1), (qx2, qy2), ...],
-                        [I(qx1,qy1), I(qx2,qy2), ...]],
+                        [I(qx1, qy1), I(qx2, qy2), ...]],
 
         [ {parameters}, 'ER', ER(pars) ],
@@ -43,4 +43,5 @@
 Precision defaults to 5 digits (relative).
 """
+from __future__ import print_function
 
 import sys
@@ -55,4 +56,13 @@
 
 def make_suite(loaders, models):
+    """
+    Construct the pyunit test suite.
+
+    *loaders* is the list of kernel drivers to use, which is one of
+    *["dll", "opencl"]*, *["dll"]* or *["opencl"]*.  For python models,
+    the python driver is always used.
+
+    *models* is the list of models to test, or *["all"]* to test all models.
+    """
 
     ModelTestCase = _hide_model_case_from_nosetests()
@@ -75,5 +85,5 @@
         # don't try to call cl kernel since it will not be
         # available in some environmentes.
-        is_py = callable(getattr(model_definition,'Iq', None))
+        is_py = callable(getattr(model_definition, 'Iq', None))
 
         if is_py:  # kernel implemented in python
@@ -111,4 +121,11 @@
 def _hide_model_case_from_nosetests():
     class ModelTestCase(unittest.TestCase):
+        """
+        Test suit for a particular model with a particular kernel driver.
+
+        The test suite runs a simple smoke test to make sure the model
+        functions, then runs the list of tests at the bottom of the model
+        description file.
+        """
         def __init__(self, test_name, definition, test_method_name,
                      platform, dtype):
@@ -123,8 +140,8 @@
         def _runTest(self):
             smoke_tests = [
-                [{},0.1,None],
-                [{},(0.1,0.1),None],
-                [{},'ER',None],
-                [{},'VR',None],
+                [{}, 0.1, None],
+                [{}, (0.1, 0.1), None],
+                [{}, 'ER', None],
+                [{}, 'VR', None],
                 ]
 
@@ -163,5 +180,5 @@
                 actual = [call_VR(model.info, pars)]
             elif isinstance(x[0], tuple):
-                Qx,Qy = zip(*x)
+                Qx, Qy = zip(*x)
                 q_vectors = [np.array(Qx), np.array(Qy)]
                 kernel = make_kernel(model, q_vectors)
@@ -179,10 +196,11 @@
                     # smoke test --- make sure it runs and produces a value
                     self.assertTrue(np.isfinite(actual_yi),
-                        'invalid f(%s): %s' % (xi, actual_yi))
+                                    'invalid f(%s): %s' % (xi, actual_yi))
                 else:
                     err = abs(yi - actual_yi)
                     nrm = abs(yi)
                     self.assertLess(err * 10**5, nrm,
-                        'f(%s); expected:%s; actual:%s' % (xi, yi, actual_yi))
+                                    'f(%s); expected:%s; actual:%s'
+                                    % (xi, yi, actual_yi))
 
     return ModelTestCase
@@ -237,5 +255,5 @@
     Run "nosetests sasmodels" on the command line to invoke it.
     """
-    tests = make_suite(['opencl','dll'],['all'])
+    tests = make_suite(['opencl', 'dll'], ['all'])
     for test_i in tests:
         yield test_i._runTest
Index: sasmodels/resolution.py
===================================================================
--- sasmodels/resolution.py	(revision 190fc2b167d3d89ca9c3b7d3128bdf8bde9f67d0)
+++ sasmodels/resolution.py	(revision 5925e909151e120132e72e91095ab06bad045240)
@@ -14,5 +14,5 @@
            "pinhole_extend_q", "slit_extend_q", "bin_edges",
            "interpolate", "linear_extrapolation", "geometric_extrapolation",
-           ]
+          ]
 
 MINIMUM_RESOLUTION = 1e-8
@@ -80,8 +80,8 @@
         # default to Perfect1D if the pinhole geometry is not defined.
         self.q, self.q_width = q, q_width
-        self.q_calc = pinhole_extend_q(q, q_width, 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.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))
 
     def apply(self, theory):
@@ -470,4 +470,11 @@
 
 def eval_form(q, form, pars):
+    """
+    Return the SAS model evaluated at *q*.
+
+    *form* is the SAS model returned from :fun:`core.load_model`.
+
+    *pars* are the parameter values to use when evaluating.
+    """
     from sasmodels import core
     kernel = core.make_kernel(form, [q])
@@ -478,4 +485,9 @@
 
 def gaussian(q, q0, dq):
+    """
+    Return the Gaussian resolution function.
+
+    *q0* is the center, *dq* is the width and *q* are the points to evaluate.
+    """
     from numpy import exp, pi
     return exp(-0.5*((q-q0)/dq)**2)/(sqrt(2*pi)*dq)
@@ -559,4 +571,7 @@
 
 class ResolutionTest(unittest.TestCase):
+    """
+    Test the resolution calculations.
+    """
 
     def setUp(self):
@@ -668,4 +683,7 @@
 
 class IgorComparisonTest(unittest.TestCase):
+    """
+    Test resolution calculations against those returned by Igor.
+    """
 
     def setUp(self):
@@ -675,5 +693,5 @@
         self.model = core.load_model(sphere, dtype='double')
 
-    def Iq_sphere(self, pars, resolution):
+    def _eval_sphere(self, pars, resolution):
         from sasmodels import core
         kernel = core.make_kernel(self.model, [resolution.q_calc])
@@ -683,5 +701,5 @@
         return result
 
-    def compare(self, q, output, answer, tolerance):
+    def _compare(self, q, output, answer, tolerance):
         #err = (output - answer)/answer
         #idx = abs(err) >= tolerance
@@ -700,6 +718,6 @@
         q, width, answer, _ = data
         resolution = Perfect1D(q)
-        output = self.Iq_sphere(pars, resolution)
-        self.compare(q, output, answer, 1e-6)
+        output = self._eval_sphere(pars, resolution)
+        self._compare(q, output, answer, 1e-6)
 
     def test_pinhole(self):
@@ -713,7 +731,7 @@
         q, q_width, answer = data
         resolution = Pinhole1D(q, q_width)
-        output = self.Iq_sphere(pars, resolution)
+        output = self._eval_sphere(pars, resolution)
         # TODO: relative error should be lower
-        self.compare(q, output, answer, 3e-4)
+        self._compare(q, output, answer, 3e-4)
 
     def test_pinhole_romberg(self):
@@ -736,14 +754,14 @@
         q_calc, tol = None, 0.01
         resolution = Pinhole1D(q, q_width, q_calc=q_calc)
-        output = self.Iq_sphere(pars, resolution)
+        output = self._eval_sphere(pars, resolution)
         if 0: # debug plot
             import matplotlib.pyplot as plt
             resolution = Perfect1D(q)
-            source = self.Iq_sphere(pars, resolution)
+            source = self._eval_sphere(pars, resolution)
             plt.loglog(q, source, '.')
             plt.loglog(q, answer, '-', hold=True)
             plt.loglog(q, output, '-', hold=True)
             plt.show()
-        self.compare(q, output, answer, tol)
+        self._compare(q, output, answer, tol)
 
     def test_slit(self):
@@ -757,9 +775,9 @@
         q, delta_qv, _, answer = data
         resolution = Slit1D(q, width=delta_qv, height=0)
-        output = self.Iq_sphere(pars, resolution)
+        output = self._eval_sphere(pars, resolution)
         # TODO: eliminate Igor test since it is too inaccurate to be useful.
         # This means we can eliminate the test data as well, and instead
         # use a generated q vector.
-        self.compare(q, output, answer, 0.5)
+        self._compare(q, output, answer, 0.5)
 
     def test_slit_romberg(self):
@@ -777,7 +795,7 @@
                                delta_qv[0], 0.)
         resolution = Slit1D(q, width=delta_qv, height=0, q_calc=q_calc)
-        output = self.Iq_sphere(pars, resolution)
+        output = self._eval_sphere(pars, resolution)
         # TODO: relative error should be lower
-        self.compare(q, output, answer, 0.025)
+        self._compare(q, output, answer, 0.025)
 
     def test_ellipsoid(self):
@@ -799,5 +817,5 @@
         # TODO: 10% is too much error; use better algorithm
         #print(np.max(abs(answer-output)/answer))
-        self.compare(q, output, answer, 0.1)
+        self._compare(q, output, answer, 0.1)
 
     #TODO: can sas q spacing be too sparse for the resolution calculation?
@@ -813,6 +831,6 @@
         q, q_width, answer = data[:, ::20] # Take every nth point
         resolution = Pinhole1D(q, q_width)
-        output = self.Iq_sphere(pars, resolution)
-        self.compare(q, output, answer, 1e-6)
+        output = self._eval_sphere(pars, resolution)
+        self._compare(q, output, answer, 1e-6)
 
 
@@ -1067,4 +1085,7 @@
 
 def demo_pinhole_1d():
+    """
+    Show example of pinhole smearing.
+    """
     q = np.logspace(-4, np.log10(0.2), 400)
     q_width = 0.1*q
@@ -1073,4 +1094,7 @@
 
 def demo_slit_1d():
+    """
+    Show example of slit smearing.
+    """
     q = np.logspace(-4, np.log10(0.2), 100)
     w = h = 0.
@@ -1083,4 +1107,7 @@
 
 def demo():
+    """
+    Run the resolution demos.
+    """
     import matplotlib.pyplot as plt
     plt.subplot(121)
Index: sasmodels/resolution2d.py
===================================================================
--- sasmodels/resolution2d.py	(revision 841753cbac7824298effa1aec88a7ee977eb257f)
+++ sasmodels/resolution2d.py	(revision 823e6204f35d6c2bf32f27a3a3c137a369a973cf)
@@ -143,11 +143,11 @@
         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))
+            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)
+            qx_res = qx + dqx*cos(dphi)
+            qy_res = qy + dqy*sin(dphi)
 
 
@@ -162,5 +162,5 @@
             theory = np.reshape(theory, (nbins, nq))
             ## Averaging with Gaussian weighting: normalization included.
-            value =np.average(theory, axis=0, weights=self.q_calc_weights)
+            value = np.average(theory, axis=0, weights=self.q_calc_weights)
             ## Return the smeared values in the range of self.index
             return value
Index: sasmodels/sasview_model.py
===================================================================
--- sasmodels/sasview_model.py	(revision 9404dd3f599fcbe72790aae8a6d176e64f292ea6)
+++ sasmodels/sasview_model.py	(revision eafc9fa8b200a3b754e7dafb960c8891361eacef)
@@ -283,5 +283,5 @@
         """
         q_vectors = [np.asarray(q) for q in args]
-        fn = self._model(self._model.make_input(q_vectors))
+        fn = self._model(q_vectors)
         pars = [self.params[v] for v in fn.fixed_pars]
         pd_pars = [self._get_weights(p) for p in fn.pd_pars]
Index: sasmodels/weights.py
===================================================================
--- sasmodels/weights.py	(revision 3c56da87f0995bdde77f4183b68e77b81584dde0)
+++ sasmodels/weights.py	(revision 5c962dff13d0a0a2f6d779be8df775f540958dbc)
@@ -22,4 +22,7 @@
 
     def get_pars(self):
+        """
+        Return the parameters to the disperser as a dictionary.
+        """
         pars = {'type': self.type}
         pars.update(self.__dict__)
@@ -28,4 +31,7 @@
     # pylint: disable=no-self-use
     def set_weights(self, values, weights):
+        """
+        Set the weights on the disperser if it is :class:`ArrayDispersion`.
+        """
         raise RuntimeError("set_weights is only available for ArrayDispersion")
 
@@ -36,5 +42,5 @@
         *center* is the center of the distribution
 
-        *lb*,*ub* are the min and max allowed values
+        *lb*, *ub* are the min and max allowed values
 
         *relative* is True if the distribution width is proportional to the
@@ -63,4 +69,11 @@
 
 class GaussianDispersion(Dispersion):
+    r"""
+    Gaussian dispersion, with 1-\ $\sigma$ width.
+
+    .. math::
+
+        w = \exp\left(-\tfrac12 (x - c)^2/\sigma^2\right)
+    """
     type = "gaussian"
     default = dict(npts=35, width=0, nsigmas=3)
@@ -72,4 +85,11 @@
 
 class RectangleDispersion(Dispersion):
+    r"""
+    Uniform dispersion, with width $\sqrt{3}\sigma$.
+
+    .. math::
+
+        w = 1
+    """
     type = "rectangle"
     default = dict(npts=35, width=0, nsigmas=1.70325)
@@ -81,18 +101,40 @@
 
 class LogNormalDispersion(Dispersion):
+    r"""
+    log Gaussian dispersion, with 1-\ $\sigma$ width.
+
+    .. math::
+
+        w = \frac{\exp\left(-\tfrac12 (\ln x - c)^2/\sigma^2\right)}{x\sigma}
+    """
     type = "lognormal"
     default = dict(npts=80, width=0, nsigmas=8)
     def _weights(self, center, sigma, lb, ub):
-        x = self._linspace(center, sigma, max(lb,1e-8), max(ub,1e-8))
-        px = np.exp(-0.5*(np.log(x)-center)**2)/sigma**2/(x*sigma)
+        x = self._linspace(center, sigma, max(lb, 1e-8), max(ub, 1e-8))
+        px = np.exp(-0.5*(np.log(x)-center)**2/sigma**2)/(x*sigma)
         return x, px
 
 
 class SchulzDispersion(Dispersion):
+    r"""
+    Schultz dispersion, with 1-\ $\sigma$ width.
+
+    .. math::
+
+        w = \frac{z^z\,R^{z-1}}{e^{Rz}\,c \Gamma(z)}
+
+    where $c$ is the center of the distribution, $R = x/c$ and $z=(c/\sigma)^2$.
+
+    This is evaluated using logarithms as
+
+    .. math::
+
+        w = \exp\left(z \ln z + (z-1)\ln R - Rz - \ln c - \ln \Gamma(z) \right)
+    """
     type = "schulz"
     default = dict(npts=80, width=0, nsigmas=8)
     def _weights(self, center, sigma, lb, ub):
-        x = self._linspace(center, sigma, max(lb,1e-8), max(ub,1e-8))
-        R= x/center
+        x = self._linspace(center, sigma, max(lb, 1e-8), max(ub, 1e-8))
+        R = x/center
         z = (center/sigma)**2
         arg = z*np.log(z) + (z-1)*np.log(R) - R*z - np.log(center) - gammaln(z)
@@ -102,4 +144,9 @@
 
 class ArrayDispersion(Dispersion):
+    r"""
+    Empirical dispersion curve.
+
+    Use :meth:`set_weights` to set $w = f(x)$.
+    """
     type = "array"
     default = dict(npts=35, width=0, nsigmas=1)
@@ -110,4 +157,7 @@
 
     def set_weights(self, values, weights):
+        """
+        Set the weights for the given x values.
+        """
         self.values = np.ascontiguousarray(values, 'd')
         self.weights = np.ascontiguousarray(weights, 'd')
@@ -117,5 +167,5 @@
         # TODO: interpolate into the array dispersion using npts
         x = center + self.values*sigma
-        idx = (x>=lb)&(x<=ub)
+        idx = (x >= lb) & (x <= ub)
         x = x[idx]
         px = self.weights[idx]
@@ -124,5 +174,5 @@
 
 # dispersion name -> disperser lookup table.
-models = dict((d.type,d) for d in (
+models = dict((d.type, d) for d in (
     GaussianDispersion, RectangleDispersion,
     ArrayDispersion, SchulzDispersion, LogNormalDispersion
@@ -149,13 +199,13 @@
     of the parameter, and false if it is an absolute width.
 
-    Returns *(value,weight)*, where *value* and *weight* are vectors.
+    Returns *(value, weight)*, where *value* and *weight* are vectors.
     """
     cls = models[disperser]
     obj = cls(n, width, nsigmas)
-    v,w = obj.get_weights(value, limits[0], limits[1], relative)
-    return v,w
+    v, w = obj.get_weights(value, limits[0], limits[1], relative)
+    return v, w
 
 # Hack to allow sasview dispersion objects to interoperate with sasmodels
-dispersers = dict((v.__name__,k) for k,v in models.items())
+dispersers = dict((v.__name__, k) for k, v in models.items())
 dispersers['DispersionModel'] = RectangleDispersion.type
 
