Index: extra/pylint.rc
===================================================================
--- extra/pylint.rc	(revision 37a725200a1edb6dbd7a029a236c3ea2cfc7522c)
+++ extra/pylint.rc	(revision 5ef0633b151126a2a8e376666e05ccbcc0341316)
@@ -6,5 +6,5 @@
 # Python code to execute, usually for sys.path manipulation such as
 # pygtk.require().
-#init-hook='import os, sys; sys.path.extend([os.getcwd(),os.path.join(os.getcwd(),"extra")])'
+#init-hook=
 
 # Profiled execution.
Index: sasmodels/__init__.py
===================================================================
--- sasmodels/__init__.py	(revision 37a725200a1edb6dbd7a029a236c3ea2cfc7522c)
+++ sasmodels/__init__.py	(revision 989005365096670e54afec254875979f35bd395a)
@@ -1,17 +1,1 @@
-"""
-sasmodels
-=========
-
-**sasmodels** is a package containing models for small angle neutron and xray
-scattering.  Models supported are the one dimensional circular average and
-two dimensional oriented patterns.  As well as the form factor calculations
-for the individual shapes **sasmodels** also provides automatic shape
-polydispersity, angular dispersion and resolution convolution.  SESANS
-patterns can be computed for any model.
-
-Models can be written in python or in C.  C models can run on the GPU if
-OpenCL drivers are available.  See :mod:`generate` for details on
-defining new models.
-"""
-
 __version__ = "0.9"
Index: sasmodels/bumps_model.py
===================================================================
--- sasmodels/bumps_model.py	(revision 37a725200a1edb6dbd7a029a236c3ea2cfc7522c)
+++ sasmodels/bumps_model.py	(revision ec7e360eb9b8ec7bb8f5006ed71720a9fa7d593b)
@@ -6,13 +6,9 @@
 arguments to set the initial value for each parameter.
 
-:class:`Experiment` combines the *Model* function with a data file loaded by
-the sasview data loader.  *Experiment* takes a *cutoff* parameter controlling
+:class:`Experiment` combines the *Model* function with a data file loaded by the
+sasview data loader.  *Experiment* takes a *cutoff* parameter controlling
 how far the polydispersity integral extends.
 
 """
-
-__all__ = [
-    "Model", "Experiment",
-    ]
 
 import warnings
@@ -24,60 +20,18 @@
 
 # CRUFT: old style bumps wrapper which doesn't separate data and model
-# pylint: disable=invalid-name
 def BumpsModel(data, model, cutoff=1e-5, **kw):
-    r"""
-    Bind a model to data, along with a polydispersity cutoff.
-
-    *data* is a :class:`data.Data1D`, :class:`data.Data2D` or
-    :class:`data.Sesans` object.  Use :func:`data.empty_data1D` or
-    :func:`data.empty_data2D` to define $q, \Delta q$ calculation
-    points for displaying the SANS curve when there is no measured data.
-
-    *model* is a runnable module as returned from :func:`core.load_model`.
-
-    *cutoff* is the polydispersity weight cutoff.
-
-    Any additional *key=value* pairs are model dependent parameters.
-
-    Returns an :class:`Experiment` object.
-
-    Note that the usual Bumps semantics is not fully supported, since
-    assigning *M.name = parameter* on the returned experiment object
-    does not set that parameter in the model.  Range setting will still
-    work as expected though.
-
-    .. deprecated:: 0.1
-        Use :class:`Experiment` instead.
-    """
     warnings.warn("Use of BumpsModel is deprecated.  Use bumps_model.Experiment instead.")
-
-    # Create the model and experiment
     model = Model(model, **kw)
     experiment = Experiment(data=data, model=model, cutoff=cutoff)
-
-    # Copy the model parameters up to the experiment object.
-    for k, v in model.parameters().items():
-        setattr(experiment, k, v)
+    for k in model._parameter_names:
+        setattr(experiment, k, getattr(model, k))
     return experiment
 
-
 def create_parameters(model_info, **kwargs):
-    """
-    Generate Bumps parameters from the model info.
-
-    *model_info* is returned from :func:`generate.model_info` on the
-    model definition module.
-
-    Any additional *key=value* pairs are initial values for the parameters
-    to the models.  Uninitialized parameters will use the model default
-    value.
-
-    Returns a dictionary of *{name: Parameter}* containing the bumps
-    parameters for each model parameter, and a dictionary of
-    *{name: str}* containing the polydispersity distribution types.
-    """
     # lazy import; this allows the doc builder and nosetests to run even
     # when bumps is not on the path.
     from bumps.names import Parameter
+
+    partype = model_info['partype']
 
     pars = {}
@@ -86,10 +40,10 @@
         value = kwargs.pop(name, default)
         pars[name] = Parameter.default(value, name=name, limits=limits)
-    for name in model_info['partype']['pd-2d']:
+    for name in partype['pd-2d']:
         for xpart, xdefault, xlimits in [
-                ('_pd', 0., limits),
-                ('_pd_n', 35., (0, 1000)),
-                ('_pd_nsigma', 3., (0, 10)),
-            ]:
+            ('_pd', 0., limits),
+            ('_pd_n', 35., (0, 1000)),
+            ('_pd_nsigma', 3., (0, 10)),
+        ]:
             xname = name + xpart
             xvalue = kwargs.pop(xname, xdefault)
@@ -97,5 +51,5 @@
 
     pd_types = {}
-    for name in model_info['partype']['pd-2d']:
+    for name in partype['pd-2d']:
         xname = name + '_pd_type'
         xvalue = kwargs.pop(xname, 'gaussian')
@@ -109,19 +63,10 @@
 
 class Model(object):
-    """
-    Bumps wrapper for a SAS model.
-
-    *model* is a runnable module as returned from :func:`core.load_model`.
-
-    *cutoff* is the polydispersity weight cutoff.
-
-    Any additional *key=value* pairs are model dependent parameters.
-    """
     def __init__(self, model, **kwargs):
         self._sasmodel = model
         pars, pd_types = create_parameters(model.info, **kwargs)
-        for k, v in pars.items():
+        for k,v in pars.items():
             setattr(self, k, v)
-        for k, v in pd_types.items():
+        for k,v in pd_types.items():
             setattr(self, k, v)
         self._parameter_names = list(pars.keys())
@@ -130,14 +75,9 @@
     def parameters(self):
         """
