Index: extra/pylint.rc
===================================================================
--- extra/pylint.rc	(revision 53d0e2453502b7ea4d37f56a426f6b219910b99b)
+++ extra/pylint.rc	(revision 3c56da87f0995bdde77f4183b68e77b81584dde0)
@@ -20,8 +20,8 @@
 # List of plugins (as comma separated values of python modules names) to load,
 # usually to register additional checkers.
-load-plugins=
+load-plugins=pylint_numpy,pylint_pyopencl
 
 # Use multiple processes to speed up Pylint.
-jobs=4
+#jobs=4
 
 # Allow loading of arbitrary C extensions. Extensions are imported into the
@@ -57,5 +57,21 @@
 
 # Disable the message(s) with the given id(s).
-disable=W0702,W0613,W0703,W0142
+# http://pylint-messages.wikidot.com/all-codes
+disable=
+    multiple-statements,
+    global-statement,
+    bare-except,
+    broad-except,
+    bad-whitespace,
+    bad-continuation,
+    too-many-statements,
+    too-many-branches,
+    too-many-locals,
+    too-many-instance-attributes,
+    too-many-arguments,
+    star-args,
+    unbalanced-tuple-unpacking,
+    locally-disabled,
+    old-style-class,
 
 [REPORTS]
@@ -98,5 +114,12 @@
 
 # Good variable names which should always be accepted, separated by a comma
-good-names=i,j,k,ex,Run,_,x,y,z,qx,qy,qz,n,q,dx,dy,dz,id,Iq,dIq,Qx,Qy,Qz
+good-names=_,
+    input,
+    i,j,k,n,x,y,z,
+    q,qx,qy,qz,
+    dt,dx,dy,dz,id,
+    Iq,dIq,Qx,Qy,Qz,
+    p,
+    ER, call_ER, VR, call_VR,
 
 # Bad variable names which should always be refused, separated by a comma
@@ -111,5 +134,5 @@
 
 # Regular expression matching correct function names
-function-rgx=[a-z_][a-z0-9_]{2,30}$
+function-rgx=[a-z_][a-z0-9_]{2,30}([123]D)?$
 
 # Naming hint for function names
@@ -172,5 +195,6 @@
 # Regular expression which should only match function or class names that do
 # not require a docstring.
-no-docstring-rgx=__.*__
+#no-docstring-rgx=__.*__
+no-docstring-rgx=_.*
 
 # Minimum line length for functions/classes that require docstrings, shorter
@@ -181,8 +205,11 @@
 
 # Maximum number of characters on a single line.
-max-line-length=100
+#max-line-length=100
+max-line-length=80
 
 # Regexp for a line that is allowed to be longer than the limit.
