Index: extra/pylint.rc
===================================================================
--- extra/pylint.rc	(revision 37a725200a1edb6dbd7a029a236c3ea2cfc7522c)
+++ extra/pylint.rc	(revision d4666ca59132213b621b9cd27fda96521894152e)
@@ -7,4 +7,5 @@
 # pygtk.require().
 #init-hook='import os, sys; sys.path.extend([os.getcwd(),os.path.join(os.getcwd(),"extra")])'
+init-hook='import sys; sys.path.append('extra')
 
 # Profiled execution.
@@ -64,4 +65,5 @@
     star-args,
     unbalanced-tuple-unpacking,
+    locally-enabled,
     locally-disabled,
 
@@ -107,10 +109,11 @@
 good-names=_,
     input,
-    i,j,k,n,x,y,z,
+    h,i,j,k,n,w,x,y,z,
     q,qx,qy,qz,
     dt,dx,dy,dz,id,
-    Iq,dIq,Qx,Qy,Qz,
+    Iq,dIq,Qx,Qy,Qz,Iqxy,
     p,
     ER, call_ER, VR, call_VR,
+    Rmax, SElength,
 
 # Bad variable names which should always be refused, separated by a comma
@@ -335,5 +338,5 @@
 
 # Maximum number of arguments for function / method
-max-args=5
+max-args=15
 
 # Argument names that match this expression will be ignored. Default to name
@@ -342,5 +345,5 @@
 
 # Maximum number of locals for function / method body
-max-locals=15
+max-locals=25
 
 # Maximum number of return / yield for function / method body
Index: sasmodels/bumps_model.py
===================================================================
--- sasmodels/bumps_model.py	(revision 37a725200a1edb6dbd7a029a236c3ea2cfc7522c)
+++ sasmodels/bumps_model.py	(revision 190fc2b167d3d89ca9c3b7d3128bdf8bde9f67d0)
@@ -12,14 +12,14 @@
 """
 
+import warnings
+
+import numpy as np
+
+from .data import plot_theory
+from .direct_model import DataMixin
+
 __all__ = [
     "Model", "Experiment",
     ]
-
-import warnings
-
-import numpy as np
-
-from .data import plot_theory
-from .direct_model import DataMixin
 
 # CRUFT: old style bumps wrapper which doesn't separate data and model
Index: sasmodels/compare.py
===================================================================
--- sasmodels/compare.py	(revision d15a908796935025260dc12a6d16932a079552d3)
+++ sasmodels/compare.py	(revision 190fc2b167d3d89ca9c3b7d3128bdf8bde9f67d0)
@@ -28,4 +28,20 @@
 
 from __future__ import print_function
+
+import sys
+import math
+from os.path import basename, dirname, join as joinpath
+import glob
+import datetime
+import traceback
+
+import numpy as np
+
+from . import core
+from . import kerneldll
+from . import generate
+from .data import plot_theory, empty_data1D, empty_data2D
+from .direct_model import DirectModel
+from .convert import revert_model, constrain_new_to_old
 
 USAGE = """
@@ -80,28 +96,8 @@
            + USAGE)
 
-
-
-import sys
-import math
-from os.path import basename, dirname, join as joinpath
-import glob
-import datetime
-import traceback
-
-import numpy as np
-
+kerneldll.ALLOW_SINGLE_PRECISION_DLLS = True
+
+# List of available models
 ROOT = dirname(__file__)
-sys.path.insert(0, ROOT)  # Make sure sasmodels is first on the path
-
-
-from . import core
-from . import kerneldll
-from . import generate
-from .data import plot_theory, empty_data1D, empty_data2D
-from .direct_model import DirectModel
-from .convert import revert_model, constrain_new_to_old
-kerneldll.ALLOW_SINGLE_PRECISION_DLLS = True
-
-# List of available models
 MODELS = [basename(f)[:-3]
           for f in sorted(glob.glob(joinpath(ROOT, "models", "[a-zA-Z]*.py")))]
@@ -117,4 +113,72 @@
         return dt.total_seconds()
 
+
+class push_seed(object):
+    """
+    Set the seed value for the random number generator.
+
+    When used in a with statement, the random number generator state is
+    restored after the with statement is complete.
+
+    :Parameters:
+
+    *seed* : int or array_like, optional
+        Seed for RandomState
+
+    :Example:
+
+    Seed can be used directly to set the seed::
+
+        >>> from numpy.random import randint
+        >>> push_seed(24)
+        <...push_seed object at...>
+        >>> print(randint(0,1000000,3))
+        [242082    899 211136]
+
+    Seed can also be used in a with statement, which sets the random
+    number generator state for the enclosed computations and restores
+    it to the previous state on completion::
+
+        >>> with push_seed(24):
+        ...    print(randint(0,1000000,3))
+        [242082    899 211136]
+
+    Using nested contexts, we can demonstrate that state is indeed
+    restored after the block completes::
+
+        >>> with push_seed(24):
+        ...    print(randint(0,1000000))
+        ...    with push_seed(24):
+        ...        print(randint(0,1000000,3))
+        ...    print(randint(0,1000000))
+        242082
+        [242082    899 211136]
+        899
+
+    The restore step is protected against exceptions in the block::
+
+        >>> with push_seed(24):
+        ...    print(randint(0,1000000))
+        ...    try:
+        ...        with push_seed(24):
+        ...            print(randint(0,1000000,3))
+        ...            raise Exception()
+        ...    except:
+        ...        print("Exception raised")
+        ...    print(randint(0,1000000))
+        242082
+        [242082    899 211136]
+        Exception raised
+        899
+    """
+    def __init__(self, seed=None):
+        self._state = np.random.get_state()
+        np.random.seed(seed)
+
+    def __enter__(self):
+        return None
+
+    def __exit__(self, *args):
+        np.random.set_state(self._state)
 
 def tic():