-        Return a dictionary of parameters objects for the parameters,
-        excluding polydispersity distribution type.
+        Return a dictionary of parameters
         """
         return dict((k, getattr(self, k)) for k in self._parameter_names)
 
     def state(self):
-        """
-        Return a dictionary of current values for all the parameters,
-        including polydispersity distribution type.
-        """
         pars = dict((k, getattr(self, k).value) for k in self._parameter_names)
         pars.update((k, getattr(self, k)) for k in self._pd_type_names)
@@ -145,18 +85,18 @@
 
 class Experiment(DataMixin):
-    r"""
-    Bumps wrapper for a SAS experiment.
+    """
+    Return a bumps wrapper for a SAS model.
 
-    *data* is a :class:`data.Data1D`, :class:`data.Data2D` or
-    :class:`data.Sesans` object.  Use :func:`data.empty_data1D` or
-    :func:`data.empty_data2D` to define $q, \Delta q$ calculation
-    points for displaying the SANS curve when there is no measured data.
+    *data* is the data to be fitted.
 
-    *model* is a :class:`Model` object.
+    *model* is the SAS model from :func:`core.load_model`.
 
     *cutoff* is the integration cutoff, which avoids computing the
     the SAS model where the polydispersity weight is low.
 
-    The resulting model can be used directly in a Bumps FitProblem call.
+    Model parameters can be initialized with additional keyword
+    arguments, or by assigning to model.parameter_name.value.
+
+    The resulting bumps model can be used directly in a FitProblem call.
     """
     def __init__(self, data, model, cutoff=1e-5):
@@ -169,13 +109,9 @@
 
     def update(self):
-        """
-        Call when model parameters have changed and theory needs to be
-        recalculated.
-        """
         self._cache = {}
 
     def numpoints(self):
         """
-        Return the number of data points
+            Return the number of points
         """
         return len(self.Iq)
@@ -188,10 +124,4 @@
 
     def theory(self):
-        """
-        Return the theory corresponding to the model parameters.
-
-        This method uses lazy evaluation, and requires model.update() to be
-        called when the parameters have changed.
-        """
         if 'theory' not in self._cache:
             pars = self.model.state()
@@ -200,16 +130,8 @@
 
     def residuals(self):
-        """
-        Return theory minus data normalized by uncertainty.
-        """
         #if np.any(self.err ==0): print("zeros in err")
         return (self.theory() - self.Iq) / self.dIq
 
     def nllf(self):
-        """
-        Return the negative log likelihood of seeing data given the model
-        parameters, up to a normalizing constant which depends on the data
-        uncertainty.
-        """
         delta = self.residuals()
         #if np.any(np.isnan(R)): print("NaN in residuals")
@@ -227,16 +149,8 @@
 
     def simulate_data(self, noise=None):
-        """
-        Generate simulated data.
-        """
         Iq = self.theory()
         self._set_data(Iq, noise)
 
     def save(self, basename):
-        """
-        Save the model parameters and data into a file.
-
-        Not Implemented.
-        """
         pass
 
Index: sasmodels/compare.py
===================================================================
--- sasmodels/compare.py	(revision d15a908796935025260dc12a6d16932a079552d3)
+++ sasmodels/compare.py	(revision 608e31e3a7d22b75d2b3b9ec919739d9345a9afb)
@@ -72,11 +72,9 @@
 # doc string so that we can display it at run time if there is an error.
 # lin