-ignore-long-lines=^\s*(# )?<?https?://\S+>?$
+#ignore-long-lines=^\s*(# )?<?https?://\S+>?$
+# No comma in the last forty characters
+ignore-long-lines=([^-,+*/]{40},?|[^"]{40}")$
 
 # Allow the body of an if to be on the same line as the test if there is no
Index: extra/pylint_numpy.py
===================================================================
--- extra/pylint_numpy.py	(revision 3c56da87f0995bdde77f4183b68e77b81584dde0)
+++ extra/pylint_numpy.py	(revision 3c56da87f0995bdde77f4183b68e77b81584dde0)
@@ -0,0 +1,12 @@
+from astroid import MANAGER
+from astroid import nodes
+
+def register(linter):
+    MANAGER.register_transform(nodes.Module, transform)
+
+def transform(module):
+    #print("processing",module.name)
+    if module.name.startswith('numpy'):
+        if module.name == 'numpy': import numpy
+        elif module.name == 'numpy.random': import numpy.random
+
Index: extra/pylint_pyopencl.py
===================================================================
--- extra/pylint_pyopencl.py	(revision 3c56da87f0995bdde77f4183b68e77b81584dde0)
+++ extra/pylint_pyopencl.py	(revision 3c56da87f0995bdde77f4183b68e77b81584dde0)
@@ -0,0 +1,10 @@
+from astroid import MANAGER
+from astroid import nodes
+
+def register(linter):
+    MANAGER.register_transform(nodes.Module, transform)
+
+def transform(module):
+    #print("processing",module.name)
+    if module.name == 'pyopencl':
+        import pyopencl
Index: extra/run-pylint.py
===================================================================
--- extra/run-pylint.py	(revision 6c2eb83ac60a11ed41f03e08a500cf29bf2a7398)
+++ extra/run-pylint.py	(revision 3c56da87f0995bdde77f4183b68e77b81584dde0)
@@ -4,7 +4,11 @@
 
 def main():
+    envpath = os.environ.get('PYTHONPATH',None)
+    path = [envpath] if envpath else []
+    path.append(abspath(dirname(__file__))) # so we can find the plugins
+    os.environ['PYTHONPATH'] = ':'.join(path)
     root = abspath(joinpath(dirname(__file__), '..'))
     os.chdir(root)
-    cmd = "pylint --rcfile extra/pylint.rc -f parseable sasmodels > pylint_violations.txt"
+    cmd = "pylint --rcfile extra/pylint.rc -f parseable sasmodels"
     status = os.system(cmd)
     sys.exit(status)
Index: sasmodels/alignment.py
===================================================================
--- sasmodels/alignment.py	(revision ff7119b4ad066090cdffae80706883cef3d80399)
+++ sasmodels/alignment.py	(revision 3c56da87f0995bdde77f4183b68e77b81584dde0)
@@ -34,11 +34,11 @@
     return view
 
-def align_data(v, dtype, alignment=128):
+def align_data(x, dtype, alignment=128):
     """
     Return a copy of an array on the alignment boundary.
     """
-    # if v is contiguous, aligned, and of the correct type then just return v
-    view = align_empty(v.shape, dtype, alignment=alignment)
-    view[:] = v
+    # if x is contiguous, aligned, and of the correct type then just return x
+    view = align_empty(x.shape, dtype, alignment=alignment)
+    view[:] = x
     return view
 
Index: sasmodels/bumps_model.py
===================================================================
--- sasmodels/bumps_model.py	(revision b89f519c356bd2e11f1aff0bcd59de0fb7a62d1c)
+++ sasmodels/bumps_model.py	(revision 3c56da87f0995bdde77f4183b68e77b81584dde0)
@@ -8,7 +8,11 @@
 # CRUFT python 2.6
 if not hasattr(datetime.timedelta, 'total_seconds'):
-    def delay(dt): return dt.days * 86400 + dt.seconds + 1e-6 * dt.microseconds
+    def delay(dt):
+        """Return number date-time delta as number seconds"""
+        return dt.days * 86400 + dt.seconds + 1e-6 * dt.microseconds
 else:
-    def delay(dt): return dt.total_seconds()
+    def delay(dt):
+        """Return number date-time delta as number seconds"""
+        return dt.total_seconds()
 
 import numpy as np
@@ -141,18 +145,21 @@
     from sas.dataloader.manipulations import Boxcut
     if half == 'right':
-        data.mask += Boxcut(x_min=-np.inf, x_max=0.0, y_min=-np.inf, y_max=np.inf)(data)
+        data.mask += \
+            Boxcut(x_min=-np.inf, x_max=0.0, y_min=-np.inf, y_max=np.inf)(data)
     if half == 'left':
-        data.mask += Boxcut(x_min=0.0, x_max=np.inf, y_min=-np.inf, y_max=np.inf)(data)
-
-
-def set_top(data, max):
-    """
-    Chop the top off the data, above *max*.
+        data.mask += \
+            Boxcut(x_min=0.0, x_max=np.inf, y_min=-np.inf, y_max=np.inf)(data)
+
+
+def set_top(data, cutoff):
+    """
+    Chop the top off the data, above *cutoff*.
     """
     from sas.dataloader.manipulations import Boxcut
-    data.mask += Boxcut(x_min=-np.inf, x_max=np.inf, y_min=-np.inf, y_max=max)(data)
-
-
-def plot_data(data, iq, vmin=None, vmax=None, view='log'):
+    data.mask += \
+        Boxcut(x_min=-np.inf, x_max=np.inf, y_min=-np.inf, y_max=cutoff)(data)
+
+
+def plot_data(data, Iq, vmin=None, vmax=None, view='log'):
     """
     Plot the target value for the data.  This could be the data itself,
@@ -161,22 +168,22 @@
     *scale* can be 'log' for log scale data, or 'linear'.
     """
-    from numpy.ma import masked_array, masked
+    from numpy.ma import masked_array
     import matplotlib.pyplot as plt
     if hasattr(data, 'qx_data'):
-        iq = iq + 0
-        valid = np.isfinite(iq)
+        Iq = Iq + 0
+        valid = np.isfinite(Iq)
         if view == 'log':
-            valid[valid] = (iq[valid] > 0)
-            iq[valid] = np.log10(iq[valid])
+            valid[valid] = (Iq[valid] > 0)
+            Iq[valid] = np.log10(Iq[valid])
         elif view == 'q4':
-            iq[valid] = iq*(data.qx_data[valid]**2+data.qy_data[valid]**2)**2
-        iq[~valid | data.mask] = 0
-        #plottable = iq
-        plottable = masked_array(iq, ~valid | data.mask)
+            Iq[valid] = Iq*(data.qx_data[valid]**2+data.qy_data[valid]**2)**2
+        Iq[~valid | data.mask] = 0
+        #plottable = Iq
+        plottable = masked_array(Iq, ~valid | data.mask)
         xmin, xmax = min(data.qx_data), max(data.qx_data)
         ymin, ymax = min(data.qy_data), max(data.qy_data)
         try:
-            if vmin is None: vmin = iq[valid & ~data.mask].min()
-            if vmax is None: vmax = iq[valid & ~data.mask].max()
+            if vmin is None: vmin = Iq[valid & ~data.mask].min()
+            if vmax is None: vmax = Iq[valid & ~data.mask].max()
         except:
             vmin, vmax = 0, 1
@@ -186,13 +193,13 @@
     else: # 1D data
         if view == 'linear' or view == 'q4':
-            #idx = np.isfinite(iq)
+            #idx = np.isfinite(Iq)
             scale = data.x**4 if view == 'q4' else 1.0
-            plt.plot(data.x, scale*iq) #, '.')
+            plt.plot(data.x, scale*Iq) #, '.')
         else:
             # Find the values that are finite and positive
-            idx = np.isfinite(iq)
-            idx[idx] = iq[idx]>0
-            iq[~idx] = np.nan
-            plt.loglog(data.x, iq)
+            idx = np.isfinite(Iq)
+            idx[idx] = (Iq[idx] > 0)
+            Iq[~idx] = np.nan
+            plt.loglog(data.x, Iq)
 
 
@@ -220,4 +227,5 @@
     plt.plot(data.x, mresid, 'x')
 
+# pylint: disable=unused-argument
 def _plot_sesans(data, theory, view):
     import matplotlib.pyplot as plt
@@ -286,12 +294,12 @@
             q = sesans.make_q(data.sample.zacceptance, data.Rmax)
             self.index = slice(None, None)
-            self.iq = data.y
-            self.diq = data.dy
+            self.Iq = data.y
+            self.dIq = data.dy
             self._theory = np.zeros_like(q)
             q_vectors = [q]
         elif self.data_type == 'Iqxy':
             self.index = (data.mask == 0) & (~np.isnan(data.data))
-            self.iq = data.data[self.index]
-            self.diq = data.err_data[self.index]
+            self.Iq = data.data[self.index]
+            self.dIq = data.err_data[self.index]
             self._theory = np.zeros_like(data.data)
             if not partype['orientation'] and not partype['magnetic']:
@@ -301,6 +309,6 @@
         elif self.data_type == 'Iq':
             self.index = (data.x >= data.qmin) & (data.x <= data.qmax) & ~np.isnan(data.y)
-            self.iq = data.y[self.index]
-            self.diq = data.dy[self.index]
+            self.Iq = data.y[self.index]
+            self.dIq = data.dy[self.index]
             self._theory = np.zeros_like(data.y)
             q_vectors = [data.x]
@@ -316,5 +324,5 @@
         pars = []
         for p in model.info['parameters']:
-            name, default, limits, ptype = p[0], p[2], p[3], p[4]
+            name, default, limits = p[0], p[2], p[3]
             value = kw.pop(name, default)
             setattr(self, name, Parameter.default(value, name=name, limits=limits))
@@ -335,5 +343,6 @@
         self._parameter_names = pars
         if kw:
-            raise TypeError("unexpected parameters: %s" % (", ".join(sorted(kw.keys()))))
+            raise TypeError("unexpected parameters: %s"
+                            % (", ".join(sorted(kw.keys()))))
         self.update()
 
@@ -345,5 +354,5 @@
             Return the number of points
         """
-        return len(self.iq)
+        return len(self.Iq)
 
     def parameters(self):
@@ -365,8 +374,8 @@
             #self._theory[:] = self._fn.eval(pars, pd_pars)
             if self.data_type == 'sesans':
-                P = sesans.hankel(self.data.x, self.data.lam * 1e-9,
-                                  self.data.sample.thickness / 10, self._fn_inputs[0],
-                                  self._theory)
-                self._cache['theory'] = P
+                result = sesans.hankel(self.data.x, self.data.lam * 1e-9,
+                                       self.data.sample.thickness / 10,
+                                       self._fn_inputs[0], self._theory)
+                self._cache['theory'] = result
             else:
                 self._cache['theory'] = self._theory
@@ -375,13 +384,13 @@
     def residuals(self):
         #if np.any(self.err ==0): print "zeros in err"
-        return (self.theory()[self.index] - self.iq) / self.diq
+        return (self.theory()[self.index] - self.Iq) / self.dIq
 
     def nllf(self):
-        R = self.residuals()
+        delta = self.residuals()
         #if np.any(np.isnan(R)): print "NaN in residuals"
-        return 0.5 * np.sum(R ** 2)
-
-    def __call__(self):
-        return 2 * self.nllf() / self.dof
+        return 0.5 * np.sum(delta ** 2)
+
+    #def __call__(self):
+    #    return 2 * self.nllf() / self.dof
 
     def plot(self, view='log'):
@@ -402,8 +411,8 @@
         print "noise", noise
         if noise is None:
-            noise = self.diq[self.index]
+            noise = self.dIq[self.index]
         else:
             noise = 0.01 * noise
-            self.diq[self.index] = noise
+            self.dIq[self.index] = noise
         y = self.theory()
         y += y * np.random.randn(*y.shape) * noise
@@ -428,10 +437,11 @@
         relative = self.model.info['partype']['pd-rel']
         limits = self.model.info['limits']
-        disperser, value, npts, width, nsigma = \
-            [getattr(self, par + ext) for ext in ('_pd_type', '', '_pd_n', '_pd', '_pd_nsigma')]
-        v, w = weights.get_weights(
+        disperser, value, npts, width, nsigma = [
+            getattr(self, par + ext)
+            for ext in ('_pd_type', '', '_pd_n', '_pd', '_pd_nsigma')]
+        value, weight = weights.get_weights(
             disperser, int(npts.value), width.value, nsigma.value,
             value.value, limits[par], par in relative)
-        return v, w / w.max()
+        return value, weight / np.sum(weight)
 
     def __getstate__(self):
@@ -442,4 +452,5 @@
 
     def __setstate__(self, state):
+        # pylint: disable=attribute-defined-outside-init
         self.__dict__ = state
 
Index: sasmodels/convert.py
===================================================================
--- sasmodels/convert.py	(revision 2a74b99fb6cc3193463c2f7ca88f81fd745ec29a)
+++ sasmodels/convert.py	(revision 3c56da87f0995bdde77f4183b68e77b81584dde0)
@@ -37,4 +37,5 @@
     Convert model from old style parameter names to new style.
     """
+    _,_ = name,pars # lint
     raise NotImplementedError
     # need to load all new models in order to determine old=>new
Index: sasmodels/core.py
===================================================================
--- sasmodels/core.py	(revision 989005365096670e54afec254875979f35bd395a)
+++ sasmodels/core.py	(revision 3c56da87f0995bdde77f4183b68e77b81584dde0)
@@ -12,5 +12,6 @@
 try:
     from .kernelcl import load_model as load_model_cl
-except Exception,exc:
+except:
+    # pylint: disable=invalid-name
     load_model_cl = None
 from .kerneldll import load_model as load_model_dll
@@ -31,8 +32,8 @@
     Return a computation kernel from the model definition and the q input.
     """
-    input = model.make_input(q_vectors)
-    return model(input)
+    model_input = model.make_input(q_vectors)
+    return model(model_input)
 
-def get_weights(kernel, pars, name):
+def get_weights(info, pars, name):
     """
     Generate the distribution for parameter *name* given the parameter values
@@ -41,15 +42,15 @@
     Searches for "name", "name_pd", "name_pd_type", "name_pd_n", "name_pd_sigma"
     """
-    relative = name in kernel.info['partype']['pd-rel']
-    limits = kernel.info['limits']
+    relative = name in info['partype']['pd-rel']
+    limits = info['limits']
     disperser = pars.get(name+'_pd_type', 'gaussian')
-    value = pars.get(name, kernel.info['defaults'][name])
+    value = pars.get(name, info['defaults'][name])
     npts = pars.get(name+'_pd_n', 0)
     width = pars.get(name+'_pd', 0.0)
     nsigma = pars.get(name+'_pd_nsigma', 3.0)
-    v,w = weights.get_weights(
+    value,weight = weights.get_weights(
         disperser, npts, width, nsigma,
         value, limits[name], relative)
-    return v,w/np.sum(w)
+    return value,weight/np.sum(weight)
 
 def dispersion_mesh(pars):
@@ -61,38 +62,38 @@
     parameter set in the vector.
     """
-    values, weights = zip(*pars)
-    if len(values) > 1:
-        values = [v.flatten() for v in np.meshgrid(*values)]
-        weights = np.vstack([v.flatten() for v in np.meshgrid(*weights)])
-        weights = np.prod(weights, axis=0)
-    return values, weights
+    value, weight = zip(*pars)
+    if len(value) > 1:
+        value = [v.flatten() for v in np.meshgrid(*value)]
+        weight = np.vstack([v.flatten() for v in np.meshgrid(*weight)])
+        weight = np.prod(weight, axis=0)
+    return value, weight
 
 def call_kernel(kernel, pars, cutoff=1e-5):
     fixed_pars = [pars.get(name, kernel.info['defaults'][name])
                   for name in kernel.fixed_pars]
-    pd_pars = [get_weights(kernel, pars, name) for name in kernel.pd_pars]
+    pd_pars = [get_weights(kernel.info, pars, name) for name in kernel.pd_pars]
     return kernel(fixed_pars, pd_pars, cutoff=cutoff)
 
-def call_ER(kernel, pars):
-    ER = kernel.info.get('ER', None)
+def call_ER(info, pars):
+    ER = info.get('ER', None)
     if ER is None:
         return 1.0
     else:
-        vol_pars = [get_weights(kernel, pars, name)
-                    for name in kernel.info['partype']['volume']]
-        values, weights = dispersion_mesh(vol_pars)
-        fv = ER(*values)
+        vol_pars = [get_weights(info, pars, name)
+                    for name in info['partype']['volume']]
+        value, weight = dispersion_mesh(vol_pars)
+        individual_radii = ER(*value)
         #print values[0].shape, weights.shape, fv.shape
-        return np.sum(weights*fv) / np.sum(weights)
+        return np.sum(weight*individual_radii) / np.sum(weight)
 
-def call_VR(kernel, pars):
-    VR = kernel.info.get('VR', None)
+def call_VR(info, pars):
+    VR = info.get('VR', None)
     if VR is None:
         return 1.0
     else:
-        vol_pars = [get_weights(kernel, pars, name)
-                    for name in kernel.info['partype']['volume']]
-        values, weights = dispersion_mesh(vol_pars)
-        whole,part = VR(*values)
-        return np.sum(weights*part)/np.sum(weights*whole)
+        vol_pars = [get_weights(info, pars, name)
+                    for name in info['partype']['volume']]
+        value, weight = dispersion_mesh(vol_pars)
+        whole,part = VR(*value)
+        return np.sum(weight*part)/np.sum(weight*whole)
 
Index: sasmodels/generate.py
===================================================================
--- sasmodels/generate.py	(revision 33e91b15304920ad5b2b99b7f595293978770dbe)
+++ sasmodels/generate.py	(revision 3c56da87f0995bdde77f4183b68e77b81584dde0)
@@ -566,10 +566,10 @@
 
 def demo_time():
+    from .models import cylinder
     import datetime
-    from .models import cylinder
     tic = datetime.datetime.now()
-    toc = lambda: (datetime.datetime.now() - tic).total_seconds()
     make(cylinder)
-    print "time:", toc()
+    toc = (datetime.datetime.now() - tic).total_seconds()
+    print "time:", toc
 
 def main():
Index: sasmodels/kernelcl.py
===================================================================
--- sasmodels/kernelcl.py	(revision c85db69ebc61ae093eb2e72b7a6bf64a0142d56d)
+++ sasmodels/kernelcl.py	(revision 3c56da87f0995bdde77f4183b68e77b81584dde0)
@@ -30,6 +30,6 @@
 try:
     import pyopencl as cl
-    context = cl.create_some_context(interactive=False)
-    del context
+    # Ask OpenCL for the default context so that we know that one exists
+    cl.create_some_context(interactive=False)
 except Exception, exc:
     warnings.warn(str(exc))
@@ -160,5 +160,5 @@
 
         if not self.context:
-            self.context = self._find_context()
+            self.context = _get_default_context()
 
         # Byte boundary for data alignment
@@ -178,18 +178,4 @@
             warnings.warn("the environment variable 'PYOPENCL_CTX' might not be set correctly")
 
-    def _find_context(self):
-        default = None
-        for platform in cl.get_platforms():
-            for device in platform.get_devices():
-                if device.type == cl.device_type.GPU:
-                    return cl.Context([device])
-                if default is None:
-                    default = device
-
-        if not default:
-            raise RuntimeError("OpenCL device not found")
-
-        return cl.Context([default])
-
     def compile_program(self, name, source, dtype):
         if name not in self.compiled:
@@ -203,4 +189,18 @@
             del self.compiled[name]
 
+def _get_default_context():
+    default = None
+    for platform in cl.get_platforms():
+        for device in platform.get_devices():
+            if device.type == cl.device_type.GPU:
+                return cl.Context([device])
+            if default is None:
+                default = device
+
+    if not default:
+        raise RuntimeError("OpenCL device not found")
+
+    return cl.Context([default])
+
 
 class GpuModel(object):
@@ -234,5 +234,6 @@
             raise TypeError("data and kernel have different types")
         if self.program is None:
-            self.program = environment().compile_program(self.info['name'], self.source, self.dtype)
+            compiler = environment().compile_program
+            self.program = compiler(self.info['name'], self.source, self.dtype)
         kernel_name = generate.kernel_name(self.info, input_value.is_2D)
         kernel = getattr(self.program, kernel_name)
@@ -303,5 +304,5 @@
     *info* is the module information
 
-    *input* is the DllInput q vectors at which the kernel should be
+    *q_input* is the DllInput q vectors at which the kernel should be
     evaluated.
 
@@ -314,10 +315,10 @@
     Call :meth:`release` when done with the kernel instance.
     """
-    def __init__(self, kernel, info, input):
-        self.input = input
+    def __init__(self, kernel, info, q_input):
+        self.q_input = q_input
         self.kernel = kernel
         self.info = info
-        self.res = np.empty(input.nq, input.dtype)
-        dim = '2d' if input.is_2D else '1d'
+        self.res = np.empty(q_input.nq, q_input.dtype)
+        dim = '2d' if q_input.is_2D else '1d'
         self.fixed_pars = info['partype']['fixed-' + dim]
         self.pd_pars = info['partype']['pd-' + dim]
@@ -327,23 +328,24 @@
         env = environment()
         self.loops_b = [cl.Buffer(env.context, mf.READ_WRITE,
-                                  2 * MAX_LOOPS * input.dtype.itemsize)
+                                  2 * MAX_LOOPS * q_input.dtype.itemsize)
                         for _ in env.queues]
         self.res_b = [cl.Buffer(env.context, mf.READ_WRITE,
-                                input.global_size[0] * input.dtype.itemsize)
+                                q_input.global_size[0] * q_input.dtype.itemsize)
                       for _ in env.queues]
 
 
     def __call__(self, fixed_pars, pd_pars, cutoff=1e-5):
-        real = np.float32 if self.input.dtype == generate.F32 else np.float64
+        real = np.float32 if self.q_input.dtype == generate.F32 else np.float64
 
         device_num = 0
         queuei = environment().queues[device_num]
         res_bi = self.res_b[device_num]
-        nq = np.uint32(self.input.nq)
+        nq = np.uint32(self.q_input.nq)
         if pd_pars:
             cutoff = real(cutoff)
             loops_N = [np.uint32(len(p[0])) for p in pd_pars]
-            loops = np.hstack(pd_pars) if pd_pars else np.empty(0, dtype=self.input.dtype)
-            loops = np.ascontiguousarray(loops.T, self.input.dtype).flatten()
+            loops = np.hstack(pd_pars) \
+                if pd_pars else np.empty(0, dtype=self.q_input.dtype)
+            loops = np.ascontiguousarray(loops.T, self.q_input.dtype).flatten()
             #print "loops",Nloops, loops
 
@@ -357,11 +359,11 @@
             loops_l = cl.LocalMemory(len(loops.data))
             #ctx = environment().context
-            #loops_bi = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=loops)
+            #loops_bi = cl.Buffer(ctx, mf.READ_ONLY|mf.COPY_HOST_PTR, hostbuf=loops)
             dispersed = [loops_bi, loops_l, cutoff] + loops_N
         else:
             dispersed = []
         fixed = [real(p) for p in fixed_pars]
-        args = self.input.q_buffers + [res_bi, nq] + dispersed + fixed
-        self.kernel(queuei, self.input.global_size, None, *args)
+        args = self.q_input.q_buffers + [res_bi, nq] + dispersed + fixed
+        self.kernel(queuei, self.q_input.global_size, None, *args)
         cl.enqueue_copy(queuei, self.res, res_bi)
 
Index: sasmodels/kerneldll.py
===================================================================
--- sasmodels/kerneldll.py	(revision f734e7d992ea6f2199f478553447ff0dbe047559)
+++ sasmodels/kerneldll.py	(revision 3c56da87f0995bdde77f4183b68e77b81584dde0)
@@ -19,5 +19,5 @@
     COMPILE = "gcc -shared -fPIC -std=c99 -O2 -Wall %(source)s -o %(output)s -lm"
 elif os.name == 'nt':
-    # make sure vcvarsall.bat is called first in order to set compiler, headers, lib paths, etc.
+    # call vcvarsall.bat before compiling to set path, headers, libs, etc.
     if "VCINSTALLDIR" in os.environ:
         # MSVC compiler is available, so use it.
@@ -25,5 +25,6 @@
         # TODO: maybe don't use randomized name for the c file
         COMPILE = "cl /nologo /Ox /MD /W3 /GS- /DNDEBUG /Tp%(source)s /openmp /link /DLL /INCREMENTAL:NO /MANIFEST /OUT:%(output)s"
-        # Can't find VCOMP90.DLL (don't know why), so remove openmp support from windows compiler build
+        # Can't find VCOMP90.DLL (don't know why), so remove openmp support
+        # from windows compiler build
         #COMPILE = "cl /nologo /Ox /MD /W3 /GS- /DNDEBUG /Tp%(source)s /link /DLL /INCREMENTAL:NO /MANIFEST /OUT:%(output)s"
     else:
@@ -123,9 +124,10 @@
         self.__dict__ = state
 
-    def __call__(self, input):
+    def __call__(self, q_input):
         if self.dll is None: self._load_dll()
-        kernel = self.Iqxy if input.is_2D else self.Iq
-        return DllKernel(kernel, self.info, input)
+        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):
         """
@@ -150,5 +152,5 @@
     *info* is the module information
 
-    *input* is the DllInput q vectors at which the kernel should be
+    *q_input* is the DllInput q vectors at which the kernel should be
     evaluated.
 
@@ -161,10 +163,10 @@
     Call :meth:`release` when done with the kernel instance.
     """
-    def __init__(self, kernel, info, input):
+    def __init__(self, kernel, info, q_input):
         self.info = info
-        self.input = input
+        self.q_input = q_input
         self.kernel = kernel
-        self.res = np.empty(input.nq, input.dtype)
-        dim = '2d' if input.is_2D else '1d'
+        self.res = np.empty(q_input.nq, q_input.dtype)
+        dim = '2d' if q_input.is_2D else '1d'
         self.fixed_pars = info['partype']['fixed-'+dim]
         self.pd_pars = info['partype']['pd-'+dim]
@@ -174,12 +176,12 @@
 
     def __call__(self, fixed_pars, pd_pars, cutoff):
-        real = np.float32 if self.input.dtype == F32 else np.float64
+        real = np.float32 if self.q_input.dtype == F32 else np.float64
 
-        nq = c_int(self.input.nq)
+        nq = c_int(self.q_input.nq)
         if pd_pars:
             cutoff = real(cutoff)
             loops_N = [np.uint32(len(p[0])) for p in pd_pars]
             loops = np.hstack(pd_pars)
-            loops = np.ascontiguousarray(loops.T, self.input.dtype).flatten()
+            loops = np.ascontiguousarray(loops.T, self.q_input.dtype).flatten()
             p_loops = loops.ctypes.data
             dispersed = [p_loops, cutoff] + loops_N
@@ -187,5 +189,5 @@
             dispersed = []
         fixed = [real(p) for p in fixed_pars]
-        args = self.input.q_pointers + [self.p_res, nq] + dispersed + fixed
+        args = self.q_input.q_pointers + [self.p_res, nq] + dispersed + fixed
         #print pars
         self.kernel(*args)
Index: sasmodels/kernelpy.py
===================================================================
--- sasmodels/kernelpy.py	(revision c85db69ebc61ae093eb2e72b7a6bf64a0142d56d)
+++ sasmodels/kernelpy.py	(revision 3c56da87f0995bdde77f4183b68e77b81584dde0)
@@ -7,9 +7,13 @@
     def __init__(self, info):
         self.info = info
-    def __call__(self, input_value):
-        kernel = self.info['Iqxy'] if input_value.is_2D else self.info['Iq']
-        return PyKernel(kernel, self.info, input_value)
+
+    def __call__(self, q_input):
+        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):
         pass
