Index: compare-new.py
===================================================================
--- compare-new.py	(revision 13d86bcf0b1f0ccf1dadb231e2eb7c390495a5e4)
+++ compare-new.py	(revision ff7119b4ad066090cdffae80706883cef3d80399)
@@ -6,5 +6,6 @@
 import numpy as np
 
-from sasmodels.core import BumpsModel, plot_data, tic, opencl_model, dll_model
+from sasmodels.bumps_model import BumpsModel, plot_data, tic
+from sasmodels.gen import opencl_model, dll_model
 
 def sasview_model(modelname, **pars):
Index: sasmodels/alignment.py
===================================================================
--- sasmodels/alignment.py	(revision 14de349c663ca9521843af50415c37aed27a9eef)
+++ sasmodels/alignment.py	(revision ff7119b4ad066090cdffae80706883cef3d80399)
@@ -3,12 +3,13 @@
 
 Some web sites say that maximizing performance for OpenCL code requires
-aligning data on certain memory boundaries.
+aligning data on certain memory boundaries.  The following functions
+provide this service:
 
-:func:`data_alignment` queries all devices in the OpenCL context, returning
-the most restrictive alignment.
+:func:`align_data` aligns an existing array, returning a new array of the
+correct alignment.
 
-:func:`align_data` aligns an existing array.
+:func:`align_empty` to create an empty array of the correct alignment.
 
-:func:`align_empty` to create a new array of the correct alignment.
+Set alignment to :func:`gpu.environment()` attribute *boundary*.
 
 Note:  This code is unused. So far, tests have not demonstrated any