@@ -179,4 +243,5 @@
         return [0, (2*v if v > 0 else 1)]
 
+
 def _randomize_one(p, v):
     """
@@ -188,4 +253,5 @@
         return np.random.uniform(*parameter_range(p, v))
 
+
 def randomize_pars(pars, seed=None):
     """
@@ -197,8 +263,8 @@
     :func:`constrain_pars` needs to be called afterward..
     """
-    np.random.seed(seed)
-    # Note: the sort guarantees order `of calls to random number generator
-    pars = dict((p, _randomize_one(p, v))
-                for p, v in sorted(pars.items()))
+    with push_seed(seed):
+        # Note: the sort guarantees order `of calls to random number generator
+        pars = dict((p, _randomize_one(p, v))
+                    for p, v in sorted(pars.items()))
     return pars
 
@@ -463,10 +529,10 @@
     if Nbase > 0:
         if Ncomp > 0: plt.subplot(131)
-        plot_theory(data, base_value, view=view, plot_data=False, limits=limits)
+        plot_theory(data, base_value, view=view, use_data=False, limits=limits)
         plt.title("%s t=%.1f ms"%(base.engine, base_time))
         #cbar_title = "log I"
     if Ncomp > 0:
         if Nbase > 0: plt.subplot(132)
-        plot_theory(data, comp_value, view=view, plot_data=False, limits=limits)
+        plot_theory(data, comp_value, view=view, use_data=False, limits=limits)
         plt.title("%s t=%.1f ms"%(comp.engine, comp_time))
         #cbar_title = "log I"
@@ -478,5 +544,5 @@
             err, errstr, errview = abs(relerr), "rel err", "log"
         #err,errstr = base/comp,"ratio"
-        plot_theory(data, None, resid=err, view=errview, plot_data=False)
+        plot_theory(data, None, resid=err, view=errview, use_data=False)
         plt.title("max %s = %.3g"%(errstr, max(abs(err))))
         #cbar_title = errstr if errview=="linear" else "log "+errstr
Index: sasmodels/compare_many.py
===================================================================
--- sasmodels/compare_many.py	(revision d15a908796935025260dc12a6d16932a079552d3)
+++ sasmodels/compare_many.py	(revision 4f2478e77cdc4e495cf403a9aeb9c6544bcdf468)
@@ -261,3 +261,5 @@
 
 if __name__ == "__main__":
+    #from .compare import push_seed
+    #with push_seed(1): main()
     main()
Index: sasmodels/core.py
===================================================================
--- sasmodels/core.py	(revision d15a908796935025260dc12a6d16932a079552d3)
+++ sasmodels/core.py	(revision 190fc2b167d3d89ca9c3b7d3128bdf8bde9f67d0)
@@ -2,8 +2,4 @@
 Core model handling routines.
 """
-__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
@@ -24,4 +20,8 @@
     HAVE_OPENCL = False
 
+__all__ = [
+    "list_models", "load_model_definition", "precompile_dll",
+    "load_model", "make_kernel", "call_kernel", "call_ER", "call_VR",
+]
 
 def list_models():
Index: sasmodels/data.py
===================================================================
--- sasmodels/data.py	(revision d15a908796935025260dc12a6d16932a079552d3)
+++ sasmodels/data.py	(revision 69ec80f641366fc1c80fed9bcfff0d283905e6ef)
@@ -90,4 +90,10 @@
         self.x, self.y, self.dx, self.dy = x, y, dx, dy
         self.dxl = None
+        self.filename = None
+        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
+        self._xaxis, self._xunit = "x", ""
+        self._yaxis, self._yunit = "y", ""
 
     def xaxis(self, label, unit):
@@ -108,7 +114,23 @@
 
 class Data2D(object):
-    def __init__(self):
+    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.q_data = np.sqrt(x**2 + y**2)
+        self.qmin = 1e-16
+        self.qmax = np.inf
         self.detector = []
         self.source = Source()
+        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, self._xunit = "x", ""
+        self._yaxis, self._yunit = "y", ""
+        self._zaxis, self._zunit = "z", ""
+        self.x_bins, self.y_bins = None, None
 
     def xaxis(self, label, unit):
@@ -139,9 +161,18 @@
 
 class Detector(object):
+    """
+    Detector attributes.
+    """
+    def __init__(self, pixel_size=(None, None), distance=None):
+        self.pixel_size = Vector(*pixel_size)
+        self.distance = distance
+
+class Source(object):
+    """
+    Beam attributes.
+    """
     def __init__(self):
-        self.pixel_size = Vector()
-
-class Source(object):
-    pass
+        self.wavelength = np.NaN
+        self.wavelength_unit = "A"
 
 
@@ -158,6 +189,4 @@
     data = Data1D(q, Iq, dx=resolution * q, dy=dIq)
     data.filename = "fake data"
-    data.qmin, data.qmax = q.min(), q.max()
-    data.mask = np.zeros(len(q), dtype='bool')
     return data
 
@@ -173,21 +202,9 @@
     if qy is None:
         qy = qx