@@ -52,5 +56,5 @@
     *info* is the module information
 
-    *input* is the DllInput q vectors at which the kernel should be
+    *q_input* is the DllInput q vectors at which the kernel should be
     evaluated.
 
@@ -63,9 +67,9 @@
     Call :meth:`release` when done with the kernel instance.
     """
-    def __init__(self, kernel, info, input):
+    def __init__(self, kernel, info, q_input):
         self.info = info
-        self.input = input
-        self.res = np.empty(input.nq, input.dtype)
-        dim = '2d' if input.is_2D else '1d'
+        self.q_input = q_input
+        self.res = np.empty(q_input.nq, q_input.dtype)
+        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
@@ -73,8 +77,10 @@
             if dim == '2d':
                 def vector_kernel(qx, qy, *args):
-                    return np.array([kernel(qxi, qyi, *args) for qxi, qyi in zip(qx, qy)])
+                    return np.array([kernel(qxi, qyi, *args)
+                                     for qxi, qyi in zip(qx, qy)])
             else:
                 def vector_kernel(q, *args):
-                    return np.array([kernel(qi, *args) for qi in q])
+                    return np.array([kernel(qi, *args)
+                                     for qi in q])
             self.kernel = vector_kernel
         else:
@@ -86,9 +92,12 @@
         # First two fixed pars are scale and background
         pars = [p[0] for p in info['parameters'][2:]]
-        offset = len(self.input.q_vectors)
-        self.args = self.input.q_vectors + [None] * len(pars)
-        self.fixed_index = np.array([pars.index(p) + offset for p in fixed_pars[2:]])
-        self.pd_index = np.array([pars.index(p) + offset for p in pd_pars])
-        self.vol_index = np.array([pars.index(p) + offset for p in vol_pars])
+        offset = len(self.q_input.q_vectors)
+        self.args = self.q_input.q_vectors + [None] * len(pars)
+        self.fixed_index = np.array([pars.index(p) + offset
+                                     for p in fixed_pars[2:]])
+        self.pd_index = np.array([pars.index(p) + offset
+                                  for p in pd_pars])
+        self.vol_index = np.array([pars.index(p) + offset
+                                   for p in vol_pars])
         try: self.theta_index = pars.index('theta') + offset
         except ValueError: self.theta_index = -1
@@ -113,5 +122,5 @@
 
     def release(self):
-        self.input = None
+        self.q_input = None
 
 def _loops(form, form_volume, cutoff, scale, background,
@@ -194,7 +203,7 @@
             args[pd_index[0]] = fast_value[fast_index]
             weight[0] = fast_weight[fast_index]
-        # This computes the weight, and if it is sufficient, calls the scattering
-        # function and adds it to the total.  If there is a volume normalization,
-        # it will also be added here.
+        # This computes the weight, and if it is sufficient, calls the
+        # scattering function and adds it to the total.  If there is a volume
+        # normalization, it will also be added here.
         # Note: make sure this is consistent with the code in PY_LOOP_BODY!!
         # Note: can precompute w1*w2*...*wn
@@ -204,8 +213,9 @@
             positive = (I >= 0.0)
 
-            # Note: can precompute spherical correction if theta_index is not the fast index
-            # Correction factor for spherical integration p(theta) I(q) sin(theta) dtheta
-            #spherical_correction = abs(sin(pi*args[theta_index])) if theta_index>=0 else 1.0
-            spherical_correction = abs(cos(pi * args[theta_index])) * pi / 2 if theta_index >= 0 else 1.0
+            # Note: can precompute spherical correction if theta_index is not
+            # the fast index. Correction factor for spherical integration
+            #spherical_correction = abs(cos(pi*args[phi_index])) if phi_index>=0 else 1.0
+            spherical_correction = (abs(cos(pi * args[theta_index])) * pi / 2
+                                    if theta_index >= 0 else 1.0)
             #spherical_correction = 1.0
             ret += w * I * spherical_correction * positive
@@ -213,8 +223,8 @@
 
             # Volume normalization.
-            # If there are "volume" polydispersity parameters, then these will be used
-            # to call the form_volume function from the user supplied kernel, and accumulate
-            # a normalized weight.
-            # Note: can precompute volume norm if the fast index is not a volume index
+            # If there are "volume" polydispersity parameters, then these
+            # will be used to call the form_volume function from the user
+            # supplied kernel, and accumulate a normalized weight.
+            # Note: can precompute volume norm if fast index is not a volume
             if form_volume:
                 vol_args = [args[index] for index in vol_index]
Index: sasmodels/model_test.py
===================================================================
--- sasmodels/model_test.py	(revision d6adfbe1e1180ee1d68c9d8ac7b8a93e0e155cca)
+++ sasmodels/model_test.py	(revision 3c56da87f0995bdde77f4183b68e77b81584dde0)
@@ -5,7 +5,7 @@
 Usage::
 
-     python -m sasmodels.model_test [opencl|dll|opencl_and_dll] model1 model2 ...
-
-     if model1 is 'all', then all except the remaining models will be tested
+    python -m sasmodels.model_test [opencl|dll|opencl_and_dll] model1 model2 ...
+
+    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),
@@ -60,7 +60,7 @@
     Example::
         >>> D = {}
-        >>> try: 
+        >>> try:
         ...    print D['hello']
-        ... except Exception,exc: 
+        ... except Exception,exc:
         ...    annotate_exception(exc, "while accessing 'D'")
         ...    raise
@@ -78,6 +78,6 @@
         except:
             exc.args = (" ".join((str(exc),msg)),)
-    
-def suite(loaders, models):
+
+def make_suite(loaders, models):
 
     ModelTestCase = _hide_model_case_from_nosetests()
@@ -100,5 +100,5 @@
             ]
         tests = smoke_tests + getattr(model_definition, 'tests', [])
-        
+
         if tests: # in case there are no smoke tests...
             #print '------'
@@ -143,38 +143,5 @@
                 model = self.loader(self.definition)
                 for test in self.tests:
-                    pars, Q, I = test
-
-                    if not isinstance(I, list):
-                        I = [I]
-                    if not isinstance(Q, list):
-                        Q = [Q]
-
-                    self.assertEqual(len(I), len(Q))
-
-                    if Q[0] == 'ER':
-                        Iq = [call_ER(kernel, pars)]
-                    elif Q[0] == 'VR':
-                        Iq = [call_VR(kernel, pars)]
-                    elif isinstance(Q[0], tuple):
-                        Qx,Qy = zip(*Q)
-                        Q_vectors = [np.array(Qx), np.array(Qy)]
-                        kernel = make_kernel(model, Q_vectors)
-                        Iq = call_kernel(kernel, pars)
-                    else:
-                        Q_vectors = [np.array(Q)]
-                        kernel = make_kernel(model, Q_vectors)
-                        Iq = call_kernel(kernel, pars)
-
-                    self.assertGreater(len(Iq), 0)
-                    self.assertEqual(len(I), len(Iq))
-
-                    for q, i, iq in zip(Q, I, Iq):
-                        if i is None:
-                            # smoke test --- make sure it runs and produces a value
-                            self.assertTrue(np.isfinite(iq), 'q:%s; not finite; actual:%s' % (q, iq))
-                        else:
-                            err = abs(i - iq)
-                            nrm = abs(i)
-                            self.assertLess(err * 10**5, nrm, 'q:%s; expected:%s; actual:%s' % (q, i, iq))
+                    self._run_one_test(model, test)
 
             except Exception,exc:
@@ -182,4 +149,42 @@
                 raise
 
+        def _run_one_test(self, model, test):
+            pars, x, y = test
+
+            if not isinstance(y, list):
+                y = [y]
+            if not isinstance(x, list):
+                x = [x]
+
+            self.assertEqual(len(y), len(x))
+
+            if x[0] == 'ER':
+                actual = [call_ER(model.info, pars)]
+            elif x[0] == 'VR':
+                actual = [call_VR(model.info, pars)]
+            elif isinstance(x[0], tuple):
+                Qx,Qy = zip(*x)
+                q_vectors = [np.array(Qx), np.array(Qy)]
+                kernel = make_kernel(model, q_vectors)
+                actual = call_kernel(kernel, pars)
+            else:
+                q_vectors = [np.array(x)]
+                kernel = make_kernel(model, q_vectors)
+                actual = call_kernel(kernel, pars)
+
+            self.assertGreater(len(actual), 0)
+            self.assertEqual(len(y), len(actual))
+
+            for xi, yi, actual_yi in zip(x, y, actual):
+                if yi is None:
+                    # smoke test --- make sure it runs and produces a value
+                    self.assertTrue(np.isfinite(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))
+
     return ModelTestCase
 
@@ -187,5 +192,5 @@
 # let nosetests sniff out the tests
 def model_tests():
-    tests = suite(['opencl','dll'],['all'])
+    tests = make_suite(['opencl','dll'],['all'])
     for test_i in tests:
         yield test_i.runTest
@@ -218,5 +223,5 @@
     #run_tests(loaders, models)
     runner = unittest.TextTestRunner()
-    result = runner.run(suite(loaders, models))
+    result = runner.run(make_suite(loaders, models))
     return 1 if result.failures or result.errors else 0
 
Index: sasmodels/models/barbell.py
===================================================================
--- sasmodels/models/barbell.py	(revision 989005365096670e54afec254875979f35bd395a)
+++ sasmodels/models/barbell.py	(revision 3c56da87f0995bdde77f4183b68e77b81584dde0)
@@ -4,8 +4,11 @@
 r"""
 