@@ -17,16 +18,10 @@
 to decide whether it is really required.
 """
-
 import numpy as np
-import pyopencl as cl
-
-def data_alignment(context):
-    """
-    Return the desired byte alignment across all devices.
-    """
-    # Note: rely on the fact that alignment will be 2^k
-    return max(d.min_data_type_align_size for d in context.devices)
 
 def align_empty(shape, dtype, alignment=128):
+    """
+    Return an empty array aligned on the alignment boundary.
+    """
     size = np.prod(shape)
     dtype = np.dtype(dtype)
@@ -40,4 +35,7 @@
 
 def align_data(v, 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)
@@ -45,10 +43,2 @@
     return view
 
-def work_group_size(context, kernel):
-    """
-    Return the desired work group size for the context.
-    """
-    max(kernel.get_work_group_info(cl.kernel_work_group_info.PREFERRED_WORK_GROUP_SIZE_MULTIPLE, d)
-        for d in context.devices)
-
-
Index: sasmodels/bumps_model.py
===================================================================
--- sasmodels/bumps_model.py	(revision ff7119b4ad066090cdffae80706883cef3d80399)
+++ sasmodels/bumps_model.py	(revision ff7119b4ad066090cdffae80706883cef3d80399)
@@ -0,0 +1,336 @@
+"""
+Sasmodels core.
+"""
+import sys, os
+import datetime
+
+import numpy as np
+
+
+def tic():
+    """
+    Timer function.
+
+    Use "toc=tic()" to start the clock and "toc()" to measure
+    a time interval.
+    """
+    then = datetime.datetime.now()
+    return lambda: (datetime.datetime.now()-then).total_seconds()
+
+
+def load_data(filename):
+    """
+    Load data using a sasview loader.
+    """
+    from sans.dataloader.loader import Loader
+    loader = Loader()
+    data = loader.load(filename)
+    if data is None:
+        raise IOError("Data %r could not be loaded"%filename)
+    return data
+
+
+def empty_data1D(q):
+    """
+    Create empty 1D data using the given *q* as the x value.
+
+    Resolutions dq/q is 5%.
+    """
+
+    from sans.dataloader.data_info import Data1D
+
+    Iq = 100*np.ones_like(q)
+    dIq = np.sqrt(Iq)
+    data = Data1D(q, Iq, dx=0.05*q, dy=dIq)
+    data.filename = "fake data"
+    data.qmin, data.qmax = q.min(), q.max()
+    return data
+
+
+def empty_data2D(qx, qy=None):
+    """
+    Create empty 2D data using the given mesh.
+
+    If *qy* is missing, create a square mesh with *qy=qx*.
+
+    Resolution dq/q is 5%.
+    """
+    from sans.dataloader.data_info import Data2D, Detector
+
+    if qy is None:
+        qy = qx
+    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
+
+    # 5% dQ/Q resolution
+    data.dqx_data = 0.05*Qx
+    data.dqy_data = 0.05*Qy
+
+    detector = Detector()
+    detector.pixel_size.x = 5 # mm
+    detector.pixel_size.y = 5 # mm
+    detector.distance = 4 # m
+    data.detector.append(detector)
+    data.xbins = qx
+    data.ybins = qy
+    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
+
+
+def set_beam_stop(data, radius, outer=None):
+    """
+    Add a beam stop of the given *radius*.  If *outer*, make an annulus.
+    """
+    from sans.dataloader.manipulations import Ringcut
+    if hasattr(data, 'qx_data'):
+        data.mask = Ringcut(0, radius)(data)
+        if outer is not None:
+            data.mask += Ringcut(outer,np.inf)(data)
+    else:
+        data.mask = (data.x>=radius)
+        if outer is not None:
+            data.mask &= (data.x<outer)
+
+
+def set_half(data, half):
+    """
+    Select half of the data, either "right" or "left".
+    """
+    from sans.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)
+    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*.
+    """
+    from sans.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, scale='log'):
+    """
+    Plot the target value for the data.  This could be the data itself,
+    the theory calculation, or the residuals.
+
+    *scale* can be 'log' for log scale data, or 'linear'.
+    """
+    from numpy.ma import masked_array, masked
+    import matplotlib.pyplot as plt
+    if hasattr(data, 'qx_data'):
+        img = masked_array(iq, data.mask)
+        if scale == 'log':
+            img[(img <= 0) | ~np.isfinite(img)] = masked
+            img = np.log10(img)
+        xmin, xmax = min(data.qx_data), max(data.qx_data)
+        ymin, ymax = min(data.qy_data), max(data.qy_data)
+        plt.imshow(img.reshape(128,128),
+                   interpolation='nearest', aspect=1, origin='upper',
+                   extent=[xmin, xmax, ymin, ymax], vmin=vmin, vmax=vmax)
+    else: # 1D data
+        if scale == 'linear':
+            idx = np.isfinite(iq)
+            plt.plot(data.x[idx], iq[idx])
+        else:
+            idx = np.isfinite(iq) & (iq>0)
+            plt.loglog(data.x[idx], iq[idx])
+
+
+def _plot_result1D(data, theory, view):
+    """
+    Plot the data and residuals for 1D data.
+    """
+    import matplotlib.pyplot as plt
+    from numpy.ma import masked_array, masked
+    #print "not a number",sum(np.isnan(data.y))
+    #data.y[data.y<0.05] = 0.5
+    mdata = masked_array(data.y, data.mask)
+    mdata[np.isnan(mdata)] = masked
+    if view is 'log':
+        mdata[mdata <= 0] = masked
+    mtheory = masked_array(theory, mdata.mask)
+    mresid = masked_array((theory-data.y)/data.dy, mdata.mask)
+
+    plt.subplot(121)
+    plt.errorbar(data.x, mdata, yerr=data.dy)
+    plt.plot(data.x, mtheory, '-', hold=True)
+    plt.yscale(view)
+    plt.subplot(122)
+    plt.plot(data.x, mresid, 'x')
+    #plt.axhline(1, color='black', ls='--',lw=1, hold=True)
+    #plt.axhline(0, color='black', lw=1, hold=True)
+    #plt.axhline(-1, color='black', ls='--',lw=1, hold=True)
+
+
+def _plot_result2D(data, theory, view):
+    """
+    Plot the data and residuals for 2D data.
+    """
+    import matplotlib.pyplot as plt
+    resid = (theory-data.data)/data.err_data
+    plt.subplot(131)
+    plot_data(data, data.data, scale=view)
+    plt.colorbar()
+    plt.subplot(132)
+    plot_data(data, theory, scale=view)
+    plt.colorbar()
+    plt.subplot(133)
+    plot_data(data, resid, scale='linear')
+    plt.colorbar()
+
+def plot_result(data, theory, view='log'):
+    """
+    Plot the data and residuals.
+    """
+    if hasattr(data, 'qx_data'):
+        _plot_result2D(data, theory, view)
+    else:
+        _plot_result1D(data, theory, view)
+
+
+class BumpsModel(object):
+    """
+    Return a bumps wrapper for a SAS model.
+
+    *data* is the data to be fitted.
+
+    *model* is the SAS model, e.g., from :func:`gen.opencl_model`.
+
+    *cutoff* is the integration cutoff, which avoids computing the
+    the SAS model where the polydispersity weight is low.
+
+    Model parameters can be initialized with additional keyword
+    arguments, or by assigning to model.parameter_name.value.
+
+    The resulting bumps model can be used directly in a FitProblem call.
+    """
+    def __init__(self, data, model, cutoff=1e-5, **kw):
+        from bumps.names import Parameter
+
+        # interpret data
+        self.data = data
+        if hasattr(data, 'qx_data'):
+            self.index = (data.mask==0) & (~np.isnan(data.data))
+            self.iq = data.data[self.index]
+            self.diq = data.err_data[self.index]
+            self._theory = np.zeros_like(data.data)
+            q_vectors = [data.qx_data, data.qy_data]
+        else:
+            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._theory = np.zeros_like(data.y)
+            q_vectors = [data.x]
+        #input = model.make_input(q_vectors)
+        input = model.make_input([v[self.index] for v in q_vectors])
+
+        # create model
+        self.fn = model(input)
+        self.cutoff = cutoff
+
+        # define bumps parameters
+        pars = []
+        extras = []
+        for p in model.info['parameters']:
+            name, default, limits, ptype = p[0], p[2], p[3], p[4]
+            value = kw.pop(name, default)
+            setattr(self, name, Parameter.default(value, name=name, limits=limits))
+            pars.append(name)
+        for name in model.info['partype']['pd-2d']:
+            for xpart,xdefault,xlimits in [
+                    ('_pd', 0, limits),
+                    ('_pd_n', 35, (0,1000)),
+                    ('_pd_nsigma', 3, (0, 10)),
+                    ('_pd_type', 'gaussian', None),
+                ]:
+                xname = name+xpart
+                xvalue = kw.pop(xname, xdefault)
+                if xlimits is not None:
+                    xvalue = Parameter.default(xvalue, name=xname, limits=xlimits)
+                    pars.append(xname)
+                setattr(self, xname, xvalue)
+        self._parameter_names = pars
+        if kw:
+            raise TypeError("unexpected parameters: %s"%(", ".join(sorted(kw.keys()))))
+        self.update()
+
+    def update(self):
+        self._cache = {}
+
+    def numpoints(self):
+        return len(self.iq)
+
+    def parameters(self):
+        return dict((k,getattr(self,k)) for k in self._parameter_names)
+
+    def theory(self):
+        if 'theory' not in self._cache:
+            pars = [getattr(self,p).value for p in self.fn.fixed_pars]
+            pd_pars = [self._get_weights(p) for p in self.fn.pd_pars]
+            #print pars
+            self._theory[self.index] = self.fn(pars, pd_pars, self.cutoff)
+            #self._theory[:] = self.fn.eval(pars, pd_pars)
+            self._cache['theory'] = self._theory
+        return self._cache['theory']
+
+    def residuals(self):
+        #if np.any(self.err ==0): print "zeros in err"
+        return (self.theory()[self.index]-self.iq)/self.diq
+
+    def nllf(self):
+        R = 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
+
+    def plot(self, view='log'):
+        plot_result(self.data, self.theory(), view=view)
+
+    def save(self, basename):
+        pass
+
+    def _get_weights(self, par):
+        from . import weights
+
+        relative = self.fn.info['partype']['pd-rel']
+        limits = self.fn.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, int(npts.value), width.value, nsigma.value,
+            value.value, limits[par], par in relative)
+        return v,w/w.max()
+
+
+def demo():
+    data = load_data('DEC07086.DAT')
+    set_beam_stop(data, 0.004)
+    plot_data(data)
+    import matplotlib.pyplot as plt; plt.show()
+
+
+if __name__ == "__main__":
+    demo()
Index: smodels/core.py
===================================================================
--- sasmodels/core.py	(revision 13d86bcf0b1f0ccf1dadb231e2eb7c390495a5e4)
+++ 	(revision )
@@ -1,325 +1,0 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
-
-import sys, os
-import datetime
-import warnings
-
-import numpy as np
-
-
-def opencl_model(kernel_module, dtype="single"):
-    from sasmodels import gen, gpu
-
-    source, info = gen.make(kernel_module)
-    ## for debugging, save source to a .cl file, edit it, and reload as model
-    #open(modelname+'.cl','w').write(source)
-    #source = open(modelname+'.cl','r').read()
-    return gpu.GpuModel(source, info, dtype)
-
-
-if sys.platform == 'darwin':
-    COMPILE = "gcc-mp-4.7 -shared -fPIC -std=c99 -fopenmp -O2 -Wall %s -o %s -lm -lgomp"
-elif os.name == 'nt':
-    COMPILE = "gcc -shared -fPIC -std=c99 -fopenmp -O2 -Wall %s -o %s -lm"
-else:
-    COMPILE = "cc -shared -fPIC -std=c99 -fopenmp -O2 -Wall %s -o %s -lm"
-DLL_PATH = "/tmp"
-
-
-def dll_path(info):
-    from os.path import join as joinpath, split as splitpath, splitext
-    basename = splitext(splitpath(info['filename'])[1])[0]
-    return joinpath(DLL_PATH, basename+'.so')
-
-
-def dll_model(kernel_module):
-    import os
-    import tempfile
-    from sasmodels import gen, dll
-
-    source, info = gen.make(kernel_module)
-    dllpath = dll_path(info)
-    if not os.path.exists(dllpath) \
-            or (os.path.getmtime(dllpath) < os.path.getmtime(info['filename'])):
-        # Replace with a proper temp file
-        srcfile = tempfile.mkstemp(suffix=".c",prefix="sas_"+info['name'])
-        open(srcfile, 'w').write(source)
-        os.system(COMPILE%(srcfile, dllpath))
-        ## comment the following to keep the generated c file
-        #os.unlink(srcfile)
-    return dll.DllModel(dllpath, info)
-
-
-TIC = None
-def tic():
-    global TIC
-    then = datetime.datetime.now()
-    TIC = lambda: (datetime.datetime.now()-then).total_seconds()
-    return TIC
-
-
-def toc():
-    return TIC()
-
-
-def load_data(filename):
-    from sans.dataloader.loader import Loader
-    loader = Loader()
-    data = loader.load(filename)
-    if data is None:
-        raise IOError("Data %r could not be loaded"%filename)
-    return data
-
-
-def empty_data2D(qx, qy=None):
-    from sans.dataloader.data_info import Data2D, Detector
-
-    if qy is None:
-        qy = qx
-    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
-
-    # 5% dQ/Q resolution
-    data.dqx_data = 0.05*Qx
-    data.dqy_data = 0.05*Qy
-
-    detector = Detector()
-    detector.pixel_size.x = 5 # mm
-    detector.pixel_size.y = 5 # mm
-    detector.distance = 4 # m
-    data.detector.append(detector)
-    data.xbins = qx
-    data.ybins = qy
-    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
-
-
-def empty_data1D(q):
-    from sans.dataloader.data_info import Data1D
-
-    Iq = 100*np.ones_like(q)
-    dIq = np.sqrt(Iq)
-    data = Data1D(q, Iq, dx=0.05*q, dy=dIq)
-    data.filename = "fake data"
-    data.qmin, data.qmax = q.min(), q.max()
-    return data
-
-
-def set_beam_stop(data, radius, outer=None):
-    from sans.dataloader.manipulations import Ringcut
-    if hasattr(data, 'qx_data'):
-        data.mask = Ringcut(0, radius)(data)
-        if outer is not None:
-            data.mask += Ringcut(outer,np.inf)(data)
-    else:
-        data.mask = (data.x>=radius)
-        if outer is not None:
-            data.mask &= (data.x<outer)
-
-
-def set_half(data, half):
-    from sans.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)
-    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):
-    from sans.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, scale='log'):
-    from numpy.ma import masked_array, masked
-    import matplotlib.pyplot as plt
-    if hasattr(data, 'qx_data'):
-        img = masked_array(iq, data.mask)
-        if scale == 'log':
-            img[(img <= 0) | ~np.isfinite(img)] = masked
-            img = np.log10(img)
-        xmin, xmax = min(data.qx_data), max(data.qx_data)
-        ymin, ymax = min(data.qy_data), max(data.qy_data)
-        plt.imshow(img.reshape(128,128),
-                   interpolation='nearest', aspect=1, origin='upper',
-                   extent=[xmin, xmax, ymin, ymax], vmin=vmin, vmax=vmax)
-    else: # 1D data
-        if scale == 'linear':
-            idx = np.isfinite(iq)
-            plt.plot(data.x[idx], iq[idx])
-        else:
-            idx = np.isfinite(iq) & (iq>0)
-            plt.loglog(data.x[idx], iq[idx])
-
-
-def plot_result2D(data, theory, view='log'):
-    import matplotlib.pyplot as plt
-    resid = (theory-data.data)/data.err_data
-    plt.subplot(131)
-    plot_data(data, data.data, scale=view)
-    plt.colorbar()
-    plt.subplot(132)
-    plot_data(data, theory, scale=view)
-    plt.colorbar()
-    plt.subplot(133)
-    plot_data(data, resid, scale='linear')
-    plt.colorbar()
-
-
-def plot_result1D(data, theory, view='log'):
-    import matplotlib.pyplot as plt
-    from numpy.ma import masked_array, masked
-    #print "not a number",sum(np.isnan(data.y))
-    #data.y[data.y<0.05] = 0.5
-    mdata = masked_array(data.y, data.mask)
-    mdata[np.isnan(mdata)] = masked
-    if view is 'log':
-        mdata[mdata <= 0] = masked
-    mtheory = masked_array(theory, mdata.mask)
-    mresid = masked_array((theory-data.y)/data.dy, mdata.mask)
-
-    plt.subplot(121)
-    plt.errorbar(data.x, mdata, yerr=data.dy)
-    plt.plot(data.x, mtheory, '-', hold=True)
-    plt.yscale(view)
-    plt.subplot(122)
-    plt.plot(data.x, mresid, 'x')
-    #plt.axhline(1, color='black', ls='--',lw=1, hold=True)
-    #plt.axhline(0, color='black', lw=1, hold=True)
-    #plt.axhline(-1, color='black', ls='--',lw=1, hold=True)
-
-
-class BumpsModel(object):
-    def __init__(self, data, model, cutoff=1e-5, **kw):
-        from bumps.names import Parameter
-        from . import gpu
-
-        # interpret data
-        self.is2D = hasattr(data,'qx_data')
-        self.data = data
-        if self.is2D:
-            self.index = (data.mask==0) & (~np.isnan(data.data))
-            self.iq = data.data[self.index]
-            self.diq = data.err_data[self.index]
-            self._theory = np.zeros_like(data.data)
-            q_vectors = [data.qx_data, data.qy_data]
-        else:
-            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._theory = np.zeros_like(data.y)
-            q_vectors = [data.x]
-        #input = model.make_input(q_vectors)
-        input = model.make_input([v[self.index] for v in q_vectors])
-
-        # create model
-        self.fn = model(input)
-        self.cutoff = cutoff
-
-        # define bumps parameters
-        pars = []
-        extras = []
-        for p in model.info['parameters']:
-            name, default, limits, ptype = p[0], p[2], p[3], p[4]
-            value = kw.pop(name, default)
-            setattr(self, name, Parameter.default(value, name=name, limits=limits))
-            pars.append(name)
-        for name in model.info['partype']['pd-2d']:
-            for xpart,xdefault,xlimits in [
-                    ('_pd', 0, limits),
-                    ('_pd_n', 35, (0,1000)),
-                    ('_pd_nsigma', 3, (0, 10)),
-                    ('_pd_type', 'gaussian', None),
-                ]:
-                xname = name+xpart
-                xvalue = kw.pop(xname, xdefault)
-                if xlimits is not None:
-                    xvalue = Parameter.default(xvalue, name=xname, limits=xlimits)
-                    pars.append(xname)
-                setattr(self, xname, xvalue)
-        self._parameter_names = pars
-        if kw:
-            raise TypeError("unexpected parameters: %s"%(", ".join(sorted(kw.keys()))))
-        self.update()
-
-    def update(self):
-        self._cache = {}
-
-    def numpoints(self):
-        return len(self.iq)
-
-    def parameters(self):
-        return dict((k,getattr(self,k)) for k in self._parameter_names)
-
-    def theory(self):
-        if 'theory' not in self._cache:
-            pars = [getattr(self,p).value for p in self.fn.fixed_pars]
-            pd_pars = [self._get_weights(p) for p in self.fn.pd_pars]
-            #print pars
-            self._theory[self.index] = self.fn(pars, pd_pars, self.cutoff)
-            #self._theory[:] = self.fn.eval(pars, pd_pars)
-            self._cache['theory'] = self._theory
-        return self._cache['theory']
-
-    def residuals(self):
-        #if np.any(self.err ==0): print "zeros in err"
-        return (self.theory()[self.index]-self.iq)/self.diq
-
-    def nllf(self):
-        R = 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
-
-    def plot(self, view='log'):
-        if self.is2D:
-            plot_result2D(self.data, self.theory(), view=view)
-        else:
-            plot_result1D(self.data, self.theory(), view=view)
-
-    def save(self, basename):
-        pass
-
-    def _get_weights(self, par):
-        from . import weights
-
-        relative = self.fn.info['partype']['pd-rel']
-        limits = self.fn.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, int(npts.value), width.value, nsigma.value,
-            value.value, limits[par], par in relative)
-        return v,w/w.max()
-
-
-def demo():
-    data = load_data('DEC07086.DAT')
-    set_beam_stop(data, 0.004)
-    plot_data(data)
-    import matplotlib.pyplot as plt; plt.show()
-
-
-if __name__ == "__main__":
-    demo()
Index: sasmodels/dll.py
===================================================================
--- sasmodels/dll.py	(revision ce27e213b076822e24ff35e97da5c4dd53eb45b0)
+++ sasmodels/dll.py	(revision ff7119b4ad066090cdffae80706883cef3d80399)
@@ -26,4 +26,6 @@
     for single and 'd', 'float64' or 'double' for double.  Double precision
     is an optional extension which may not be available on all devices.