+    # 5% dQ/Q resolution
     Qx, Qy = np.meshgrid(qx, qy)
     Qx, Qy = Qx.flatten(), Qy.flatten()
     Iq = 100 * np.ones_like(Qx)
     dIq = np.sqrt(Iq)
-    mask = np.ones(len(Iq), dtype='bool')
-
-    data = Data2D()
-    data.filename = "fake data"
-    data.qx_data = Qx
-    data.qy_data = Qy
-    data.data = Iq
-    data.err_data = dIq
-    data.mask = mask
-    data.qmin = 1e-16
-    data.qmax = np.inf
-
-    # 5% dQ/Q resolution
     if resolution != 0:
         # https://www.ncnr.nist.gov/staff/hammouda/distance_learning/chapter_15.pdf
@@ -197,24 +214,19 @@
         # radial (which instead it should be inverse).
         Q = np.sqrt(Qx**2 + Qy**2)
-        data.dqx_data = resolution * Q
-        data.dqy_data = resolution * Q
+        dqx = resolution * Q
+        dqy = resolution * Q
     else:
-        data.dqx_data = data.dqy_data = None
-
-    detector = Detector()
-    detector.pixel_size.x = 5 # mm
-    detector.pixel_size.y = 5 # mm
-    detector.distance = 4 # m
-    data.detector.append(detector)
+        dqx = dqy = None
+
+    data = Data2D(x=Qx, y=Qy, z=Iq, dx=dqx, dy=dqy, dz=dIq)
     data.x_bins = qx
     data.y_bins = qy
+    data.filename = "fake data"
+
+    # pixel_size in mm, distance in m
+    detector = Detector(pixel_size=(5, 5), distance=4)
+    data.detector.append(detector)
     data.source.wavelength = 5 # angstroms
     data.source.wavelength_unit = "A"
-    data.Q_unit = "1/A"
-    data.I_unit = "1/cm"
-    data.q_data = np.sqrt(Qx ** 2 + Qy ** 2)
-    data.xaxis("Q_x", "A^{-1}")
-    data.yaxis("Q_y", "A^{-1}")
-    data.zaxis("Intensity", r"\text{cm}^{-1}")
     return data
 
@@ -228,19 +240,19 @@
     # do not repeat.
     if hasattr(data, 'lam'):
-        _plot_result_sesans(data, None, None, plot_data=True, limits=limits)
+        _plot_result_sesans(data, None, None, use_data=True, limits=limits)
     elif hasattr(data, 'qx_data'):
-        _plot_result2D(data, None, None, view, plot_data=True, limits=limits)
+        _plot_result2D(data, None, None, view, use_data=True, limits=limits)
     else:
-        _plot_result1D(data, None, None, view, plot_data=True, limits=limits)
+        _plot_result1D(data, None, None, view, use_data=True, limits=limits)
 
 
 def plot_theory(data, theory, resid=None, view='log',
-                plot_data=True, limits=None):
+                use_data=True, limits=None):
     if hasattr(data, 'lam'):
-        _plot_result_sesans(data, theory, resid, plot_data=True, limits=limits)
+        _plot_result_sesans(data, theory, resid, use_data=True, limits=limits)
     elif hasattr(data, 'qx_data'):
-        _plot_result2D(data, theory, resid, view, plot_data, limits=limits)
+        _plot_result2D(data, theory, resid, view, use_data, limits=limits)
     else:
-        _plot_result1D(data, theory, resid, view, plot_data, limits=limits)
+        _plot_result1D(data, theory, resid, view, use_data, limits=limits)
 
 
@@ -251,5 +263,4 @@
         except:
             traceback.print_exc()
-            pass
 
     return wrapper
@@ -257,5 +268,5 @@
 
 @protect