-Calculates the scattering from a barbell-shaped cylinder (This model simply becomes the DumBellModel when the length of
-the cylinder, *L*, is set to zero). That is, a sphereocylinder with spherical end caps that have a radius larger than
-that of the cylinder and the center of the end cap radius lies outside of the cylinder. All dimensions of the BarBell
-are considered to be monodisperse. See the diagram for the details of the geometry and restrictions on parameter values.
+Calculates the scattering from a barbell-shaped cylinder (This model simply
+becomes the DumBellModel when the length of the cylinder, *L*, is set to zero).
+That is, a sphereocylinder with spherical end caps that have a radius larger
+than that of the cylinder and the center of the end cap radius lies outside
+of the cylinder. All dimensions of the BarBell are considered to be
+monodisperse. See the diagram for the details of the geometry and restrictions
+on parameter values.
 
 Definition
@@ -18,8 +21,10 @@
 .. image:: img/barbell_geometry.jpg
 
-where *r* is the radius of the cylinder. All other parameters are as defined in the diagram.
+where *r* is the radius of the cylinder. All other parameters are as defined
+in the diagram.
 
 Since the end cap radius
-*R* >= *r* and by definition for this geometry *h* < 0, *h* is then defined by *r* and *R* as
+*R* >= *r* and by definition for this geometry *h* < 0, *h* is then
+defined by *r* and *R* as
 
 *h* = -1 \* sqrt(*R*\ :sup:`2` - *r*\ :sup:`2`)