+
+    Call :meth:`release` when done with the kernel.
     """
     def __init__(self, dllpath, info):
@@ -69,4 +71,7 @@
         return DllInput(q_vectors)
 
+    def release(self):
+        pass # TODO: should release the dll
+
 
 class DllInput(object):
@@ -100,8 +105,25 @@
 
 class DllKernel(object):
+    """
+    Callable SAS kernel.
+
+    *kernel* is the DllKernel object to call.
+
+    *info* is the module information
+
+    *input* is the DllInput q vectors at which the kernel should be
+    evaluated.
+
+    The resulting call method takes the *pars*, a list of values for
+    the fixed parameters to the kernel, and *pd_pars*, a list of (value,weight)
+    vectors for the polydisperse parameters.  *cutoff* determines the
+    integration limits: any points with combined weight less than *cutoff*
+    will not be calculated.
+
+    Call :meth:`release` when done with the kernel instance.
+    """
     def __init__(self, kernel, info, input):
         self.input = input
         self.kernel = kernel
-        self.info = info
         self.res = np.empty(input.nq, input.dtype)
         dim = '2d' if input.is_2D else '1d'
Index: sasmodels/gen.py
===================================================================
--- sasmodels/gen.py	(revision a7684e5745716d713f91b1ebc1027632a3ff9d36)
+++ sasmodels/gen.py	(revision ff7119b4ad066090cdffae80706883cef3d80399)
@@ -55,4 +55,35 @@
 them with the desired polydispersity.
 
+The available kernel parameters are defined as a list, with each parameter
+defined as a sublist with the following elements:
+
+    *name* is the name that will be used in the call to the kernel
+    function and the name that will be displayed to the user.  Names
+    should be lower case, with words separated by underscore.  If
+    acronyms are used, the whole acronym should be upper case.
+
+    *units* should be one of *degrees* for angles, *Ang* for lengths,
+    *1e-6/Ang^2* for SLDs.
+
+    *default value* will be the initial value for  the model when it
+    is selected, or when an initial value is not otherwise specified.
+
+    [*lb*, *ub*] are the hard limits on the parameter value, used to limit
+    the polydispersity density function.  In the fit, the parameter limits
+    given to the fit are the limits  on the central value of the parameter.
+    If there is polydispersity, it will evaluate parameter values outside
+    the fit limits, but not outside the hard limits specified in the model.
+    If there are no limits, use +/-inf imported from numpy.
+
+    *type* indicates how the parameter will be used.  "volume" parameters
+    will be used in all functions.  "orientation" parameters will be used
+    in *Iqxy* and *Imagnetic*.  "magnetic* parameters will be used in
+    *Imagnetic* only.  If *type* is the empty string, the parameter will
+    be used in all of *Iq*, *Iqxy* and *Imagnetic*.
+
+    *description* is a short description of the parameter.  This will
+    be displayed in the parameter table and used as a tool tip for the
+    parameter value in the user interface.
+
 The kernel module must set variables defining the kernel meta data:
 
@@ -65,7 +96,10 @@
     while the model parameters are being edited.
 
-    *parameters* is a list of parameters.  Each parameter is itself a
-    list containing *name*, *units*, *default value*,
-    [*lower bound*, *upper bound*], *type* and *description*.
+    *parameters* is the list of parameters.  Parameters in the kernel
+    functions must appear in the same order as they appear in the
+    parameters list.  Two additional parameters, *scale* and *background*
+    are added to the beginning of the parameter list.  They will show up
+    in the documentation as model parameters, but they are never sent to
+    the kernel functions.
 
     *source* is the list of C-99 source files that must be joined to
@@ -78,4 +112,29 @@
     *VR* is a python function defining the volume ratio.  If it is not
     present, the volume ratio is 1.
+
+An *info* dictionary is constructed from the kernel meta data and
+returned to the caller.  It includes the additional fields
+
+
+The model evaluator, function call sequence consists of q inputs and the return vector,
+followed by the loop value/weight vector, followed by the values for
+the non-polydisperse parameters, followed by the lengths of the
+polydispersity loops.  To construct the call for 1D models, the
+categories *fixed-1d* and *pd-1d* list the names of the parameters
+of the non-polydisperse and the polydisperse parameters respectively.
+Similarly, *fixed-2d* and *pd-2d* provide parameter names for 2D models.
+The *pd-rel* category is a set of those parameters which give
+polydispersitiy as a portion of the value (so a 10% length dispersity
+would use a polydispersity value of 0.1) rather than absolute
+dispersity such as an angle plus or minus 15 degrees.
+
+The *volume* category lists the volume parameters in order for calls
+to volume within the kernel (used for volume normalization) and for
+calls to ER and VR for effective radius and volume ratio respectively.
+
+The *orientation* and *magnetic* categories list the orientation and
+magnetic parameters.  These are used by the sasview interface.  The
+blank category is for parameters such as scale which don't have any
+other marking.
 
 The doc string at the start of the kernel module will be used to
@@ -86,43 +145,14 @@
 file names and extensions.
 
-Parameters are defined as follows:
-
-    *name* is the name that will be used in the call to the kernel
-    function and the name that will be displayed to the user.  Names
-    should be lower case, with words separated by underscore.  If
-    acronyms are used, the whole acronym should be upper case.
-
-    *units* should be one of *degrees* for angles, *Ang* for lengths,
-    *1e-6/Ang^2* for SLDs.
-
-    *default value* will be the initial value for  the model when it
-    is selected, or when an initial value is not otherwise specified.
-
-    *limits* are the valid bounds of the parameter, used to limit the
-    polydispersity density function.   In the fit, the parameter limits
-    given to the fit are the limits  on the central value of the parameter.
-    If there is polydispersity, it will evaluate parameter values outside
-    the fit limits, but not outside the hard limits specified in the model.
-    If there are no limits, use +/-inf imported from numpy.
-
-    *type* indicates how the parameter will be used.  "volume" parameters
-    will be used in all functions.  "orientation" parameters will be used
-    in *Iqxy* and *Imagnetic*.  "magnetic* parameters will be used in
-    *Imagnetic* only.  If *type* is none, the parameter will be used in
-    all of *Iq*, *Iqxy* and *Imagnetic*.  This will probably be a
-    is the empty string if the parameter is used in all model calculations,
-    it is "volu
-
-    *description* is a short description of the parameter.  This will
-    be displayed in the parameter table and used as a tool tip for the
-    parameter value in the user interface.
 
 The function :func:`make` loads the metadata from the module and returns
-the kernel source.  The function :func:`doc` extracts
+the kernel source.  The function :func:`doc` extracts the doc string
+and adds the parameter table to the top.  The function :func:`sources`
+returns a list of files required by the model.
 """
 
 # TODO: identify model files which have changed since loading and reload them.
 