-def _plot_result1D(data, theory, resid, view, plot_data, limits=None):
+def _plot_result1D(data, theory, resid, view, use_data, limits=None):
     """
     Plot the data and residuals for 1D data.
@@ -264,19 +275,16 @@
     from numpy.ma import masked_array, masked
 
-    plot_theory = theory is not None
-    plot_resid = resid is not None
-
-    if data.y is None:
-        plot_data = False
+    use_data = use_data and data.y is not None
+    use_theory = theory is not None
+    use_resid = resid is not None
+    num_plots = (use_data or use_theory) + use_resid
+
     scale = data.x**4 if view == 'q4' else 1.0
 
-    if plot_data or plot_theory:
-        if plot_resid:
-            plt.subplot(121)
-
+    if use_data or use_theory:
         #print(vmin, vmax)
         all_positive = True
         some_present = False
-        if plot_data:
+        if use_data:
             mdata = masked_array(data.y, data.mask.copy())
             mdata[~np.isfinite(mdata)] = masked
@@ -288,5 +296,5 @@
 
 
-        if plot_theory:
+        if use_theory:
             mtheory = masked_array(theory, data.mask.copy())
             mtheory[~np.isfinite(mtheory)] = masked
@@ -299,4 +307,7 @@
         if limits is not None:
             plt.ylim(*limits)
+
+        if num_plots > 1:
+            plt.subplot(1, num_plots, 1)
         plt.xscale('linear' if not some_present else view)
         plt.yscale('linear'
@@ -306,11 +317,11 @@
         plt.ylabel('$I(q)$')
 
-    if plot_resid:
-        if plot_data or plot_theory:
-            plt.subplot(122)
-
+    if use_resid:
         mresid = masked_array(resid, data.mask.copy())
         mresid[~np.isfinite(mresid)] = masked
         some_present = (mresid.count() > 0)
+
+        if num_plots > 1:
+            plt.subplot(1, num_plots, (use_data or use_theory) + 1)
         plt.plot(data.x/10, mresid, '-')
         plt.xlabel("$q$/nm$^{-1}$")
@@ -320,15 +331,15 @@
 
 @protect
-def _plot_result_sesans(data, theory, resid, plot_data, limits=None):
+def _plot_result_sesans(data, theory, resid, use_data, limits=None):
     import matplotlib.pyplot as plt
-    if data.y is None:
-        plot_data = False
-    plot_theory = theory is not None
-    plot_resid = resid is not None
-
-    if plot_data or plot_theory:
-        if plot_resid:
-            plt.subplot(121)
-        if plot_data:
+    use_data = use_data and data.y is not None
+    use_theory = theory is not None
+    use_resid = resid is not None
+    num_plots = (use_data or use_theory) + use_resid
+
+    if use_data or use_theory:
+        if num_plots > 1:
+            plt.subplot(1, num_plots, 1)
+        if use_data:
             plt.errorbar(data.x, data.y, yerr=data.dy)
         if theory is not None:
@@ -340,7 +351,6 @@
 
     if resid is not None:
-        if plot_data or plot_theory:
-            plt.subplot(122)
-
+        if num_plots > 1:
+            plt.subplot(1, num_plots, (use_data or use_theory) + 1)
         plt.plot(data.x, resid, 'x')
         plt.xlabel('spin echo length (nm)')
@@ -349,36 +359,36 @@
 
 @protect
-def _plot_result2D(data, theory, resid, view, plot_data, limits=None):
+def _plot_result2D(data, theory, resid, view, use_data, limits=None):
     """
     Plot the data and residuals for 2D data.
     """
     import matplotlib.pyplot as plt
-    if data.data is None:
-        plot_data = False
-    plot_theory = theory is not None
-    plot_resid = resid is not None
+    use_data = use_data and data.data is not None
+    use_theory = theory is not None
+    use_resid = resid is not None
+    num_plots = use_data + use_theory + use_resid
 
     # Put theory and data on a common colormap scale
-    if limits is None:
-        vmin, vmax = np.inf, -np.inf
-        if plot_data:
-            target = data.data[~data.mask]
-            datamin = target[target>0].min() if view == 'log' else target.min()
-            datamax = target.max()
-            vmin = min(vmin, datamin)
-            vmax = max(vmax, datamax)
-        if plot_theory:
-            theorymin = theory[theory>0].min() if view=='log' else theory.min()
-            theorymax = theory.max()
-            vmin = min(vmin, theorymin)
-            vmax = max(vmax, theorymax)
-    else:
+    vmin, vmax = np.inf, -np.inf
+    if use_data:
+        target = data.data[~data.mask]
+        datamin = target[target > 0].min() if view == 'log' else target.min()
+        datamax = target.max()
+        vmin = min(vmin, datamin)
+        vmax = max(vmax, datamax)
+    if use_theory:
+        theorymin = theory[theory > 0].min() if view == 'log' else theory.min()
+        theorymax = theory.max()
+        vmin = min(vmin, theorymin)
+        vmax = max(vmax, theorymax)
+
+    # Override data limits from the caller
+    if limits is not None:
         vmin, vmax = limits
 
-    if plot_data:
-        if plot_theory and plot_resid:
-            plt.subplot(131)
-        elif plot_theory or plot_resid:
-            plt.subplot(121)
+    # Plot data
+    if use_data:
+        if num_plots > 1:
+            plt.subplot(1, num_plots, 1)
         _plot_2d_signal(data, target, view=view, vmin=vmin, vmax=vmax)
         plt.title('data')
@@ -386,11 +396,8 @@
         h.set_label('$I(q)$')
 
-    if plot_theory:
-        if plot_data and plot_resid:
-            plt.subplot(132)
-        elif plot_data:
-            plt.subplot(122)
-        elif plot_resid:
-            plt.subplot(121)
+    # plot theory
+    if use_theory:
+        if num_plots > 1:
+            plt.subplot(1, num_plots, use_data+1)
         _plot_2d_signal(data, theory, view=view, vmin=vmin, vmax=vmax)
         plt.title('theory')
@@ -400,12 +407,8 @@
                     else '$I(q)$')
 
-    #if plot_data or plot_theory:
-    #    plt.colorbar()
-
-    if plot_resid:
-        if plot_data and plot_theory:
-            plt.subplot(133)
-        elif plot_data or plot_theory:
-            plt.subplot(122)
+    # plot resid
+    if use_resid:
+        if num_plots > 1:
+            plt.subplot(1, num_plots, use_data+use_theory+1)
         _plot_2d_signal(data, resid, view='linear')
         plt.title('residuals')
Index: sasmodels/generate.py
===================================================================
--- sasmodels/generate.py	(revision e66c9f9bddfbf0bb75d5d64484195eaf6802ca82)
+++ sasmodels/generate.py	(revision 190fc2b167d3d89ca9c3b7d3128bdf8bde9f67d0)
@@ -197,6 +197,4 @@
 # TODO: identify model files which have changed since loading and reload them.
 
-__all__ = ["make", "doc", "sources", "convert_type"]
-
 import sys
 from os.path import abspath, dirname, join as joinpath, exists, basename, \
@@ -206,4 +204,7 @@
 
 import numpy as np
+
+__all__ = ["make", "doc", "sources", "convert_type"]
+
 C_KERNEL_TEMPLATE_PATH = joinpath(dirname(__file__), 'kernel_template.c')
 
@@ -216,5 +217,4 @@
     F128 = None
 
-
 # Scale and background, which are parameters common to every form factor
 COMMON_PARAMETERS = [
@@ -222,5 +222,4 @@
     ["background", "1/cm", 0, [0, np.inf], "", "Source background"],
     ]
-
 
 # Conversion from units defined in the parameter table for each model
@@ -266,4 +265,7 @@
 
 def format_units(units):
+    """
+    Convert units into ReStructured Text format.
+    """
     return "string" if isinstance(units, list) else RST_UNITS.get(units, units)
 
@@ -335,4 +337,7 @@
     """
     Convert code from double precision to the desired type.