@@ -46,7 +51,9 @@
              \over QR\sin\theta \left(1-t^2\right)^{1/2}}
 
-The < > brackets denote an average of the structure over all orientations. <*A* :sup:`2`\ *(q)*> is then the form
-factor, *P(q)*. The scale factor is equivalent to the volume fraction of cylinders, each of volume, *V*. Contrast is
-the difference of scattering length densities of the cylinder and the surrounding solvent.
+The < > brackets denote an average of the structure over all orientations.
+<*A* :sup:`2`\ *(q)*> is then the form factor, *P(q)*. The scale factor is
+equivalent to the volume fraction of cylinders, each of volume, *V*. Contrast
+is the difference of scattering length densities of the cylinder and the
+surrounding solvent.
 
 The volume of the barbell is
@@ -69,5 +76,6 @@
         \left( 4R^3 6R^2h - 2h^3 + 3r^2L \right)^{-1}
 
-**The requirement that** *R* >= *r* **is not enforced in the model!** It is up to you to restrict this during analysis.
+**The requirement that** *R* >= *r* **is not enforced in the model!** It is
+up to you to restrict this during analysis.
 
 This example dataset is produced by running the Macro PlotBarbell(),
@@ -79,6 +87,7 @@
 *Figure. 1D plot using the default values (w/256 data point).*
 
-For 2D data: The 2D scattering intensity is calculated similar to the 2D cylinder model. For example, for
-|theta| = 45 deg and |phi| = 0 deg with default values for other parameters
+For 2D data: The 2D scattering intensity is calculated similar to the 2D
+cylinder model. For example, for |theta| = 45 deg and |phi| = 0 deg with
+default values for other parameters
 
 .. image:: img/barbell_2d.jpg
@@ -102,5 +111,5 @@
 
 """
-from numpy import pi, inf
+from numpy import inf
 
 name = "barbell"
Index: sasmodels/models/bcc.py
===================================================================
--- sasmodels/models/bcc.py	(revision 989005365096670e54afec254875979f35bd395a)
+++ sasmodels/models/bcc.py	(revision 3c56da87f0995bdde77f4183b68e77b81584dde0)
@@ -80,5 +80,6 @@
 .. image:: img/bcc_1d.jpg
 
-*Figure. 1D plot in the linear scale using the default values (w/200 data point).*
+*Figure. 1D plot in the linear scale using the default values
+(w/200 data point).*
 
 The 2D (Anisotropic model) is based on the reference below where $I(q)$ is
@@ -103,5 +104,5 @@
 """
 
-from numpy import pi, inf
+from numpy import inf
 
 name = "bcc_paracrystal"
@@ -120,5 +121,5 @@
     [ "radius", "Ang",  40, [0, inf], "volume","Particle radius" ],
     [ "sld", "1e-6/Ang^2", 4, [-inf,inf], "", "Particle scattering length density" ],
-    [ "solvent_sld", "1e-6/Ang^2", 1, [-inf,inf], "","Solvent scattering length density" ],
+    [ "solvent_sld", "1e-6/Ang^2", 1, [-inf,inf], "", "Solvent scattering length density" ],
     [ "theta", "degrees", 60, [-inf, inf], "orientation","In plane angle" ],
     [ "phi", "degrees", 60, [-inf, inf], "orientation","Out of plane angle" ],
Index: sasmodels/models/broad_peak.py
===================================================================
--- sasmodels/models/broad_peak.py	(revision a5d0d008778409dd90a86aa3db6037187ed3ff57)
+++ sasmodels/models/broad_peak.py	(revision 3c56da87f0995bdde77f4183b68e77b81584dde0)
@@ -6,6 +6,6 @@
 layered structures, etc.
 