-__all__ = ["make, doc"]
+__all__ = ["make, doc", "sources"]
 
 import os.path
@@ -158,8 +188,12 @@
     ]
 
+# Minimum width for a default value (this is shorter than the column header
+# width, so will be ignored).
 PARTABLE_VALUE_WIDTH = 10
 
 # Header included before every kernel.
-# This makes sure that the appropriate math constants are defined, and
+# This makes sure that the appropriate math constants are defined, and does
+# whatever is required to make the kernel compile as pure C rather than
+# as an OpenCL kernel.
 KERNEL_HEADER = """\
 // GENERATED CODE --- DO NOT EDIT ---
@@ -231,5 +265,5 @@
     }
 
-# Generic kernel template for opencl/openmp.
+# Generic kernel template for the polydispersity loop.
 # This defines the opencl kernel that is available to the host.  The same
 # structure is used for Iq and Iqxy kernels, so extra flexibility is needed
@@ -333,4 +367,7 @@
 
 def kernel_name(info, is_2D):
+    """
+    Name of the exported kernel symbol.
+    """
     return info['name'] + "_" + ("Iqxy" if is_2D else "Iq")
 
@@ -431,4 +468,7 @@
 
 def make_partable(info):
+    """
+    Generate the parameter table to include in the sphinx documentation.
+    """
     pars = info['parameters']
     column_widths = [
@@ -467,8 +507,21 @@
     raise ValueError("%r not found in %s"%(filename, search_path))
 
-def make_model(search_path, info):
+def sources(info):
+    """
+    Return a list of the sources file paths for the module.
+    """
+    from os.path import abspath, dirname, join as joinpath
+    search_path = [ dirname(info['filename']),
+                    abspath(joinpath(dirname(__file__),'models')) ]
+    return [_search(search_path) for f in info['source']]
+
+def make_model(info):
+    """
+    Generate the code for the kernel defined by info, using source files
+    found in the given search path.
+    """
     kernel_Iq = make_kernel(info, is_2D=False)
     kernel_Iqxy = make_kernel(info, is_2D=True)
-    source = [open(_search(search_path, f)).read() for f in info['source']]
+    source = [open(f).read() for f in sources(info)]
     kernel = "\n\n".join([KERNEL_HEADER]+source+[kernel_Iq, kernel_Iqxy])
     return kernel
@@ -479,25 +532,4 @@
 
     Returns a dictionary of categories.
-
-    The function call sequence consists of q inputs and the return vector,
-    followed by the loop value/weight vector, followed by the values for
-    the non-polydisperse parameters, followed by the lengths of the
-    polydispersity loops.  To construct the call for 1D models, the
-    categories *fixed-1d* and *pd-1d* list the names of the parameters
-    of the non-polydisperse and the polydisperse parameters respectively.
-    Similarly, *fixed-2d* and *pd-2d* provide parameter names for 2D models.
-    The *pd-rel* category is a set of those parameters which give
-    polydispersitiy as a portion of the value (so a 10% length dispersity
-    would use a polydispersity value of 0.1) rather than absolute
-    dispersity such as an angle plus or minus 15 degrees.
-
-    The *volume* category lists the volume parameters in order for calls
-    to volume within the kernel (used for volume normalization) and for
-    calls to ER and VR for effective radius and volume ratio respectively.
-
-    The *orientation* and *magnetic* categories list the orientation and
-    magnetic parameters.  These are used by the sasview interface.  The
-    blank category is for parameters such as scale which don't have any
-    other marking.
     """
     partype = {
@@ -535,5 +567,5 @@
     """
     # TODO: allow Iq and Iqxy to be defined in python
-    from os.path import abspath, dirname, join as joinpath
+    from os.path import abspath
     #print kernelfile
     info = dict(
@@ -550,7 +582,5 @@
     info['partype'] = categorize_parameters(info['parameters'])
 
-    search_path = [ dirname(info['filename']),
-                    abspath(joinpath(dirname(__file__),'models')) ]
-    source = make_model(search_path, info)
+    source = make_model(info)
 
     return source, info
@@ -565,4 +595,66 @@
                  doc=kernel_module.__doc__)
     return DOC_HEADER%subst
+
+# Compiler platform details
+if sys.platform == 'darwin':
+    COMPILE = "gcc-mp-4.7 -shared -fPIC -std=c99 -fopenmp -O2 -Wall %s -o %s -lm -lgomp"
+elif os.name == 'nt':
+    COMPILE = "gcc -shared -fPIC -std=c99 -fopenmp -O2 -Wall %s -o %s -lm"
+else:
+    COMPILE = "cc -shared -fPIC -std=c99 -fopenmp -O2 -Wall %s -o %s -lm"
+DLL_PATH = "/tmp"
+
+
+def dll_path(info):
+    """
+    Path to the compiled model defined by *info*.
+    """
+    from os.path import join as joinpath, split as splitpath, splitext
+    basename = splitext(splitpath(info['filename'])[1])[0]
+    return joinpath(DLL_PATH, basename+'.so')
+
+
+def dll_model(kernel_module, dtype=None):
+    """
+    Load the compiled model defined by *kernel_module*.
+
+    Recompile if any files are newer than the model file.
+
+    *dtype* is ignored.  Compiled files are always double.
+
+    The DLL is not loaded until the kernel is called so models an
+    be defined without using too many resources.
+    """
+    import tempfile
+    from sasmodels import dll
+
+    source, info = make(kernel_module)
+    source_files = sources(info) + [info['filename']]
+    newest = max(os.path.getmtime(f) for f in source_files)
+    dllpath = dll_path(info)
+    if not os.path.exists(dllpath) or os.path.getmtime(dllpath)<newest:
+        # Replace with a proper temp file
+        srcfile = tempfile.mkstemp(suffix=".c",prefix="sas_"+info['name'])
+        open(srcfile, 'w').write(source)
+        os.system(COMPILE%(srcfile, dllpath))
+        ## comment the following to keep the generated c file
+        os.unlink(srcfile)
+    return dll.DllModel(dllpath, info)
+
+
+def opencl_model(kernel_module, dtype="single"):
+    """
+    Load the OpenCL model defined by *kernel_module*.
+
+    Access to the OpenCL device is delayed until the kernel is called
+    so models can be defined without using too many resources.
+    """
+    from sasmodels import gpu
+
+    source, info = make(kernel_module)
+    ## for debugging, save source to a .cl file, edit it, and reload as model
+    #open(modelname+'.cl','w').write(source)
+    #source = open(modelname+'.cl','r').read()
+    return gpu.GpuModel(source, info, dtype)
 
 
Index: sasmodels/gpu.py
===================================================================
--- sasmodels/gpu.py	(revision 1780d5917f12726485d1b9b4d0343ac8454b6e34)
+++ sasmodels/gpu.py	(revision ff7119b4ad066090cdffae80706883cef3d80399)
@@ -242,4 +242,22 @@
 
 class GpuKernel(object):
+    """
+    Callable SAS kernel.
+
+    *kernel* is the GpuKernel object to call.
+
+    *info* is the module information
+
+    *input* is the DllInput q vectors at which the kernel should be
+    evaluated.
+
+    The resulting call method takes the *pars*, a list of values for
+    the fixed parameters to the kernel, and *pd_pars*, a list of (value,weight)
+    vectors for the polydisperse parameters.  *cutoff* determines the
+    integration limits: any points with combined weight less than *cutoff*
+    will not be calculated.
+
+    Call :meth:`release` when done with the kernel instance.
+    """
     def __init__(self, kernel, info, input):
         self.input = input
Index: smodels/jsonutil.py
===================================================================
--- sasmodels/jsonutil.py	(revision 14de349c663ca9521843af50415c37aed27a9eef)
+++ 	(revision )
@@ -1,172 +1,0 @@
-"""
-Read a relaxed JSON file.
-
-Relaxed JSON allows comments (introduced by // and going to the end
-of the line), optional quotes on key names in dictionaries and optional
-trailing commas at the ends of lists and dictionaries.  It also strips
-the leading characters up to the first '{'.  Multiline strings can be
-formatted using "\n\" at the end of each line.
-
-If the file contains e.g., "var _ = {...}", then it can be edited with
-a JavaScript aware editor such as VJET for eclipse, and it will be easier
-to locate errors and adjust formatting.
-"""
-import re
-import json
-from contextlib import contextmanager
-
-try:
-    from collections import OrderedDict
-except:
-    from ordered_dict import OrderedDict
-
-
-_LEADING_TEXT = re.compile(r'^.*?[{]', re.DOTALL)
-_LINE_CONTINUATION = re.compile(r'\\\s*\n')
-_TRAILING_COMMENT = re.compile(r'(?P<comment>\s*//.*?)\n')
-_MULTILINE_COMMENT = re.compile(r'(?P<comment>/\*.*?\*/)', re.DOTALL)
-_UNQUOTED_FIELDNAME = re.compile(r'(?P<prefix>[,{]\s*)(?P<key>[^\s,{}:"]+)(?P<tail>\s*:)')
-_TRAILING_COMMA = re.compile(r',(?P<tail>\s*[]}])')
-
-def relaxed_load(path, **kw):
-    return relaxed_loads(open(path).read(), **kw)
-
-def relaxed_loads(text, **kw):
-    """
-    Parse and return a relaxed JSON string.
-    """
-    ordered = kw.pop('ordered', False)
-    if ordered:  kw['object_pairs_hook'] = OrderedDict
-    # TODO: need a little state machine that performs the translation so that
-    # TODO: line and column numbers are preserved, and so that we can have
-    # TODO: http:// in a string (instead of it being treated like a comment). 
-    #print "== raw text\n", text
-    text = _LINE_CONTINUATION.sub('', text)
-    #print "== joined lines\n", text
-    text = _TRAILING_COMMENT.sub(r'\n', text)
-    text = _MULTILINE_COMMENT.sub(r'', text)
-    #print "== stripped comments\n", text
-    text = _LEADING_TEXT.sub('{', text)
-    #print "== trimmed text\n", text
-    text = _UNQUOTED_FIELDNAME.sub(r'\g<prefix>"\g<key>"\g<tail>', text)
-    #print "== quoted field names\n", text
-    text = _TRAILING_COMMA.sub(r'\g<tail>', text)
-    #print "== processed text\n", text
-    try:
-        obj = json.loads(text, object_hook=decode_dict_as_str, **kw)
-    except ValueError, e:
-        msg = [str(e)]
-        M = re.findall('line ([0-9]*) column ([0-9]*)', msg[0])
-        if M:
-            line,col = int(M[0][0]), int(M[0][1])
-            lines = text.split("\n")
-            if line>=2: msg.append(lines[line-2])
-            if line>=1: msg.append(lines[line-1])
-            msg.append(" "*(col-1) + "^")
-            if line<len(lines): msg.append(lines[line])
-        msg = "\n".join(msg)
-        raise e.__class__(msg)
-    return obj
-
-@contextmanager
-def float_format(formatstr='.15g'):
-    """
-    Allow the float format to be changed for a json encoding action.
-
-    This is a context manager, and should be used for example as::
-
-        >>> with float_format('.2g'):
-        >>>    print json.dumps(sqrt(2))
-        1.41 
-    """
-    formatter = json.encoder.FLOAT_REPR
-    json.encoder.FLOAT_REPR = lambda o: format(o, formatstr)
-    yield
-    json.encoder.FLOAT_REPR = formatter
-
-def numpy_encoder(o):
-    """
-    JSON encoder for numpy data.
-
-    To automatically convert numpy data to lists when writing a datastream
-    use json.dumps(object, default=numpy_json).
-    """
-    try:
-        return o.tolist()
-    except AttributeError:
-        raise TypeError
-
-
-def _decode_list(lst):
-    newlist = []
-    for i in lst:
-        if isinstance(i, unicode):
-            i = i.encode('utf-8')
-        elif isinstance(i, list):
-            i = _decode_list(i)
-        newlist.append(i)
-    return newlist
-
-def decode_dict_as_str(dct):
-    newdict = {}
-    for k, v in dct.iteritems():
-        if isinstance(k, unicode):
-            k = k.encode('utf-8')
-        if isinstance(v, unicode):
-            v = v.encode('utf-8')
-        elif isinstance(v, list):
-            v = _decode_list(v)
-        newdict[k] = v
-    return newdict
-
-
-def test():
-    """
-    Verify that the translation from pseudo-JSON to JSON works.
-    """
-    good = """\
-// This is a source definition with no errors
-var entry = {
-field : { // A comment about the field
-  "field" : "te\\
-x\\  
-t",
-  other$field : 56,
-  },
-/*
-multiline comment
-*/
-secondfield : {
-  content: ["string", "string"], /* a second comment */
-  content: [{name:"good", value:3, URL:"http:\\/\\/my.url.com"},]
-  },
-}
-"""
-    broken = """\
-// This is a source definition with a missing comma
-{
-field : { // A comment about the field
-  field : "te\\   
-x\\  
-t"
-  other$field : 56,
-  },
-/*
-multiline comment
-*/
-secondfield : {
-  content: ["string", "string"], /* a second comment */
-  content: [{name:"good", value:3},]
-  },
-}
-"""
-    result = relaxed_loads(good)
-    assert result['field']['field'] == "text"
-    assert result['field']['other$field'] == 56
-    assert result['secondfield']['content'][0]['name'] == 'good'
-    try: relaxed_loads(broken)
-    except ValueError, _: pass
-    else: raise Exception("No exception raised in broken")
-
-if __name__ == "__main__":
-    test()
Index: sasmodels/models/README.rst
===================================================================
--- sasmodels/models/README.rst	(revision a7684e5745716d713f91b1ebc1027632a3ff9d36)
+++ sasmodels/models/README.rst	(revision ff7119b4ad066090cdffae80706883cef3d80399)
@@ -33,2 +33,5 @@
 will mean extending the data handler to support multiple cross sections
 in the same data set.
+
+Need to write code to generate the polydispersity loops in python for
+kernels that are only implemented in python.
Index: sasmodels/models/cylinder.c
===================================================================
--- sasmodels/models/cylinder.c	(revision a7684e5745716d713f91b1ebc1027632a3ff9d36)
+++ sasmodels/models/cylinder.c	(revision ff7119b4ad066090cdffae80706883cef3d80399)
@@ -14,12 +14,15 @@
     real length)
 {
-    const real h = REAL(0.5)*length;
+    const real halflength = REAL(0.5)*length;
     real summ = REAL(0.0);
+    // real lower=0, upper=M_PI_2;
     for (int i=0; i<76 ;i++) {
-        //const real zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
+        // translate a point in [-1,1] to a point in [lower,upper]
+        //const real zi = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0;
         const real zi = REAL(0.5)*(Gauss76Z[i]*M_PI_2 + M_PI_2);
-        summ += Gauss76Wt[i] * CylKernel(q, radius, h, zi);
+        summ += Gauss76Wt[i] * CylKernel(q, radius, halflength, zi);
     }
-    //const real form = (uplim-lolim)/2.0*summ;
+    // translate dx in [-1,1] to dx in [lower,upper]
+    //const real form = (upper-lower)/2.0*summ;
     const real form = summ * M_PI_4;
 
Index: sasmodels/models/cylinder_clone.c
===================================================================
--- sasmodels/models/cylinder_clone.c	(revision 32c160aab45a38358aa878cf4622c2cec2d87725)
+++ sasmodels/models/cylinder_clone.c	(revision ff7119b4ad066090cdffae80706883cef3d80399)
@@ -45,5 +45,5 @@
     // # The following correction factor exists in sasview, but it can't be
     // # right, so we are leaving it out for now.
-    // const real correction = fabs(cn)*M_PI_2;
+    const real spherical_integration = fabs(cn)*M_PI_2;
     const real q = sqrt(qx*qx+qy*qy);
     const real cos_val = cn*cos(cyl_phi*M_PI_180)*(qx/q) + sn*(qy/q);
@@ -65,4 +65,4 @@
     // The additional volume factor is for polydisperse volume normalization.
     const real s = (sldCyl - sldSolv) * form_volume(radius, length);
-    return REAL(1.0e8) * form * s * s; // * correction;
+    return REAL(1.0e8) * form * s * s * spherical_integration;
 }
Index: sasmodels/models/lib/cylkernel.c
===================================================================
--- sasmodels/models/lib/cylkernel.c	(revision 14de349c663ca9521843af50415c37aed27a9eef)
+++ sasmodels/models/lib/cylkernel.c	(revision ff7119b4ad066090cdffae80706883cef3d80399)
@@ -1,16 +1,19 @@
 // Compute the form factor for a cylinder
-// qq is the q-value for the calculation (1/A)
+//     F^2(q) sin a
+// given
+//     F(q) = sin(Q L/2 cos a)/(Q L/2 cos a) * 2 J1(Q R sin a) / (Q R sin a)
+// q is the q-value for the calculation (1/A)
 // radius is the radius of the cylinder (A)
-// h is the HALF-LENGTH of the cylinder = L/2 (A)
-real CylKernel(real q, real radius, real h, real theta);
-real CylKernel(real q, real radius, real h, real theta)
+// halflength is the HALF-LENGTH of the cylinder = L/2 (A)
+real CylKernel(real q, real radius, real halflength, real theta);
+real CylKernel(real q, real radius, real halflength, real theta)
 {
     real sn, cn;
     SINCOS(theta, sn, cn);
     const real besarg = q*radius*sn;
-    const real siarg = q*h*cn;
+    const real siarg = q*halflength*cn;
     // lim_{x->0} J1(x)/x = 1/2,  lim_{x->0} sin(x)/x = 1
     const real bj = (besarg == REAL(0.0) ? REAL(0.5) : J1(besarg)/besarg);
     const real si = (siarg == REAL(0.0) ? REAL(1.0) : sin(siarg)/siarg);
-    return REAL(4.0)*sin(theta)*bj*bj*si*si;
+    return REAL(4.0)*sn*bj*bj*si*si;
 }
Index: sasmodels/sasview_model.py
===================================================================
--- sasmodels/sasview_model.py	(revision a7684e5745716d713f91b1ebc1027632a3ff9d36)
+++ sasmodels/sasview_model.py	(revision ff7119b4ad066090cdffae80706883cef3d80399)
@@ -1,10 +1,24 @@
 import math
 from copy import deepcopy
+import warnings
 
 import numpy as np
 
+try:
+    import pyopencl
+    from .gen import opencl_model as load_model
+except ImportError:
+    warnings.warn("OpenCL not available --- using ctypes instead")
+    from .gen import dll_model as load_model
+
+
 def make_class(kernel_module, dtype='single'):
-    from .core import opencl_model
-    model =  opencl_model(kernel_module, dtype=dtype)
+    """
+    Load the sasview model defined in *kernel_module*.
+    :param kernel_module:
+    :param dtype:
+    :return:
+    """
+    model =  load_model(kernel_module, dtype=dtype)
     def __init__(self, multfactor=1):
         SasviewModel.__init__(self, model)
@@ -242,4 +256,12 @@
 
     def calculate_Iq(self, *args):
+        """
+        Calculate Iq for one set of q with the current parameters.
+
+        If the model is 1D, use *q*.  If 2D, use *qx*, *qy*.
+
+        This should NOT be used for fitting since it copies the *q* vectors
+        to the card for each evaluation.
+        """
         q_vectors = [np.asarray(q) for q in args]
         fn = self._model(self._model.make_input(q_vectors))
Index: sasmodels/weights.py
===================================================================
--- sasmodels/weights.py	(revision 1780d5917f12726485d1b9b4d0343ac8454b6e34)
+++ sasmodels/weights.py	(revision ff7119b4ad066090cdffae80706883cef3d80399)
@@ -1,2 +1,5 @@
+"""
+SAS distributions for polydispersity.
+"""
 # TODO: include dispersion docs with the disperser models
 from math import sqrt
@@ -116,4 +119,5 @@
 
 
+# dispersion name -> disperser lookup table.
 models = dict((d.type,d) for d in (
     GaussianDispersion, RectangleDispersion,
@@ -123,4 +127,24 @@
 
 def get_weights(disperser, n, width, nsigmas, value, limits, relative):
+    """
+    Return the set of values and weights for a polydisperse parameter.
+
+    *disperser* is the name of the disperser.
+
+    *n* is the number of points in the weight vector.
+
+    *width* is the width of the disperser distribution.
+
+    *nsigmas* is the number of sigmas to span for the dispersion convolution.
+
+    *value* is the value of the parameter in the model.
+
+    *limits* is [lb, ub], the lower and upper bound of the weight value.
+
+    *relative* is true if *width* is defined in proportion to the value
+    of the parameter, and false if it is an absolute width.
+
+    Returns *(value,weight)*, where *value* and *weight* are vectors.
+    """
     cls = models[disperser]
     obj = cls(n, width, nsigmas)