+
+    Floating point constants are tagged with 'f' for single precision or 'L'
+    for long double precision.
     """
     if dtype == F16:
@@ -350,4 +355,8 @@
 
 def _convert_type(source, type_name, constant_flag):
+    """
+    Replace 'double' with *type_name* in *source*, tagging floating point
+    constants with *constant_flag*.
+    """
     # Convert double keyword to float/long double/half.
     # Accept an 'n' # parameter for vector # values, where n is 2, 4, 8 or 16.
@@ -625,4 +634,7 @@
                             %re.escape(string.punctuation))
 def _convert_section_titles_to_boldface(lines):
+    """
+    Do the actual work of identifying and converting section headings.
+    """
     prior = None
     for line in lines:
@@ -642,6 +654,13 @@
         yield prior
 
-def convert_section_titles_to_boldface(string):
-    return "\n".join(_convert_section_titles_to_boldface(string.split('\n')))
+def convert_section_titles_to_boldface(s):
+    """
+    Use explicit bold-face rather than section headings so that the table of
+    contents is not polluted with section names from the model documentation.
+
+    Sections are identified as the title line followed by a line of punctuation
+    at least as long as the title line.
+    """
+    return "\n".join(_convert_section_titles_to_boldface(s.split('\n')))
 
 def doc(kernel_module):
@@ -666,4 +685,7 @@
 
 def demo_time():
+    """
+    Show how long it takes to process a model.
+    """
     from .models import cylinder
     import datetime
@@ -674,4 +696,7 @@
 
 def main():
+    """
+    Program which prints the source produced by the model.
+    """
     if len(sys.argv) <= 1:
         print("usage: python -m sasmodels.generate modelname")
Index: sasmodels/models/lamellarCailleHG.py
===================================================================
--- sasmodels/models/lamellarCailleHG.py	(revision d18f8a8e10d65f0ace6fdc9360abd60dfcd6dcce)
+++ sasmodels/models/lamellarCailleHG.py	(revision d4666ca59132213b621b9cd27fda96521894152e)
@@ -92,25 +92,25 @@
 
 parameters = [
-              #   [ "name", "units", default, [lower, upper], "type",
-              #     "description" ],
-              [ "tail_length", "Ang",  10, [0, inf], "volume",
-                "Tail thickness" ],
-              [ "head_length", "Ang",  2, [0, inf], "volume",
-                "head thickness" ],
-              [ "Nlayers", "",  30, [0, inf], "",
-                "Number of layers" ],
-              [ "spacing", "Ang", 40., [0.0,inf], "volume",
-                "d-spacing of Caille S(Q)" ],
-              [ "Caille_parameter", "", 0.001, [0.0,0.8], "",
-                "Caille parameter" ],
-              [ "sld", "1e-6/Ang^2", 0.4, [-inf,inf], "",
-                "Tail scattering length density" ],
-              [ "head_sld", "1e-6/Ang^2", 2.0, [-inf,inf], "",
-                "Head scattering length density" ],
-              [ "solvent_sld", "1e-6/Ang^2", 6, [-inf,inf], "",
-                "Solvent scattering length density" ],
+    #   [ "name", "units", default, [lower, upper], "type",
+    #     "description" ],
+    ["tail_length", "Ang", 10, [0, inf], "volume",
+     "Tail thickness"],
+    ["head_length", "Ang", 2, [0, inf], "volume",
+     "head thickness"],
+    ["Nlayers", "", 30, [0, inf], "",
+     "Number of layers"],
+    ["spacing", "Ang", 40., [0.0, inf], "volume",
+     "d-spacing of Caille S(Q)"],
+    ["Caille_parameter", "", 0.001, [0.0, 0.8], "",
+     "Caille parameter"],
+    ["sld", "1e-6/Ang^2", 0.4, [-inf, inf], "",
+     "Tail scattering length density"],
+    ["head_sld", "1e-6/Ang^2", 2.0, [-inf, inf], "",
+     "Head scattering length density"],
+    ["solvent_sld", "1e-6/Ang^2", 6, [-inf, inf], "",
+     "Solvent scattering length density"],
     ]
 
-source = [ "lamellarCailleHG_kernel.c"]
+source = ["lamellarCailleHG_kernel.c"]
 
 # No volume normalization despite having a volume parameter
@@ -128,15 +128,17 @@
 
 demo = dict(
-            scale=1, background=0,
-            Nlayers=20,
-            spacing=200., Caille_parameter=0.05,
-            tail_length=15,head_length=10,
-            #sld=-1, head_sld=4.0, solvent_sld=6.0,
-            sld=-1, head_sld=4.1, solvent_sld=6.0,
-            tail_length_pd= 0.1, tail_length_pd_n=20,
-            head_length_pd= 0.05, head_length_pd_n=30,
-            spacing_pd= 0.2, spacing_pd_n=40
-           )
+    scale=1, background=0,
+    Nlayers=20, spacing=200., Caille_parameter=0.05,
+    tail_length=15, head_length=10,
+    #sld=-1, head_sld=4.0, solvent_sld=6.0,
+    sld=-1, head_sld=4.1, solvent_sld=6.0,
+    tail_length_pd=0.1, tail_length_pd_n=20,
+    head_length_pd=0.05, head_length_pd_n=30,
+    spacing_pd=0.2, spacing_pd_n=40,
+    )
 
 oldname = 'LamellarPSHGModel'
-oldpars = dict(tail_length='deltaT',head_length='deltaH',Nlayers='n_plates',Caille_parameter='caille', sld='sld_tail', head_sld='sld_head',solvent_sld='sld_solvent')
+oldpars = dict(
+    tail_length='deltaT', head_length='deltaH', Nlayers='n_plates',
+    Caille_parameter='caille', sld='sld_tail', head_sld='sld_head',
+    solvent_sld='sld_solvent')
Index: sasmodels/models/lib/sph_j1c.c
===================================================================
--- sasmodels/models/lib/sph_j1c.c	(revision 9c461c713e9961727442c844be3e6354d6548c2f)
+++ sasmodels/models/lib/sph_j1c.c	(revision a3e78c304bba2ffdbd8ad1b8da217f22e084346c)
@@ -1,4 +1,4 @@
 /**
-* Spherical Bessel function j1(x)/x
+* Spherical Bessel function 3*j1(x)/x
 *
 * Used for low q to avoid cancellation error.
@@ -8,4 +8,9 @@
 * requires the first 3 terms.  Double precision requires the 4th term.
 * The fifth term is not needed, and is commented out.
+* Taylor expansion:
+*      1.0 + q2*(-3./30. + q2*(3./840.))+ q2*(-3./45360. + q2*(3./3991680.))))
+* Expression returned from Herbie (herbie.uwpise.org/demo):
+*      const double t = ((1. + 3.*q2*q2/5600.) - q2/20.);
+*      return t*t;
 */
 
@@ -18,8 +23,23 @@
     SINCOS(q, sin_q, cos_q);
 
-    const double bessel = (q < 1.e-1) ?
-        1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360.))) // + q2*(3./3991680.)))
+    const double bessel = (q < 0.384038453352533)
+        ? (1.0 + q2*(-3./30. + q2*(3./840.)))
         : 3.0*(sin_q/q - cos_q)/q2;
 
     return bessel;
+
+ /*
+    // Code to test various expressions
+    if (sizeof(q2) > 4) {
+        return 3.0*(sin_q/q - cos_q)/q2;
+    } else if (q < 0.384038453352533) {
+        //const double t = ((1. + 3.*q2*q2/5600.) - q2/20.); return t*t;
+        return 1.0 + q2*q2*(3./840.) - q2*(3./30.);
+        //return 1.0 + q2*(-3./30. + q2*(3./840.));
+        //return 1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360.)));
+        //return 1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360. + q2*(3./3991680.))));
+    } else {
+        return 3.0*(sin_q/q - cos_q)/q2;
+    }
+*/
 }
Index: sasmodels/models/pearl_necklace.py
===================================================================
--- sasmodels/models/pearl_necklace.py	(revision cf404cbf1344577110785a29f1390cbabc7969b2)
+++ sasmodels/models/pearl_necklace.py	(revision 841753cbac7824298effa1aec88a7ee977eb257f)
@@ -1,7 +1,7 @@
 r"""