-The d-spacing corresponding to the broad peak is a characteristic distance 
-between the scattering inhomogeneities (such as in lamellar, cylindrical, or 
+The d-spacing corresponding to the broad peak is a characteristic distance
+between the scattering inhomogeneities (such as in lamellar, cylindrical, or
 spherical morphologies, or for bicontinuous structures).
 
@@ -44,6 +44,5 @@
 """
 
-import numpy as np
-from numpy import pi, inf, sin, cos, sqrt, exp, log, fabs
+from numpy import inf, sqrt
 
 name = "broad_peak"
@@ -82,20 +81,13 @@
 
 
-#def form_volume():
-#    return 1
-
 def Iq(q, porod_scale, porod_exp, lorentz_scale, lorentz_length, peak_pos, lorentz_exp):
-    inten = porod_scale/pow(q,porod_exp) + lorentz_scale/(1.0 \
-        + pow((fabs(q-peak_pos)*lorentz_length),lorentz_exp))
-    return inten  
-
-# FOR VECTORIZED VERSION, UNCOMMENT THE NEXT LINE
-Iq.vectorized = True
+    inten = (porod_scale/q**porod_exp + lorentz_scale
+        / (1.0 + (abs(q-peak_pos)*lorentz_length)**lorentz_exp))
+    return inten
+Iq.vectorized = True  # Iq accepts an array of Q values
 
 def Iqxy(qx, qy, *args):
     return Iq(sqrt(qx**2 + qy**2), *args)
-
-# FOR VECTORIZED VERSION, UNCOMMENT THE NEXT LINE
-Iqxy.vectorized = True
+Iqxy.vectorized = True # Iqxy accepts an array of Qx, Qy values
 
 
@@ -105,6 +97,7 @@
     lorentz_scale=10,lorentz_length=50, peak_pos=0.1, lorentz_exp=2,
     )
+
 oldname = "BroadPeakModel"
-oldpars = dict(porod_scale='scale_p', porod_exp='exponent_p', 
-        lorentz_scale='scale_l', lorentz_length='length_l', peak_pos='q_peak', 
+oldpars = dict(porod_scale='scale_p', porod_exp='exponent_p',
+        lorentz_scale='scale_l', lorentz_length='length_l', peak_pos='q_peak',
         lorentz_exp='exponent_l', scale=None)
Index: sasmodels/models/ellipsoid.py
===================================================================
--- sasmodels/models/ellipsoid.py	(revision a5d0d008778409dd90a86aa3db6037187ed3ff57)
+++ sasmodels/models/ellipsoid.py	(revision 3c56da87f0995bdde77f4183b68e77b81584dde0)
@@ -116,5 +116,5 @@
 """
 
-from numpy import pi, inf
+from numpy import inf
 
 name = "ellipsoid"
Index: sasmodels/models/fcc.py
===================================================================
--- sasmodels/models/fcc.py	(revision 989005365096670e54afec254875979f35bd395a)
+++ sasmodels/models/fcc.py	(revision 3c56da87f0995bdde77f4183b68e77b81584dde0)
@@ -3,7 +3,9 @@
 #note - calculation requires double precision
 r"""
-Calculates the scattering from a **face-centered cubic lattice** with paracrystalline distortion. Thermal vibrations
-are considered to be negligible, and the size of the paracrystal is infinitely large. Paracrystalline distortion is
-assumed to be isotropic and characterized by a Gaussian distribution.
+Calculates the scattering from a **face-centered cubic lattice** with
+paracrystalline distortion. Thermal vibrations are considered to be
+negligible, and the size of the paracrystal is infinitely large.
+Paracrystalline distortion is assumed to be isotropic and characterized by
+a Gaussian distribution.
 
 The returned value is scaled to units of |cm^-1|\ |sr^-1|, absolute scale.
@@ -16,17 +18,20 @@
 .. image:: img/image158.jpg
 
-where *scale* is the volume fraction of spheres, *Vp* is the volume of the primary particle, *V(lattice)* is a volume
-correction for the crystal structure, *P(q)* is the form factor of the sphere (normalized), and *Z(q)* is the
-paracrystalline structure factor for a face-centered cubic structure.
+where *scale* is the volume fraction of spheres, *Vp* is the volume of
+the primary particle, *V(lattice)* is a volume correction for the crystal
+structure, *P(q)* is the form factor of the sphere (normalized), and *Z(q)*
+is the paracrystalline structure factor for a face-centered cubic structure.
 
-Equation (1) of the 1990 reference is used to calculate *Z(q)*, using equations (23)-(25) from the 1987 paper for
-*Z1*\ , *Z2*\ , and *Z3*\ .
+Equation (1) of the 1990 reference is used to calculate *Z(q)*, using
+equations (23)-(25) from the 1987 paper for *Z1*\ , *Z2*\ , and *Z3*\ .
 
-The lattice correction (the occupied volume of the lattice) for a face-centered cubic structure of particles of radius
-*R* and nearest neighbor separation *D* is
+The lattice correction (the occupied volume of the lattice) for a
+face-centered cubic structure of particles of radius *R* and nearest
+neighbor separation *D* is
 
 .. image:: img/image159.jpg
 
-The distortion factor (one standard deviation) of the paracrystal is included in the calculation of *Z(q)*
+The distortion factor (one standard deviation) of the paracrystal is
+included in the calculation of *Z(q)*
 
 .. image:: img/image160.jpg
@@ -42,16 +47,19 @@
 .. image:: img/image162.jpg
 
-where for a face-centered cubic lattice *h*\ , *k*\ , *l* all odd or all even are allowed and reflections where
-*h*\ , *k*\ , *l* are mixed odd/even are forbidden. Thus the peak positions correspond to (just the first 5)
+where for a face-centered cubic lattice *h*\ , *k*\ , *l* all odd or all
+even are allowed and reflections where *h*\ , *k*\ , *l* are mixed odd/even
+are forbidden. Thus the peak positions correspond to (just the first 5)
 
 .. image:: img/image163.jpg
 
-**NB: The calculation of** *Z(q)* **is a double numerical integral that must be carried out with a high density of**
-**points to properly capture the sharp peaks of the paracrystalline scattering.** So be warned that the calculation is
-SLOW. Go get some coffee. Fitting of any experimental data must be resolution smeared for any meaningful fit. This
-makes a triple integral. Very, very slow. Go get lunch!
+**NB: The calculation of** *Z(q)* **is a double numerical integral that
+must be carried out with a high density of** **points to properly capture
+the sharp peaks of the paracrystalline scattering.** So be warned that the
+calculation is SLOW. Go get some coffee. Fitting of any experimental data
+must be resolution smeared for any meaningful fit. This makes a triple
+integral. Very, very slow. Go get lunch!
 
-This example dataset is produced using 200 data points, *qmin* = 0.01 |Ang^-1|, *qmax* = 0.1 |Ang^-1| and the above
-default values.
+This example dataset is produced using 200 data points, *qmin* = 0.01 |Ang^-1|,
+*qmax* = 0.1 |Ang^-1| and the above default values.
 
 .. image:: img/image164.jpg
@@ -59,7 +67,8 @@
 *Figure. 1D plot in the linear scale using the default values (w/200 data point).*
 
-The 2D (Anisotropic model) is based on the reference below where *I(q)* is approximated for 1d scattering. Thus the
-scattering pattern for 2D may not be accurate. Note that we are not responsible for any incorrectness of the 2D model
-computation.
+The 2D (Anisotropic model) is based on the reference below where *I(q)* is
+approximated for 1d scattering. Thus the scattering pattern for 2D may not
+be accurate. Note that we are not responsible for any incorrectness of the
+2D model computation.
 
 .. image:: img/image165.gif
@@ -78,5 +87,5 @@
 """
 
-from numpy import pi, inf
+from numpy import inf
 
 name = "fcc_paracrystal"
@@ -118,4 +127,3 @@
 # names and the target sasview model name.
 oldname='FCCrystalModel'
-oldpars=dict(sld='sldSph',
-             solvent_sld='sldSolv')
+oldpars=dict(sld='sldSph', solvent_sld='sldSolv')
Index: sasmodels/models/gaussian_peak.py
===================================================================
--- sasmodels/models/gaussian_peak.py	(revision a5d0d008778409dd90a86aa3db6037187ed3ff57)
+++ sasmodels/models/gaussian_peak.py	(revision 3c56da87f0995bdde77f4183b68e77b81584dde0)
@@ -21,5 +21,5 @@
 """
 
-from numpy import pi, inf
+from numpy import inf
 
 name = "gaussian_peak"
Index: sasmodels/models/hardsphere.py
===================================================================
--- sasmodels/models/hardsphere.py	(revision a5d0d008778409dd90a86aa3db6037187ed3ff57)
+++ sasmodels/models/hardsphere.py	(revision 3c56da87f0995bdde77f4183b68e77b81584dde0)
@@ -32,5 +32,5 @@
 """
 
-from numpy import pi, inf
+from numpy import inf
 
 name = "hardsphere"
Index: sasmodels/models/lamellar.py
===================================================================
--- sasmodels/models/lamellar.py	(revision bfb195e4f4fec755fb7d02be1dd464a4085e0bb3)
+++ sasmodels/models/lamellar.py	(revision 3c56da87f0995bdde77f4183b68e77b81584dde0)
@@ -47,5 +47,5 @@
 """
 
-from numpy import pi, inf
+from numpy import inf
 
 name = "lamellar"
Index: sasmodels/models/lamellarCaille.py
===================================================================
--- sasmodels/models/lamellarCaille.py	(revision bfb195e4f4fec755fb7d02be1dd464a4085e0bb3)
+++ sasmodels/models/lamellarCaille.py	(revision 3c56da87f0995bdde77f4183b68e77b81584dde0)
@@ -70,5 +70,5 @@
 also in J. Phys. Chem. B, 105, (2001) 11081-11088
 """
-from numpy import pi, inf
+from numpy import inf
 
 name = "lamellarPS"
Index: sasmodels/models/lamellarCailleHG.py
===================================================================
--- sasmodels/models/lamellarCailleHG.py	(revision 12c810f78dbf238faaad5e049615baec671280ae)
+++ sasmodels/models/lamellarCailleHG.py	(revision 3c56da87f0995bdde77f4183b68e77b81584dde0)
@@ -74,5 +74,5 @@
 also in J. Phys. Chem. B, 105, (2001) 11081-11088
 """
-from numpy import pi, inf
+from numpy import inf
 
 name = "lamellarCailleHG"
Index: sasmodels/models/lamellarFFHG.py
===================================================================
--- sasmodels/models/lamellarFFHG.py	(revision bfb195e4f4fec755fb7d02be1dd464a4085e0bb3)
+++ sasmodels/models/lamellarFFHG.py	(revision 3c56da87f0995bdde77f4183b68e77b81584dde0)
@@ -1,6 +1,7 @@
 # Note: model title and parameter table are inserted automatically
 r"""
-This model provides the scattering intensity, *I(q)*, for a lyotropic lamellar phase where a random distribution in
-solution are assumed. The SLD of the head region is taken to be different from the SLD of the tail region.
+This model provides the scattering intensity, *I(q)*, for a lyotropic lamellar
+phase where a random distribution in solution are assumed. The SLD of the head
+region is taken to be different from the SLD of the tail region.
 
 *2.1.31.1. Definition*
@@ -16,8 +17,10 @@
 .. image:: img/lamellarFFHG_.jpg
 
-where |delta|\ T = tail length (or *tail_length*), |delta|\ H = head thickness (or *h_thickness*),
-|drho|\ H = SLD(headgroup) - SLD(solvent), and |drho|\ T = SLD(tail) - SLD(solvent).
+where |delta|\ T = tail length (or *tail_length*), |delta|\ H = head thickness
+(or *h_thickness*), |drho|\ H = SLD(headgroup) - SLD(solvent),
+and |drho|\ T = SLD(tail) - SLD(solvent).
 
-The 2D scattering intensity is calculated in the same way as 1D, where the *q* vector is defined as
+The 2D scattering intensity is calculated in the same way as 1D, where
+the *q* vector is defined as
 
 .. math::
@@ -25,18 +28,7 @@
     Q = \sqrt{Q_x^2 + Q_y^2}
 
-The returned value is in units of |cm^-1|, on absolute scale. In the parameters, *sld_tail* = SLD of the tail group,
-and *sld_head* = SLD of the head group.
-
-==============  ========  =============
-Parameter name  Units     Default value
-==============  ========  =============
-background      |cm^-1|   0.0
-head_sld        |Ang^-2|  3e-06
-scale            None      1
-solvent_sld     |Ang^-2|  6e-06
-head_length     |Ang|     10
-tail_length     |Ang|     15
-sld (tail)      |Ang^-2|  0
-==============  ========  =============
+The returned value is in units of |cm^-1|, on absolute scale. In the
+parameters, *sld_tail* = SLD of the tail group, and *sld_head* = SLD
+of the head group.
 
 .. image:: img/lamellarFFHG_138.jpg
@@ -44,6 +36,6 @@
 *Figure. 1D plot using the default values (w/1000 data point).*
 
-Our model uses the form factor calculations implemented in a c-library provided by the NIST Center for Neutron Research
-(Kline, 2006).
+Our model uses the form factor calculations implemented in a C library
+provided by the NIST Center for Neutron Research (Kline, 2006).
 
 REFERENCE
@@ -56,5 +48,5 @@
 """
 
-from numpy import pi, inf
+from numpy import inf
 
 name = "lamellar_FFHG"
@@ -96,18 +88,20 @@
 Iq = """
     const double qsq = q*q;
- 	const double drh = head_sld - solvent_sld;
-	const double drt = sld	- solvent_sld;		//correction 13FEB06 by L.Porcar
-	const double qT = q*tail_length;
-	double Pq, inten;
-	Pq = drh*(sin(q*(head_length+tail_length))-sin(qT)) + drt*sin(qT);
-	Pq *= Pq;
-	Pq *= 4.0/(qsq);
-	
-	inten = 2.0e-4*M_PI*Pq/qsq;	
-	
-	return inten /= 2.0*(head_length+tail_length);	//normalize by the bilayer thickness	
-		
+    const double drh = head_sld - solvent_sld;
+    const double drt = sld - solvent_sld;    //correction 13FEB06 by L.Porcar
+    const double qT = q*tail_length;
+    double Pq, inten;
+    Pq = drh*(sin(q*(head_length+tail_length))-sin(qT)) + drt*sin(qT);
+    Pq *= Pq;
+    Pq *= 4.0/(qsq);
+
+    inten = 2.0e-4*M_PI*Pq/qsq;
+
+    // normalize by the bilayer thickness
+    inten /= 2.0*(head_length+tail_length);
+
+    return inten;
     """
-	
+
 Iqxy = """
     return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);