-__doc__ = (__doc__  # pylint: disable=redefined-builtin
-           + """
+__doc__ = __doc__ + """
 Program description
 -------------------
 
-"""
-           + USAGE)
+""" + USAGE
 
 
@@ -283,6 +281,5 @@
             theory = lambda: smearer.get_value()
         else:
-            theory = lambda: model.evalDistribution([data.qx_data[index],
-                                                     data.qy_data[index]])
+            theory = lambda: model.evalDistribution([data.qx_data[index], data.qy_data[index]])
     elif smearer is not None:
         theory = lambda: smearer(model.evalDistribution(data.x))
@@ -419,6 +416,5 @@
         try:
             base_value, base_time = time_calculation(base, pars, Nbase)
-            print("%s t=%.1f ms, intensity=%.0f"
-                  % (base.engine, base_time, sum(base_value)))
+            print("%s t=%.1f ms, intensity=%.0f"%(base.engine, base_time, sum(base_value)))
         except ImportError:
             traceback.print_exc()
@@ -430,6 +426,5 @@
         try:
             comp_value, comp_time = time_calculation(comp, pars, Ncomp)
-            print("%s t=%.1f ms, intensity=%.0f"
-                  % (comp.engine, comp_time, sum(comp_value)))
+            print("%s t=%.1f ms, intensity=%.0f"%(comp.engine, comp_time, sum(comp_value)))
         except ImportError:
             traceback.print_exc()
@@ -438,11 +433,12 @@
     # Compare, but only if computing both forms
     if Nbase > 0 and Ncomp > 0:
+        #print("speedup %.2g"%(comp_time/base_time))
+        #print("max |base/comp|", max(abs(base_value/comp_value)), "%.15g"%max(abs(base_value)), "%.15g"%max(abs(comp_value)))
+        #comp *= max(base_value/comp_value)
         resid = (base_value - comp_value)
         relerr = resid/comp_value
-        _print_stats("|%s-%s|"
-                     % (base.engine, comp.engine) + (" "*(3+len(comp.engine))),
+        _print_stats("|%s-%s|"%(base.engine, comp.engine) + (" "*(3+len(comp.engine))),
                      resid)
-        _print_stats("|(%s-%s)/%s|"
-                     % (base.engine, comp.engine, comp.engine),
+        _print_stats("|(%s-%s)/%s|"%(base.engine, comp.engine, comp.engine),
                      relerr)
 
@@ -595,5 +591,5 @@
     invalid = [o[1:] for o in flags
                if o[1:] not in NAME_OPTIONS
-               and not any(o.startswith('-%s='%t) for t in VALUE_OPTIONS)]
+                   and not any(o.startswith('-%s='%t) for t in VALUE_OPTIONS)]
     if invalid:
         print("Invalid options: %s"%(", ".join(invalid)))
@@ -601,5 +597,4 @@
 
 
-    # pylint: disable=bad-whitespace
     # Interpret the flags
     opts = {
@@ -656,5 +651,4 @@
         elif arg == '-sasview': engines.append(arg[1:])
         elif arg == '-edit':    opts['explore'] = True
-    # pylint: enable=bad-whitespace
 
     if len(engines) == 0:
@@ -681,9 +675,9 @@
     presets = {}
     for arg in values:
-        k, v = arg.split('=', 1)
+        k,v = arg.split('=',1)
         if k not in pars:
             # extract base name without polydispersity info
             s = set(p.split('_pd')[0] for p in pars)
-            print("%r invalid; parameters are: %s"%(k, ", ".join(sorted(s))))
+            print("%r invalid; parameters are: %s"%(k,", ".join(sorted(s))))
             sys.exit(1)
         presets[k] = float(v) if not k.endswith('type') else v
@@ -703,5 +697,5 @@
 
     # Create the computational engines
-    data, _ = make_data(opts)
+    data, _index = make_data(opts)
     if n1:
         base = make_engine(model_definition, data, engines[0], opts['cutoff'])
@@ -713,5 +707,4 @@
         comp = None
 
-    # pylint: disable=bad-whitespace
     # Remember it all
     opts.update({
@@ -725,12 +718,15 @@
         'engines'   : [base, comp],
     })
-    # pylint: enable=bad-whitespace
 
     return opts
 
+def main():
+    opts = parse_opts()
+    if opts['explore']:
+        explore(opts)
+    else:
+        compare(opts)
+
 def explore(opts):
-    """
-    Explore the model using the Bumps GUI.
-    """
     import wx
     from bumps.names import FitProblem
@@ -738,20 +734,17 @@
 
     problem = FitProblem(Explore(opts))
-    is_mac = "cocoa" in wx.version()
+    isMac = "cocoa" in wx.version()
     app = wx.App()
     frame = AppFrame(parent=None, title="explore")
-    if not is_mac: frame.Show()
+    if not isMac: frame.Show()
     frame.panel.set_model(model=problem)
     frame.panel.Layout()
     frame.panel.aui.Split(0, wx.TOP)
-    if is_mac: frame.Show()
+    if isMac: frame.Show()
     app.MainLoop()
 
 class Explore(object):
     """
-    Bumps wrapper for a SAS model comparison.
-
-    The resulting object can be used as a Bumps fit problem so that
-    parameters can be adjusted in the GUI, with plots updated on the fly.
+    Return a bumps wrapper for a SAS model comparison.
     """
     def __init__(self, opts):
@@ -794,5 +787,4 @@
         Return cost.
         """
-        # pylint: disable=no-self-use
         return 0.  # No nllf
 
@@ -812,14 +804,4 @@
 
 
-def main():
-    """
-    Main program.
-    """
-    opts = parse_opts()
-    if opts['explore']:
-        explore(opts)
-    else:
-        compare(opts)
-
 if __name__ == "__main__":
     main()
Index: sasmodels/compare_many.py
===================================================================
--- sasmodels/compare_many.py	(revision d15a908796935025260dc12a6d16932a079552d3)
+++ sasmodels/compare_many.py	(revision 6458608d36e65099c9e90af3b188b8a10cadb474)
@@ -1,15 +1,3 @@
 #!/usr/bin/env python
-"""
-Program to compare results from many random parameter sets for a given model.
-
-The result is a comma separated value (CSV) table that can be redirected
-from standard output into a file and loaded into a spreadsheet.
-
-The models are compared for each parameter set and if the difference is
-greater than expected for that precision, the parameter set is labeled
-as bad and written to the output, along with the random seed used to
-generate that parameter value.  This seed can be used with :mod:`compare`
-to reload and display the details of the model.
-"""
 from __future__ import print_function
 
@@ -26,38 +14,16 @@
 
 def calc_stats(target, value, index):
-    """
-    Calculate statistics between the target value and the computed value.
-
-    *target* and *value* are the vectors being compared, with the
-    difference normalized by *target* to get relative error.  Only
-    the elements listed in *index* are used, though index may be
-    and empty slice defined by *slice(None, None)*.
-
-    Returns:
-
-        *maxrel* the maximum relative difference
-
-        *rel95* the relative difference with the 5% biggest differences ignored
-
-        *maxabs* the maximum absolute difference for the 5% biggest differences
-
-        *maxval* the maximum value for the 5% biggest differences
-    """
     resid = abs(value-target)[index]
     relerr = resid/target[index]
-    sorted_rel_index = np.argsort(relerr)
+    srel = np.argsort(relerr)
     #p90 = int(len(relerr)*0.90)
     p95 = int(len(relerr)*0.95)
     maxrel = np.max(relerr)
-    rel95 = relerr[sorted_rel_index[p95]]
-    maxabs = np.max(resid[sorted_rel_index[p95:]])
-    maxval = np.max(value[sorted_rel_index[p95:]])
-    return maxrel, rel95, maxabs, maxval
+    rel95 = relerr[srel[p95]]
+    maxabs = np.max(resid[srel[p95:]])
+    maxval = np.max(value[srel[p95:]])
+    return maxrel,rel95,maxabs,maxval
 
 def print_column_headers(pars, parts):
-    """
-    Generate column headers for the differences and for the parameters,
-    and print them to standard output.
-    """
     stats = list('Max rel err|95% rel err|Max abs err above 90% rel|Max value above 90% rel'.split('|'))
     groups = ['']
@@ -70,5 +36,4 @@
     print(','.join('"%s"'%c for c in columns))
 
-# Target 'good' value for various precision levels.
 PRECISION = {
     'fast': 1e-3,
@@ -83,30 +48,13 @@
 def compare_instance(name, data, index, N=1, mono=True, cutoff=1e-5,
                      base='sasview', comp='double'):
-    r"""
-    Compare the model under different calculation engines.
-
-    *name* is the name of the model.
-
-    *data* is the data object giving $q, \Delta q$ calculation points.
-
-    *index* is the active set of points.
-
-    *N* is the number of comparisons to make.
-
-    *cutoff* is the polydispersity weight cutoff to make the calculation
-    a little bit faster.
-
-    *base* and *comp* are the names of the calculation engines to compare.
-    """
-
-    is_2d = hasattr(data, 'qx_data')
+    is2D = hasattr(data, 'qx_data')
     model_definition = core.load_model_definition(name)
     pars = get_demo_pars(model_definition)
     header = ('\n"Model","%s","Count","%d","Dimension","%s"'
-              % (name, N, "2D" if is_2d else "1D"))
+              % (name, N, "2D" if is2D else "1D"))
     if not mono: header += ',"Cutoff",%g'%(cutoff,)
     print(header)
 
-    if is_2d:
+    if is2D:
         info = generate.make_info(model_definition)
         partype = info['partype']
@@ -121,8 +69,4 @@
     # declarations are not available in python 2.7.
     def try_model(fn, pars):
-        """
-        Return the model evaluated at *pars*.  If there is an exception,
-        print it and return NaN of the right shape.
-        """
         try:
             result = fn(**pars)
@@ -138,8 +82,4 @@
         return result
     def check_model(pars):
-        """
-        Run the two calculators against *pars*, returning statistics
-        on the differences.  See :func:`calc_stats` for the list of stats.
-        """
         base_value = try_model(calc_base, pars)
         comp_value = try_model(calc_comp, pars)
@@ -168,5 +108,5 @@
         good = [True]
         columns = check_model(pars_i)
-        columns += [v for _, v in sorted(pars_i.items())]
+        columns += [v for _,v in sorted(pars_i.items())]
         if first:
             labels = [" vs. ".join((calc_base.engine, calc_comp.engine))]
@@ -181,21 +121,12 @@
 
 def print_usage():
-    """
-    Print the command usage string.
-    """
     print("usage: compare_many.py MODEL COUNT (1dNQ|2dNQ) (CUTOFF|mono) (single|double|quad)")
 
 
 def print_models():
-    """
-    Print the list of available models in columns.
-    """
     print(columnize(MODELS, indent="  "))
 
 
 def print_help():
-    """
-    Print usage string, the option description and the list of available models.
-    """
     print_usage()
     print("""\
@@ -227,8 +158,5 @@
 
 def main():
-    """
-    Main program.
-    """
-    if len(sys.argv) not in (6, 7):
+    if len(sys.argv) not in (6,7):
         print_help()
         sys.exit(1)
@@ -254,5 +182,5 @@
 
     data, index = make_data({'qmax':1.0, 'is2d':is2D, 'nq':Nq, 'res':0.,
-                             'accuracy': 'Low', 'view':'log'})
+                              'accuracy': 'Low', 'view':'log'})
     model_list = [model] if model != "all" else MODELS
     for model in model_list:
Index: sasmodels/convert.py
===================================================================
--- sasmodels/convert.py	(revision d15a908796935025260dc12a6d16932a079552d3)
+++ sasmodels/convert.py	(revision 29da2138b6792800b3b961f991743c21f14a1abd)
@@ -49,5 +49,5 @@
     """
     return dict((p, (v*1e6 if p.endswith('sld') else v))
-                for p, v in pars.items())
+                for p,v in pars.items())
 
 def convert_model(name, pars):
@@ -55,5 +55,5 @@
     Convert model from old style parameter names to new style.
     """
-    _, _ = name, pars # lint
+    _,_ = name,pars # lint
     raise NotImplementedError
     # need to load all new models in order to determine old=>new
@@ -67,5 +67,5 @@
     """
     return dict((p, (v*1e-6 if p.endswith('sld') else v))
-                for p, v in pars.items())
+                for p,v in pars.items())
 
 def _remove_pd(pars, key, name):
Index: sasmodels/core.py
===================================================================
--- sasmodels/core.py	(revision d15a908796935025260dc12a6d16932a079552d3)
+++ sasmodels/core.py	(revision cde11f0f857d1861f37ded681777699e3aa91738)
@@ -2,8 +2,6 @@
 Core model handling routines.
 """
-__all__ = [
-    "list_models", "load_model_definition", "precompile_dll",
-    "load_model", "make_kernel", "call_kernel", "call_ER", "call_VR",
-    ]
+__all__ = ["list_models", "load_model_definition",  "precompile_dll",
+           "load_model", "make_kernel", "call_kernel", "call_ER", "call_VR" ]
 
 from os.path import basename, dirname, join as joinpath
@@ -63,7 +61,4 @@
 
 def isstr(s):
-    """
-    Return True if *s* is a string-like object.
-    """
     try: s + ''
     except: return False
@@ -104,5 +99,5 @@
     # source = open(info['name']+'.cl','r').read()
 
-    if (platform == "dll"
+    if (platform=="dll"
             or not HAVE_OPENCL
             or not kernelcl.environment().has_type(dtype)):
@@ -133,5 +128,5 @@
     width = pars.get(name+'_pd', 0.0)
     nsigma = pars.get(name+'_pd_nsigma', 3.0)
-    value, weight = weights.get_weights(
+    value,weight = weights.get_weights(
         disperser, npts, width, nsigma, value, limits, relative)
     return value, weight / np.sum(weight)
@@ -200,5 +195,5 @@
                     for name in info['partype']['volume']]
         value, weight = dispersion_mesh(vol_pars)
-        whole, part = VR(*value)
+        whole,part = VR(*value)
         return np.sum(weight*part)/np.sum(weight*whole)
 
Index: sasmodels/data.py
===================================================================
--- sasmodels/data.py	(revision d15a908796935025260dc12a6d16932a079552d3)
+++ sasmodels/data.py	(revision 013adb7ddc1b713ccb5eff4052d37a298cc674e4)
@@ -284,5 +284,5 @@
                 mdata[mdata <= 0] = masked
             plt.errorbar(data.x/10, scale*mdata, yerr=data.dy, fmt='.')
-            all_positive = all_positive and (mdata > 0).all()
+            all_positive = all_positive and (mdata>0).all()
             some_present = some_present or (mdata.count() > 0)
 
@@ -292,7 +292,7 @@
             mtheory[~np.isfinite(mtheory)] = masked
             if view is 'log':
-                mtheory[mtheory <= 0] = masked
+                mtheory[mtheory<=0] = masked
             plt.plot(data.x/10, scale*mtheory, '-', hold=True)
-            all_positive = all_positive and (mtheory > 0).all()
+            all_positive = all_positive and (mtheory>0).all()
             some_present = some_present or (mtheory.count() > 0)
 
@@ -396,5 +396,5 @@
         plt.title('theory')
         h = plt.colorbar()
-        h.set_label(r'$\log_{10}I(q)$' if view == 'log'
+        h.set_label(r'$\log_{10}I(q)$' if view=='log'
                     else r'$q^4 I(q)$' if view == 'q4'
                     else '$I(q)$')
@@ -411,5 +411,5 @@
         plt.title('residuals')
         h = plt.colorbar()
-        h.set_label(r'$\Delta I(q)$')
+        h.set_label('$\Delta I(q)$')
 
 
Index: sasmodels/resolution.py
===================================================================
--- sasmodels/resolution.py	(revision fdc538a09f8fc03aca7f7563f14b6fc5ca0d178f)
+++ sasmodels/resolution.py	(revision 5258859dcf2d6cbd1ba2ea63ce85ff792c2056bb)
@@ -135,5 +135,5 @@
     """
     #print("apply shapes", theory.shape, weight_matrix.shape)
-    Iq = np.dot(theory[None, :], weight_matrix)
+    Iq = np.dot(theory[None,:], weight_matrix)
     #print("result shape",Iq.shape)
     return Iq.flatten()
@@ -153,8 +153,8 @@
     # neither trapezoid nor Simpson's rule improved the accuracy.
     edges = bin_edges(q_calc)
-    edges[edges < 0.0] = 0.0 # clip edges below zero
-    G = erf((edges[:, None] - q[None, :]) / (sqrt(2.0)*q_width)[None, :])
+    edges[edges<0.0] = 0.0 # clip edges below zero
+    G = erf( (edges[:,None] - q[None,:]) / (sqrt(2.0)*q_width)[None,:] )
     weights = G[1:] - G[:-1]
-    weights /= np.sum(weights, axis=0)[None, :]
+    weights /= np.sum(weights, axis=0)[None,:]
     return weights
 
@@ -287,5 +287,5 @@
     # The current algorithm is a midpoint rectangle rule.
     q_edges = bin_edges(q_calc) # Note: requires q > 0
-    q_edges[q_edges < 0.0] = 0.0 # clip edges below zero
+    q_edges[q_edges<0.0] = 0.0 # clip edges below zero
     weights = np.zeros((len(q), len(q_calc)), 'd')
 
@@ -306,10 +306,10 @@
             #print(qi - h, qi + h)
             #print(in_x + abs_x)
-            weights[i, :] = (in_x + abs_x) * np.diff(q_edges) / (2*h)
+            weights[i,:] = (in_x + abs_x) * np.diff(q_edges) / (2*h)
         else:
             L = n_height
             for k in range(-L, L+1):
-                weights[i, :] += _q_perp_weights(q_edges, qi+k*h/L, w)
-            weights[i, :] /= 2*L + 1
+                weights[i,:] += _q_perp_weights(q_edges, qi+k*h/L, w)
+            weights[i,:] /= 2*L + 1
 
     return weights.T
@@ -358,5 +358,5 @@
     log-scaled data.
     """
-    if len(x) < 2 or (np.diff(x) < 0).any():
+    if len(x) < 2 or (np.diff(x)<0).any():
         raise ValueError("Expected bins to be an increasing set")
     edges = np.hstack([
@@ -373,13 +373,13 @@
     """
     step = np.diff(q)
-    index = step > max_step
+    index = step>max_step
     if np.any(index):
         inserts = []
-        for q_i, step_i in zip(q[:-1][index], step[index]):
+        for q_i,step_i in zip(q[:-1][index],step[index]):
             n = np.ceil(step_i/max_step)
-            inserts.extend(q_i + np.arange(1, n)*(step_i/n))
+            inserts.extend(q_i + np.arange(1,n)*(step_i/n))
         # Extend a couple of fringes beyond the end of the data
-        inserts.extend(q[-1] + np.arange(1, 8)*max_step)
-        q_calc = np.sort(np.hstack((q, inserts)))
+        inserts.extend(q[-1] + np.arange(1,8)*max_step)
+        q_calc = np.sort(np.hstack((q,inserts)))
     else:
         q_calc = q
@@ -399,10 +399,10 @@
     if q_min < 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
+        n_low = np.ceil((q[0]-q_min) / (q[1]-q[0])) if q[1]>q[0] else 15
         q_low = np.linspace(q_min, q[0], n_low+1)[:-1]
     else:
         q_low = []
     if q_max > q[-1]:
-        n_high = np.ceil((q_max-q[-1]) / (q[-1]-q[-2])) if q[-1] > q[-2] else 15
+        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:]
     else:
@@ -452,5 +452,5 @@
         if q_min < 0: q_min = q[0]*MIN_Q_SCALE_FOR_NEGATIVE_Q_EXTRAPOLATION
         n_low = log_delta_q * (log(q[0])-log(q_min))
-        q_low = np.logspace(log10(q_min), log10(q[0]), np.ceil(n_low)+1)[:-1]
+        q_low  = np.logspace(log10(q_min), log10(q[0]), np.ceil(n_low)+1)[:-1]
     else:
         q_low = []
@@ -489,5 +489,5 @@
     that make it slow to evaluate but give it good accuracy.
     """
-    from scipy.integrate import romberg
+    from scipy.integrate import romberg, dblquad
 
     if any(k not in form.info['defaults'] for k in pars.keys()):
@@ -520,15 +520,16 @@
             result[i] = r/(2*h)
         else:
-            w_grid = np.linspace(0, w, 21)[None, :]
-            h_grid = np.linspace(-h, h, 23)[:, None]
+            w_grid = np.linspace(0, w, 21)[None,:]
+            h_grid = np.linspace(-h, h, 23)[:,None]
             u = sqrt((qi+h_grid)**2 + w_grid**2)
             Iu = np.interp(u, q_calc, Iq)
             #print(np.trapz(Iu, w_grid, axis=1))
-            Is = np.trapz(np.trapz(Iu, w_grid, axis=1), h_grid[:, 0])
+            Is = np.trapz(np.trapz(Iu, w_grid, axis=1), h_grid[:,0])
             result[i] = Is / (2*h*w)
-            # from scipy.integrate import dblquad
-            # r, err = dblquad(_int_wh, -h, h, lambda h: 0., lambda h: w,
-            #                  args=(qi,))
-            # result[i] = r/(w*2*h)
+            """
+            r, err = dblquad(_int_wh, -h, h, lambda h: 0., lambda h: w,
+                             args=(qi,))
+            result[i] = r/(w*2*h)
+            """
 
     # r should be [float, ...], but it is [array([float]), array([float]),...]
@@ -552,7 +553,7 @@
 
     _fn = lambda q, q0, dq: eval_form(q, form, pars)*gaussian(q, q0, dq)
-    r = [romberg(_fn, max(qi-nsigma*dqi, 1e-10*q[0]), qi+nsigma*dqi,
-                 args=(qi, dqi), divmax=100, vec_func=True, tol=0, rtol=1e-8)
-         for qi, dqi in zip(q, q_width)]
+    r = [romberg(_fn, max(qi-nsigma*dqi,1e-10*q[0]), qi+nsigma*dqi, args=(qi, dqi),
+                 divmax=100, vec_func=True, tol=0, rtol=1e-8)
+         for qi,dqi in zip(q,q_width)]
     return np.asarray(r).flatten()
 
@@ -594,8 +595,6 @@
         theory = self.Iq(resolution.q_calc)
         output = resolution.apply(theory)
-        answer = [
-            9.0618, 8.6402, 8.1187, 7.1392, 6.1528,
-            5.5555, 4.5584, 3.5606, 2.5623, 2.0000,
-            ]
+        answer = [ 9.0618, 8.6402, 8.1187, 7.1392, 6.1528,
+                   5.5555, 4.5584, 3.5606, 2.5623, 2.0000 ]
         np.testing.assert_allclose(output, answer, atol=1e-4)
 
@@ -609,8 +608,6 @@
         theory = 1000*self.Iq(resolution.q_calc**4)
         output = resolution.apply(theory)
-        answer = [
-            8.85785, 8.43012, 7.92687, 6.94566, 6.03660,
-            5.40363, 4.40655, 3.40880, 2.41058, 2.00000,
-            ]
+        answer = [ 8.85785, 8.43012, 7.92687, 6.94566, 6.03660,
+                   5.40363, 4.40655, 3.40880, 2.41058, 2.00000 ]
         np.testing.assert_allclose(output, answer, atol=1e-4)
 
@@ -623,7 +620,5 @@
         theory = self.Iq(resolution.q_calc)
         output = resolution.apply(theory)
-        answer = [
-            11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0,
-            ]
+        answer = [ 11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0 ]
         np.testing.assert_allclose(output, answer, atol=1e-4)
 
@@ -637,7 +632,5 @@
         theory = self.Iq(resolution.q_calc)
         output = resolution.apply(theory)
-        answer = [
-            11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0,
-            ]
+        answer = [ 11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0 ]
         np.testing.assert_allclose(output, answer, atol=1e-4)
 
@@ -659,9 +652,7 @@
         theory = 12.0-1000.0*resolution.q_calc
         output = resolution.apply(theory)
-        answer = [
-            10.44785079, 9.84991299, 8.98101708,
-            7.99906585, 6.99998311, 6.00001689,
-            5.00093415, 4.01898292, 3.15008701, 2.55214921,
-            ]
+        answer = [ 10.44785079, 9.84991299, 8.98101708,
+                  7.99906585, 6.99998311, 6.00001689,
+                  5.00093415, 4.01898292, 3.15008701, 2.55214921]
         np.testing.assert_allclose(output, answer, atol=1e-8)
 
@@ -792,5 +783,5 @@
             }
         form = load_model('ellipsoid', dtype='double')
-        q = np.logspace(log10(4e-5), log10(2.5e-2), 68)
+        q = np.logspace(log10(4e-5),log10(2.5e-2), 68)
         width, height = 0.117, 0.
         resolution = Slit1D(q, width=width, height=height)
@@ -1080,5 +1071,5 @@
     #h = 0.0277790
     resolution = Slit1D(q, w, h)
-    _eval_demo_1d(resolution, title="(%g,%g) Slit Resolution"%(w, h))
+    _eval_demo_1d(resolution, title="(%g,%g) Slit Resolution"%(w,h))
 
 def demo():
Index: sasmodels/sesans.py
===================================================================
--- sasmodels/sesans.py	(revision 384d114757de20d9c246bd823485afb6835ec978)
+++ sasmodels/sesans.py	(revision 3c56da87f0995bdde77f4183b68e77b81584dde0)
@@ -12,25 +12,40 @@
 import numpy as np
 from numpy import pi, exp
+
 from scipy.special import jv as besselj
 
-def make_q(q_max, Rmax):
-    r"""
-    Return a $q$ vector suitable for SESANS covering from $2\pi/ (10 R_{\max})$
-    to $q_max$.
-    """
+def make_q(q_zmax, Rmax):
     q_min = dq = 0.1 * 2*pi / Rmax
-    return np.arange(q_min, q_max, dq)
+    #q_min = 0.00003
+    return np.arange(q_min, q_zmax, dq)
+
+# TODO: dead code; for now the call to the hankel transform happens in BumpsModel
+class SesansCalculator:
+    def __init__(self, kernel, q_zmax, Rmax, SElength, wavelength, thickness):
+        self._set_kernel(kernel, q_zmax, Rmax)
+        self.SElength = SElength
+        self.wavelength = wavelength
+        self.thickness = thickness
+
+    def _set_kernel(self, kernel, q_zmax, Rmax):
+        kernel_input = kernel.make_input([make_q(q_zmax, Rmax)])
+        self.sans_calculator = kernel(kernel_input)
+
+    def __call__(self, pars, pd_pars, cutoff=1e-5):
+        Iq = self.sans_calculator(pars, pd_pars, cutoff)
+        P = hankel(self.SElength, self.wavelength, self.thickness, self.q, Iq)
+        self.Iq = Iq
+        return P
 
 def hankel(SElength, wavelength, thickness, q, Iq):
-    r"""
+    """
     Compute the expected SESANS polarization for a given SANS pattern.
 
-    Uses the hankel transform followed by the exponential.  The values for *zz*
-    (or spin echo length, or delta), wavelength and sample thickness should
-    come from the dataset.  $q$ should be chosen such that the oscillations
-    in $I(q)$ are well sampled (e.g., $5 \cdot 2 \pi/d_{\max}$).
+    Uses the hankel transform followed by the exponential.  The values
+    for zz (or spin echo length, or delta), wavelength and sample thickness
+    information should come from the dataset.  *q* should be chosen such
+    that the oscillations in *I(q)* are well sampled (e.g., 5*2*pi/d_max).
 
-    *SElength* [A] is the set of $z$ points at which to compute the
-    Hankel transform
+    *SElength* [A] is the set of z points at which to compute the hankel transform
 
     *wavelength* [m]  is the wavelength of each individual point *zz*
@@ -38,20 +53,15 @@
     *thickness* [cm] is the sample thickness.
 
-    *q* [A$^{-1}$] is the set of $q$ points at which the model has been
-    computed. These should be equally spaced.
+    *q* [A^{-1}] is the set of q points at which the model has been computed.
+    These should be equally spaced.
 
-    *I* [cm$^{-1}$] is the value of the SANS model at *q*
+    *I* [cm^{-1}] is the value of the SANS model at *q*
     """
     G = np.zeros(len(SElength), 'd')
     for i in range(len(SElength)):
-        integr = besselj(0, q*SElength[i])*Iq*q
+        integr = besselj(0,q*SElength[i])*Iq*q
         G[i] = np.sum(integr)
-
-    # [m^-1] step size in q, needed for integration
-    dq=(q[1]-q[0])*1e10
-
-    # integration step, convert q into [m**-1] and 2 pi circle integration
-    G *= dq*1e10*2*pi
-
+    dq=(q[1]-q[0])*1e10   # [m^-1] step size in q, needed for integration
+    G *= dq*1e10*2*pi # integr step, conver q into [m**-1] and 2 pi circle integr
     P = exp(thickness*wavelength**2/(4*pi**2)*(G-G[0]))
 