-This model provides the form factor for a pearl necklace composed of two 
-elements: *N* pearls (homogeneous spheres of radius *R*) freely jointed by *M* 
+This model provides the form factor for a pearl necklace composed of two
+elements: *N* pearls (homogeneous spheres of radius *R*) freely jointed by *M*
 rods (like strings - with a total mass *Mw* = *M* \* *m*\ :sub:`r` + *N* \* *m*\
-:sub:`s`, and the string segment length (or edge separation) *l* 
+:sub:`s`, and the string segment length (or edge separation) *l*
 (= *A* - 2\ *R*)). *A* is the center-to-center pearl separation distance.
 
@@ -13,5 +13,5 @@
 ----------
 
-The output of the scattering intensity function for the PearlNecklaceModel is 
+The output of the scattering intensity function for the PearlNecklaceModel is
 given by (Schweins, 2004)
 
@@ -37,8 +37,8 @@
     \beta(q) &= \frac{\int_{qR}^{q(A-R)}\frac{sin(t)}{t}dt}{ql}
 
-where the mass *m*\ :sub:`i` is (SLD\ :sub:`i` - SLD\ :sub:`solvent`) \* 
+where the mass *m*\ :sub:`i` is (SLD\ :sub:`i` - SLD\ :sub:`solvent`) \*
 (volume of the *N* pearls/rods). *V* is the total volume of the necklace.
 
-The 2D scattering intensity is the same as $P(q)$ above, regardless of the 
+The 2D scattering intensity is the same as $P(q)$ above, regardless of the
 orientation of the *q* vector.
 
@@ -54,5 +54,5 @@
 REFERENCE
 