@@ -118,11 +112,13 @@
 
 demo = dict(
-        scale=1, background=0,
-		tail_length=15,head_length=10,
-        sld=0.4, head_sld=3.0, solvent_sld=6.0,
-		tail_length_pd= 0.2, tail_length_pd_n=40,
-		head_length_pd= 0.01, head_length_pd_n=40,
-          )
+    scale=1, background=0,
+    tail_length=15,head_length=10,
+    sld=0.4, head_sld=3.0, solvent_sld=6.0,
+    tail_length_pd= 0.2, tail_length_pd_n=40,
+    head_length_pd= 0.01, head_length_pd_n=40,
+    )
+
 oldname = 'LamellarFFHGModel'
-oldpars = dict(head_length='h_thickness', sld='sld_tail', head_sld='sld_head', solvent_sld='sld_solvent')
+oldpars = dict(head_length='h_thickness', sld='sld_tail',
+               head_sld='sld_head', solvent_sld='sld_solvent')
 
Index: sasmodels/models/lamellarPC.py
===================================================================
--- sasmodels/models/lamellarPC.py	(revision a5d0d008778409dd90a86aa3db6037187ed3ff57)
+++ sasmodels/models/lamellarPC.py	(revision 3c56da87f0995bdde77f4183b68e77b81584dde0)
@@ -1,7 +1,9 @@
 # Note: model title and parameter table are inserted automatically
 r"""
-This model calculates the scattering from a stack of repeating lamellar structures. The stacks of lamellae (infinite
-in lateral dimension) are treated as a paracrystal to account for the repeating spacing. The repeat distance is further
-characterized by a Gaussian polydispersity. **This model can be used for large multilamellar vesicles.**
+This model calculates the scattering from a stack of repeating lamellar
+structures. The stacks of lamellae (infinite in lateral dimension) are
+treated as a paracrystal to account for the repeating spacing. The repeat
+distance is further characterized by a Gaussian polydispersity. **This model
+can be used for large multilamellar vesicles.**
 
 *2.1.33.1. Definition*
@@ -11,18 +13,23 @@
 .. image:: img/image145.jpg
 
-The form factor of the bilayer is approximated as the cross section of an infinite, planar bilayer of thickness *t*
+The form factor of the bilayer is approximated as the cross section of an
+infinite, planar bilayer of thickness *t*
 
 .. image:: img/image146.jpg
 
-Here, the scale factor is used instead of the mass per area of the bilayer (*G*). The scale factor is the volume
-fraction of the material in the bilayer, *not* the total excluded volume of the paracrystal. *Z*\ :sub:`N`\ *(q)*
-describes the interference effects for aggregates consisting of more than one bilayer. The equations used are (3-5)
+Here, the scale factor is used instead of the mass per area of the
+bilayer (*G*). The scale factor is the volume fraction of the material in
+the bilayer, *not* the total excluded volume of the paracrystal.
+*Z*\ :sub:`N`\ *(q)* describes the interference effects for aggregates
+consisting of more than one bilayer. The equations used are (3-5)
 from the Bergstrom reference below.
 
-Non-integer numbers of stacks are calculated as a linear combination of the lower and higher values
+Non-integer numbers of stacks are calculated as a linear combination of
+the lower and higher values
 
 .. image:: img/image147.jpg
 
-The 2D scattering intensity is the same as 1D, regardless of the orientation of the *q* vector which is defined as
+The 2D scattering intensity is the same as 1D, regardless of the orientation
+of the *q* vector which is defined as
 
 .. math::
@@ -30,5 +37,6 @@
     Q = \sqrt{Q_x^2 + Q_y^2}
 
-The parameters of the model are *Nlayers* = no. of layers, and *pd_spacing* = polydispersity of spacing.
+The parameters of the model are *Nlayers* = no. of layers, and
+*pd_spacing* = polydispersity of spacing.
 
 ==============  ========  =============
@@ -49,14 +57,15 @@
 *Figure. 1D plot using the default values above (w/20000 data point).*
 
-Our model uses the form factor calculations implemented in a c-library provided by the NIST Center for Neutron Research
-(Kline, 2006).
+Our model uses the form factor calculations implemented in a C library
+provided by the NIST Center for Neutron Research (Kline, 2006).
 
 REFERENCE
 
-M Bergstrom, J S Pedersen, P Schurtenberger, S U Egelhaaf, *J. Phys. Chem. B*, 103 (1999) 9888-9897
+M Bergstrom, J S Pedersen, P Schurtenberger, S U Egelhaaf,
+*J. Phys. Chem. B*, 103 (1999) 9888-9897
 
 """
 
-from numpy import pi, inf
+from numpy import inf
 
 name = "lamellarPC"
@@ -90,5 +99,5 @@
     ]
 
-	
+
 source = [ "lamellarPC_kernel.c"]
 
@@ -105,11 +114,15 @@
 
 demo = dict(
-        scale=1, background=0,
-		thickness=33,Nlayers=20,spacing=250,spacing_polydisp=0.2,
-        sld=1.0, solvent_sld=6.34,
-		thickness_pd= 0.2, thickness_pd_n=40
-         )
+    scale=1, background=0,
+    thickness=33, Nlayers=20, spacing=250, spacing_polydisp=0.2,
+    sld=1.0, solvent_sld=6.34,
+    thickness_pd= 0.2, thickness_pd_n=40
+    )
+
 oldname = 'LamellarPCrystalModel'
-oldpars = dict(spacing_polydisp='pd_spacing', sld='sld_layer',solvent_sld='sld_solvent')
+oldpars = dict(
+    spacing_polydisp='pd_spacing', sld='sld_layer',
+    solvent_sld='sld_solvent'
+    )
 
 
Index: sasmodels/models/sphere.py
===================================================================
--- sasmodels/models/sphere.py	(revision a5d0d008778409dd90a86aa3db6037187ed3ff57)
+++ sasmodels/models/sphere.py	(revision 3c56da87f0995bdde77f4183b68e77b81584dde0)
@@ -58,5 +58,5 @@
 """
 
-from numpy import pi, inf
+from numpy import inf
 
 name = "sphere"
Index: sasmodels/models/spherepy.py
===================================================================
--- sasmodels/models/spherepy.py	(revision a5d0d008778409dd90a86aa3db6037187ed3ff57)
+++ sasmodels/models/spherepy.py	(revision 3c56da87f0995bdde77f4183b68e77b81584dde0)
@@ -59,5 +59,5 @@
 
 import numpy as np
-from numpy import pi, inf, sin, cos, sqrt, exp, log
+from numpy import pi, inf, sin, cos, sqrt, log
 
 name = "sphere"
Index: sasmodels/models/stickyhardsphere.py
===================================================================
--- sasmodels/models/stickyhardsphere.py	(revision 7e224c2f323be37a409d509cfa32db53733a321a)
+++ sasmodels/models/stickyhardsphere.py	(revision 3c56da87f0995bdde77f4183b68e77b81584dde0)
@@ -1,16 +1,19 @@
 # Note: model title and parameter table are inserted automatically
-r"""This calculates the interparticle structure factor for a hard sphere fluid with
-a narrow attractive well. A perturbative solution of the Percus-Yevick closure is used.
-The strength of the attractive well is described in terms of "stickiness"
-as defined below. The returned value is a dimensionless structure factor, *S(q)*.
+r"""
+This calculates the interparticle structure factor for a hard sphere fluid
+with a narrow attractive well. A perturbative solution of the Percus-Yevick
+closure is used. The strength of the attractive well is described in terms
+of "stickiness" as defined below. The returned value is a dimensionless
+structure factor, *S(q)*.
 
-The perturb (perturbation parameter), |epsilon|, should be held between 0.01 and 0.1.
-It is best to hold the perturbation parameter fixed and let the "stickiness" vary to
-adjust the interaction strength. The stickiness, |tau|, is defined in the equation
-below and is a function of both the perturbation parameter and the interaction strength.
-|tau| and |epsilon| are defined in terms of the hard sphere diameter (|sigma| = 2\*\ *R*\ ),
-the width of the square well, |bigdelta| (same units as *R*), and the depth of the well,
-*Uo*, in units of kT. From the definition, it is clear that smaller |tau| means stronger
-attraction.
+The perturb (perturbation parameter), |epsilon|, should be held between 0.01
+and 0.1. It is best to hold the perturbation parameter fixed and let
+the "stickiness" vary to adjust the interaction strength. The stickiness,
+|tau|, is defined in the equation below and is a function of both the
+perturbation parameter and the interaction strength. |tau| and |epsilon|
+are defined in terms of the hard sphere diameter (|sigma| = 2\*\ *R*\ ), the
+width of the square well, |bigdelta| (same units as *R*), and the depth of
+the well, *Uo*, in units of kT. From the definition, it is clear that
+smaller |tau| means stronger attraction.
 
 .. image:: img/stickyhardsphere_228.PNG
@@ -20,24 +23,27 @@
 .. image:: img/stickyhardsphere_229.PNG
 
-The Percus-Yevick (PY) closure was used for this calculation, and is an adequate closure
-for an attractive interparticle potential. This solution has been compared to Monte Carlo
-simulations for a square well fluid, with good agreement.
+The Percus-Yevick (PY) closure was used for this calculation, and is an
+adequate closure for an attractive interparticle potential. This solution
+has been compared to Monte Carlo simulations for a square well fluid, with
+good agreement.
 
-The true particle volume fraction, |phi|, is not equal to *h*, which appears in most of
-the reference. The two are related in equation (24) of the reference. The reference also
-describes the relationship between this perturbation solution and the original sticky hard
-sphere (or adhesive sphere) model by Baxter.
+The true particle volume fraction, |phi|, is not equal to *h*, which appears
+in most of the reference. The two are related in equation (24) of the
+reference. The reference also describes the relationship between this
+perturbation solution and the original sticky hard sphere (or adhesive
+sphere) model by Baxter.
 
-NB: The calculation can go haywire for certain combinations of the input parameters,
-producing unphysical solutions - in this case errors are reported to the command window and
-the *S(q)* is set to -1 (so it will disappear on a log-log plot). Use tight bounds to keep
-the parameters to values that you know are physical (test them) and keep nudging them until
+NB: The calculation can go haywire for certain combinations of the input
+parameters, producing unphysical solutions - in this case errors are
+reported to the command window and the *S(q)* is set to -1 (so it will
+disappear on a log-log plot). Use tight bounds to keep the parameters to
+values that you know are physical (test them) and keep nudging them until
 the optimization does not hit the constraints.
 
-In sasview the effective radius will be calculated from the parameters used in the
-form factor P(Q) that this S(Q) is combined with.
+In sasview the effective radius will be calculated from the parameters
+used in the form factor P(Q) that this S(Q) is combined with.
 
-For 2D data: The 2D scattering intensity is calculated in the same way as 1D, where
-the *q* vector is defined as
+For 2D data: The 2D scattering intensity is calculated in the same way
+as 1D, where the *q* vector is defined as
 
 .. math::
@@ -99,67 +105,67 @@
 
 Iq = """