-R Schweins and K Huber, *Particle Scattering Factor of Pearl Necklace Chains*, 
+R Schweins and K Huber, *Particle Scattering Factor of Pearl Necklace Chains*,
 *Macromol. Symp.* 211 (2004) 25-42 2004
 """
@@ -79,19 +79,19 @@
 
 #             ["name", "units", default, [lower, upper], "type","description"],
-parameters = [["radius", "Angstrom", 80.0, [0, inf], "volume", 
+parameters = [["radius", "Angstrom", 80.0, [0, inf], "volume",
                "Mean radius of the chained spheres"],
-              ["edge_separation", "Angstrom", 350.0, [0, inf], "volume", 
+              ["edge_separation", "Angstrom", 350.0, [0, inf], "volume",
                "Mean separation of chained particles"],
-              ["string_thickness", "Angstrom", 2.5, [0, inf], "volume", 
+              ["string_thickness", "Angstrom", 2.5, [0, inf], "volume",
                "Thickness of the chain linkage"],
-              ["number_of_pearls", "none", 3, [0, inf], "volume", 
+              ["number_of_pearls", "none", 3, [0, inf], "volume",
                "Mean number of pearls in each necklace"],
-              ["sld", "Angstrom^2", 1.0, [-inf, inf], "", 
+              ["sld", "Angstrom^2", 1.0, [-inf, inf], "",
                "Scattering length density of the chained spheres"],
-              ["string_sld", "Angstrom^2", 1.0, [-inf, inf], "", 
+              ["string_sld", "Angstrom^2", 1.0, [-inf, inf], "",
                "Scattering length density of the chain linkage"],
-              ["solvent_sld", "Angstrom^2", 6.3, [-inf, inf], "", 
+              ["solvent_sld", "Angstrom^2", 6.3, [-inf, inf], "",
                "Scattering length density of the solvent"],
-              ]
+             ]
 
 source = ["lib/Si.c", "pearl_necklace.c"]
@@ -110,5 +110,5 @@
 # names and the target sasview model name.
 oldname = 'PearlNecklaceModel'
-oldpars = dict(scale='scale',background='background',radius='radius',
+oldpars = dict(scale='scale', background='background', radius='radius',
                number_of_pearls='num_pearls', solvent_sld='sld_solv',
                string_thickness='thick_string', sld='sld_pearl',
Index: sasmodels/models/power_law.py
===================================================================
--- sasmodels/models/power_law.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
+++ sasmodels/models/power_law.py	(revision 841753cbac7824298effa1aec88a7ee977eb257f)
@@ -12,9 +12,9 @@
     I(q) = \text{scale} \cdot q^{-\text{power}} + \text{background}
 
-Note the minus sign in front of the exponent. The exponent *power* 
+Note the minus sign in front of the exponent. The exponent *power*
 should therefore be entered as a **positive** number for fitting.
 
-Also note that unlike many other models, *scale* in this model 
-is NOT explicitly related to a volume fraction. Be careful if 
+Also note that unlike many other models, *scale* in this model
+is NOT explicitly related to a volume fraction. Be careful if
 combining this model with other models.
 
@@ -31,12 +31,12 @@
 from numpy import inf, sqrt
 
-name =  "power_law"
+name = "power_law"
 title = "Simple power law with a flat background"
 
-description = """\
-        Evaluates the function
-        I(q) = scale * q^(-power) + background
-        NB: enter power as a positive number!
-        """
+description = """
+    Evaluates the function
+    I(q) = scale * q^(-power) + background
+    NB: enter power as a positive number!
+    """
 category = "shape-independent"
 
@@ -45,5 +45,6 @@
 
 # NB: Scale and Background are implicit parameters on every model
-def Iq(q,power):
+def Iq(q, power):
+    # pylint: disable=missing-docstring
     inten = (q**-power)
     return inten
@@ -51,4 +52,5 @@
 
 def Iqxy(qx, qy, *args):
+    # pylint: disable=missing-docstring
     return Iq(sqrt(qx ** 2 + qy ** 2), *args)
 Iqxy.vectorized = True # Iqxy accepts an array of qx, qy values
@@ -64,4 +66,5 @@
 
 tests = [
-        [ {'scale': 1.0, 'power': 4.0, 'background' : 0.0}, [0.0106939, 0.469418], [7.64644e+07, 20.5949]]
-        ]
+    [{'scale': 1.0, 'power': 4.0, 'background' : 0.0},
+     [0.0106939, 0.469418], [7.64644e+07, 20.5949]],
+    ]
Index: sasmodels/resolution.py
===================================================================
--- sasmodels/resolution.py	(revision fdc538a09f8fc03aca7f7563f14b6fc5ca0d178f)
+++ sasmodels/resolution.py	(revision 190fc2b167d3d89ca9c3b7d3128bdf8bde9f67d0)
@@ -5,4 +5,8 @@
 """
 from __future__ import division
+
+from scipy.special import erf
+from numpy import sqrt, log, log10
+import numpy as np
 
 __all__ = ["Resolution", "Perfect1D", "Pinhole1D", "Slit1D",
@@ -11,8 +15,4 @@
            "interpolate", "linear_extrapolation", "geometric_extrapolation",
            ]
-
-from scipy.special import erf
-from numpy import sqrt, log, log10
-import numpy as np
 
 MINIMUM_RESOLUTION = 1e-8
Index: sasmodels/resolution2d.py
===================================================================
--- sasmodels/resolution2d.py	(revision 9404dd3f599fcbe72790aae8a6d176e64f292ea6)
+++ sasmodels/resolution2d.py	(revision 841753cbac7824298effa1aec88a7ee977eb257f)
@@ -2,5 +2,5 @@
 #This software was developed by the University of Tennessee as part of the
 #Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
-#project funded by the US National Science Foundation. 
+#project funded by the US National Science Foundation.
 #See the license text in license.txt
 """
@@ -19,5 +19,5 @@
 ## Defaults
 NR = {'xhigh':10, 'high':5, 'med':5, 'low':3}
-NPHI ={'xhigh':20, 'high':12, 'med':6, 'low':4}
+NPHI = {'xhigh':20, 'high':12, 'med':6, 'low':4}
 
 class Pinhole2D(Resolution):
@@ -25,10 +25,10 @@
     Gaussian Q smearing class for SAS 2d data
     """
-     
+
     def __init__(self, data=None, index=None,
                  nsigma=NSIGMA, accuracy='Low', coords='polar'):
         """
         Assumption: equally spaced bins in dq_r, dq_phi space.
-        
+
         :param data: 2d data used to set the smearing parameters
         :param index: 1d array with len(data) to define the range
@@ -42,8 +42,8 @@
         ## number of bins in r axis for over-sampling
         self.nr = NR[accuracy.lower()]
-        ## number of bins in phi axis for over-sampling 
+        ## number of bins in phi axis for over-sampling
         self.nphi = NPHI[accuracy.lower()]
         ## maximum nsigmas
-        self.nsigma= nsigma
+        self.nsigma = nsigma
         self.coords = coords
         self._init_data(data, index)
@@ -85,7 +85,7 @@
     def _calc_res(self):
         """
-        Over sampling of r_nbins times phi_nbins, calculate Gaussian weights, 
+        Over sampling of r_nbins times phi_nbins, calculate Gaussian weights,
         then find smeared intensity
-        """    
+        """
         nr, nphi = self.nr, self.nphi
         # Total number of bins = # of bins
@@ -143,5 +143,5 @@
         if self.coords == 'polar':
             q_r = sqrt(qx**2 + qy**2)
-            qx_res = ( (dqx*cos(dphi) + q_r) * 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) +
@@ -167,70 +167,2 @@
         else:
             return theory
-
-"""
-if __name__ == '__main__':
-    ## Test w/ 2D linear function
-    x = 0.001*np.arange(1, 11)
-    dx = np.ones(len(x))*0.0003
-    y = 0.001*np.arange(1, 11)
-    dy = np.ones(len(x))*0.001
-    z = np.ones(10)
-    dz = sqrt(z)
-
-    from sas.dataloader import Data2D
-    #for i in range(10): print(i, 0.001 + i*0.008/9.0)
-    #for i in range(100): print(i, int(math.floor( (i/ (100/9.0)) )))
-    out = Data2D()
-    out.data = z
-    out.qx_data = x
-    out.qy_data = y
-    out.dqx_data = dx
-    out.dqy_data = dy
-    out.q_data = sqrt(dx * dx + dy * dy)
-    index = np.ones(len(x), dtype = bool)
-    out.mask = index
-    from sas.models.LineModel import LineModel
-    model = LineModel()
-    model.setParam("A", 0)
-
-    smear = Smearer2D(out, model, index)
-    #smear.set_accuracy('Xhigh')
-    value = smear.get_value()
-    ## All data are ones, so the smeared should also be ones.
-    print("Data length =", len(value))
-    print(" 2D linear function, I = 0 + 1*qy")
-    text = " Gaussian weighted averaging on a 2D linear function will "
-    text += "provides the results same as without the averaging."
-    print(text)
-    print("qx_data", "qy_data", "I_nonsmear", "I_smeared")
-    for ind in range(len(value)):
-        print(x[ind], y[ind], model.evalDistribution([x, y])[ind], value[ind])
-
-
-if __name__ == '__main__':
-    ## Another Test w/ constant function
-    x = 0.001*np.arange(1,11)
-    dx = np.ones(len(x))*0.001
-    y = 0.001*np.arange(1,11)
-    dy = np.ones(len(x))*0.001
-    z = np.ones(10)
-    dz = sqrt(z)
-
-    from DataLoader import Data2D
-    #for i in range(10): print(i, 0.001 + i*0.008/9.0)
-    #for i in range(100): print(i, int(math.floor( (i/ (100/9.0)) )))
-    out = Data2D()
-    out.data = z
-    out.qx_data = x
-    out.qy_data = y
-    out.dqx_data = dx
-    out.dqy_data = dy
-    index = np.ones(len(x), dtype = bool)
-    out.mask = index
-    from sas.models.Constant import Constant
-    model = Constant()
-
-    value = Smearer2D(out,model,index).get_value()
-    ## All data are ones, so the smeared values should also be ones.
-    print("Data length =",len(value), ", Data=",value)
-"""
Index: sasmodels/sesans.py
===================================================================
--- sasmodels/sesans.py	(revision 384d114757de20d9c246bd823485afb6835ec978)
+++ sasmodels/sesans.py	(revision 190fc2b167d3d89ca9c3b7d3128bdf8bde9f67d0)
@@ -43,11 +43,11 @@
     *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
-        G[i] = np.sum(integr)
+    G = np.zeros_like(SElength, 'd')
+    for i, SElength_i in enumerate(SElength):
+        integral = besselj(0, q*SElength_i)*Iq*q
+        G[i] = np.sum(integral)
 
     # [m^-1] step size in q, needed for integration
-    dq=(q[1]-q[0])*1e10
+    dq = (q[1]-q[0])*1e10
 
     # integration step, convert q into [m**-1] and 2 pi circle integration