-	double onemineps,eta;
-	double sig,aa,etam1,etam1sq,qa,qb,qc,radic;
-	double lam,lam2,test,mu,alpha,beta;
-	double kk,k2,k3,ds,dc,aq1,aq2,aq3,aq,bq1,bq2,bq3,bq,sq;
+    double onemineps,eta;
+    double sig,aa,etam1,etam1sq,qa,qb,qc,radic;
+    double lam,lam2,test,mu,alpha,beta;
+    double kk,k2,k3,ds,dc,aq1,aq2,aq3,aq,bq1,bq2,bq3,bq,sq;
 
-	onemineps = 1.0-perturb;
-	eta = volfraction/onemineps/onemineps/onemineps;
+    onemineps = 1.0-perturb;
+    eta = volfraction/onemineps/onemineps/onemineps;
 
-	sig = 2.0 * effect_radius;
-	aa = sig/onemineps;
-	etam1 = 1.0 - eta;
-	etam1sq=etam1*etam1;
-	//C
-	//C  SOLVE QUADRATIC FOR LAMBDA
-	//C
-	qa = eta/12.0;
-	qb = -1.0*(stickiness + eta/etam1);
-	qc = (1.0 + eta/2.0)/etam1sq;
-	radic = qb*qb - 4.0*qa*qc;
-	if(radic<0) {
-		//if(x>0.01 && x<0.015)
-	 	//	Print "Lambda unphysical - both roots imaginary"
-	 	//endif
-		return(-1.0);
-	}
-	//C   KEEP THE SMALLER ROOT, THE LARGER ONE IS UNPHYSICAL
-	lam = (-1.0*qb-sqrt(radic))/(2.0*qa);
-	lam2 = (-1.0*qb+sqrt(radic))/(2.0*qa);
-	if(lam2<lam) {
-		lam = lam2;
-	}
-	test = 1.0 + 2.0*eta;
-	mu = lam*eta*etam1;
-	if(mu>test) {
-		//if(x>0.01 && x<0.015)
-		// Print "Lambda unphysical mu>test"
-		//endif
-		return(-1.0);
-	}
-	alpha = (1.0 + 2.0*eta - mu)/etam1sq;
-	beta = (mu - 3.0*eta)/(2.0*etam1sq);
-	//C
-	//C   CALCULATE THE STRUCTURE FACTOR
-	//C
-	kk = q*aa;
-	k2 = kk*kk;
-	k3 = kk*k2;
-	SINCOS(kk,ds,dc);
-	//ds = sin(kk);
-	//dc = cos(kk);
-	aq1 = ((ds - kk*dc)*alpha)/k3;
-	aq2 = (beta*(1.0-dc))/k2;
-	aq3 = (lam*ds)/(12.0*kk);
-	aq = 1.0 + 12.0*eta*(aq1+aq2-aq3);
-	//
-	bq1 = alpha*(0.5/kk - ds/k2 + (1.0 - dc)/k3);
-	bq2 = beta*(1.0/kk - ds/k2);
-	bq3 = (lam/12.0)*((1.0 - dc)/kk);
-	bq = 12.0*eta*(bq1+bq2-bq3);
-	//
-	sq = 1.0/(aq*aq +bq*bq);
+    sig = 2.0 * effect_radius;
+    aa = sig/onemineps;
+    etam1 = 1.0 - eta;
+    etam1sq=etam1*etam1;
+    //C
+    //C  SOLVE QUADRATIC FOR LAMBDA
+    //C
+    qa = eta/12.0;
+    qb = -1.0*(stickiness + eta/etam1);
+    qc = (1.0 + eta/2.0)/etam1sq;
+    radic = qb*qb - 4.0*qa*qc;
+    if(radic<0) {
+        //if(x>0.01 && x<0.015)
+        //	Print "Lambda unphysical - both roots imaginary"
+        //endif
+        return(-1.0);
+    }
+    //C   KEEP THE SMALLER ROOT, THE LARGER ONE IS UNPHYSICAL
+    lam = (-1.0*qb-sqrt(radic))/(2.0*qa);
+    lam2 = (-1.0*qb+sqrt(radic))/(2.0*qa);
+    if(lam2<lam) {
+        lam = lam2;
+    }
+    test = 1.0 + 2.0*eta;
+    mu = lam*eta*etam1;
+    if(mu>test) {
+        //if(x>0.01 && x<0.015)
+        // Print "Lambda unphysical mu>test"
+        //endif
+        return(-1.0);
+    }
+    alpha = (1.0 + 2.0*eta - mu)/etam1sq;
+    beta = (mu - 3.0*eta)/(2.0*etam1sq);
+    //C
+    //C   CALCULATE THE STRUCTURE FACTOR
+    //C
+    kk = q*aa;
+    k2 = kk*kk;
+    k3 = kk*k2;
+    SINCOS(kk,ds,dc);
+    //ds = sin(kk);
+    //dc = cos(kk);
+    aq1 = ((ds - kk*dc)*alpha)/k3;
+    aq2 = (beta*(1.0-dc))/k2;
+    aq3 = (lam*ds)/(12.0*kk);
+    aq = 1.0 + 12.0*eta*(aq1+aq2-aq3);
+    //
+    bq1 = alpha*(0.5/kk - ds/k2 + (1.0 - dc)/k3);
+    bq2 = beta*(1.0/kk - ds/k2);
+    bq3 = (lam/12.0)*((1.0 - dc)/kk);
+    bq = 12.0*eta*(bq1+bq2-bq3);
+    //
+    sq = 1.0/(aq*aq +bq*bq);
 
-	return(sq);
+    return(sq);
 """
 
@@ -176,4 +182,2 @@
             stickiness=0.2, effect_radius_pd=0.1, effect_radius_pd_n=40)
 
-
-
Index: sasmodels/models/triaxial_ellipsoid.py
===================================================================
--- sasmodels/models/triaxial_ellipsoid.py	(revision a5d0d008778409dd90a86aa3db6037187ed3ff57)
+++ sasmodels/models/triaxial_ellipsoid.py	(revision 3c56da87f0995bdde77f4183b68e77b81584dde0)
@@ -90,5 +90,5 @@
 """
 
-from numpy import pi, inf
+from numpy import inf
 
 name = "triaxial_ellipsoid"
Index: sasmodels/sasview_model.py
===================================================================
--- sasmodels/sasview_model.py	(revision de0c4baaebb72d2cb7c4db56d90e0fc508888f7d)
+++ sasmodels/sasview_model.py	(revision 3c56da87f0995bdde77f4183b68e77b81584dde0)
@@ -69,5 +69,5 @@
         self.dispersion = dict()
         partype = model.info['partype']
-        for name, units, default, limits, ptype, description in model.info['parameters']:
+        for name, units, default, limits, _, _ in model.info['parameters']:
             self.params[name] = default
             self.details[name] = [units] + limits
@@ -120,4 +120,5 @@
 
 
+    # pylint: disable=no-self-use
     def getProfile(self):
         """
@@ -213,4 +214,5 @@
         """
         if isinstance(x, (list, tuple)):
+            # pylint: disable=unpacking-non-sequence
             q, phi = x
             return self.calculate_Iq([q * math.cos(phi)],
@@ -263,5 +265,6 @@
 
 
-        :param qdist: ndarray of scalar q-values or list [qx,qy] where qx,qy are 1D ndarrays
+        :param qdist: ndarray of scalar q-values or list [qx,qy]
+        where qx,qy are 1D ndarrays
         """
         if isinstance(qdist, (list, tuple)):
@@ -279,5 +282,6 @@
 
         else:
-            raise TypeError("evalDistribution expects q or [qx, qy], not %r" % type(qdist))
+            raise TypeError("evalDistribution expects q or [qx, qy], not %r"
+                            % type(qdist))
 
     def calculate_Iq(self, *args):
@@ -382,7 +386,7 @@
         limits = self._model.info['limits']
         dis = self.dispersion[par]
-        v, w = weights.get_weights(
+        value, weight = weights.get_weights(
             dis['type'], dis['npts'], dis['width'], dis['nsigmas'],
             self.params[par], limits[par], par in relative)
-        return v, w / w.max()
-
+        return value, weight / np.sum(weight)
+
Index: sasmodels/sesans.py
===================================================================
--- sasmodels/sesans.py	(revision c97724e7592a12c8cc1d1da48b84c88fb8a97093)
+++ sasmodels/sesans.py	(revision 3c56da87f0995bdde77f4183b68e77b81584dde0)
@@ -22,13 +22,13 @@
 # TODO: dead code; for now the call to the hankel transform happens in BumpsModel
 class SesansCalculator:
-    def __init__(self, sans_kernel, q_zmax, Rmax, SElength, wavelength, thickness):
-        self._set_kernel(sans_kernel, q_zmax, Rmax)
+    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, sans_kernel, q_zmax, Rmax):
-        input = sans_kernel.make_input([make_q(q_zmax, Rmax)])
-        self.sans_calculator = sans_kernel(input)
+    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):
Index: sasmodels/weights.py
===================================================================
--- sasmodels/weights.py	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
+++ sasmodels/weights.py	(revision 3c56da87f0995bdde77f4183b68e77b81584dde0)
@@ -26,4 +26,5 @@
         return pars
 
+    # pylint: disable=no-self-use
     def set_weights(self, values, weights):
         raise RuntimeError("set_weights is only available for ArrayDispersion")
