Index: sasmodels/dll.py
===================================================================
--- sasmodels/dll.py	(revision cb6ecf4c19ecd47d43632a230eb29c56fba61a7c)
+++ sasmodels/dll.py	(revision b3f6bc3b9123ab59a035e8360b67cc7d380ea695)
@@ -11,4 +11,5 @@
 
 from . import generate
+from .pykernel import PyInput, PyKernel
 
 from .generate import F32, F64
@@ -116,6 +117,11 @@
 
     def __call__(self, input):
+        # Support pure python kernel call
+        if input.is_2D and callable(self.info['Iqxy']):
+            return PyKernel(self.info['Iqxy'], self.info, input)
+        elif not input.is_2D and callable(self.info['Iq']):
+            return PyKernel(self.info['Iq'], self.info, 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)
@@ -125,10 +131,9 @@
         Make q input vectors available to the model.
 
-        This only needs to be done once for all models that operate on the
-        same input.  So for example, if you are adding two different models
-        together to compare to a data set, then only one model needs to
-        needs to call make_input, so long as the models have the same dtype.
+        Note that each model needs its own q vector even if the case of
+        mixture models because some models may be OpenCL, some may be
+        ctypes and some may be pure python.
         """
-        return DllInput(q_vectors)
+        return PyInput(q_vectors, dtype=F64)
 
     def release(self):
@@ -136,38 +141,9 @@
 
 
-class DllInput(object):
-    """
-    Make q data available to the gpu.
-
-    *q_vectors* is a list of q vectors, which will be *[q]* for 1-D data,
-    and *[qx, qy]* for 2-D data.  Internally, the vectors will be reallocated
-    to get the best performance on OpenCL, which may involve shifting and
-    stretching the array to better match the memory architecture.  Additional
-    points will be evaluated with *q=1e-3*.
-
-    *dtype* is the data type for the q vectors. The data type should be
-    set to match that of the kernel, which is an attribute of
-    :class:`GpuProgram`.  Note that not all kernels support double
-    precision, so even if the program was created for double precision,
-    the *GpuProgram.dtype* may be single precision.
-
-    Call :meth:`release` when complete.  Even if not called directly, the
-    buffer will be released when the data object is freed.
-    """
-    def __init__(self, q_vectors):
-        self.nq = q_vectors[0].size
-        self.dtype = np.dtype('double')
-        self.is_2D = (len(q_vectors) == 2)
-        self.q_vectors = [np.ascontiguousarray(q,self.dtype) for q in q_vectors]
-        self.q_pointers = [q.ctypes.data for q in q_vectors]
-
-    def release(self):
-        self.q_vectors = []
-
 class DllKernel(object):
     """
     Callable SAS kernel.
 
-    *kernel* is the DllKernel object to call.
+    *kernel* is the c function to call.
 
     *info* is the module information
Index: sasmodels/generate.py
===================================================================
--- sasmodels/generate.py	(revision cb6ecf4c19ecd47d43632a230eb29c56fba61a7c)
+++ sasmodels/generate.py	(revision b3f6bc3b9123ab59a035e8360b67cc7d380ea695)
@@ -354,4 +354,17 @@
 """
 
+
+
+##########################################################
+#                                                        #
+#   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
+#   !!                                              !!   #
+#   !!  KEEP THIS CODE CONSISTENT WITH PYKERNEL.PY  !!   #
+#   !!                                              !!   #
+#   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
+#                                                        #
+##########################################################
+
+
 # Polydispersity loop body.
 # This computes the weight, and if it is sufficient, calls the scattering
@@ -387,6 +400,6 @@
 # to call the form_volume function from the user supplied kernel, and accumulate
 # a normalized weight.
-VOLUME_NORM="""const double vol_weight = %(weight)s;
-    vol += vol_weight*form_volume(%(pars)s);
+VOLUME_NORM="""const double vol_weight = %(vol_weight)s;
+    vol += vol_weight*form_volume(%(vol_pars)s);
     norm_vol += vol_weight;\
 """
@@ -490,6 +503,6 @@
     if vol_pars:
         subst = {
-            'weight': "*".join(p+"_w" for p in vol_pars),
-            'pars': ", ".join(vol_pars),
+            'vol_weight': "*".join(p+"_w" for p in vol_pars),
+            'vol_pars': ", ".join(vol_pars),
             }
         volume_norm = VOLUME_NORM%subst
@@ -637,6 +650,6 @@
             }
         source.append(WORK_FUNCTION%subst)
-    kernel_Iq = make_kernel(info, is_2D=False)
-    kernel_Iqxy = make_kernel(info, is_2D=True)
+    kernel_Iq = make_kernel(info, is_2D=False) if not callable(info['Iq']) else ""
+    kernel_Iqxy = make_kernel(info, is_2D=True) if not callable(info['Iqxy']) else ""
     kernel = "\n\n".join([KERNEL_HEADER]+source+[kernel_Iq, kernel_Iqxy])
     return kernel
@@ -718,15 +731,15 @@
 def demo_time():
     import datetime
+    from .models import cylinder
+    toc = lambda: (datetime.datetime.now()-tic).total_seconds()
     tic = datetime.datetime.now()
-    toc = lambda: (datetime.datetime.now()-tic).total_seconds()
-    path = os.path.dirname("__file__")
-    doc, c = make_model(os.path.join(path, "models", "cylinder.c"))
+    source, info = make(cylinder)
     print "time:",toc()
 
 def demo():
-    from os.path import join as joinpath, dirname
-    c, info, doc = make_model(joinpath(dirname(__file__), "models", "cylinder.c"))
-    #print doc
-    #print c
+    from .models import cylinder
+    source, info = make(cylinder)
+    #print doc(cylinder)
+    print source
 
 if __name__ == "__main__":
Index: sasmodels/gpu.py
===================================================================
--- sasmodels/gpu.py	(revision cb6ecf4c19ecd47d43632a230eb29c56fba61a7c)
+++ sasmodels/gpu.py	(revision b3f6bc3b9123ab59a035e8360b67cc7d380ea695)
@@ -24,8 +24,10 @@
 """
 import numpy as np
+
 import pyopencl as cl
 from pyopencl import mem_flags as mf
 
 from . import generate
+from .pykernel import PyInput, PyKernel
 
 F64_DEFS = """\
@@ -185,4 +187,10 @@
 
     def __call__(self, input):
+        # Support pure python kernel call
+        if input.is_2D and callable(self.info['Iqxy']):
+            return PyKernel(self.info['Iqxy'], self.info, input)
+        elif not input.is_2D and callable(self.info['Iq']):
+            return PyKernel(self.info['Iq'], self.info, input)
+
         if self.dtype != input.dtype:
             raise TypeError("data and kernel have different types")
@@ -202,14 +210,15 @@
         Make q input vectors available to the model.
 
-        This only needs to be done once for all models that operate on the
-        same input.  So for example, if you are adding two different models
-        together to compare to a data set, then only one model needs to
-        needs to call make_input, so long as the models have the same dtype.
+        Note that each model needs its own q vector even if the case of
+        mixture models because some models may be OpenCL, some may be
+        ctypes and some may be pure python.
         """
-        # Note: the weird interface, where make_input doesn't care which
-        # model calls it, allows us to ask the model to define the data
-        # and the caller does not need to know if it is opencl or ctypes.
-        # The returned data object is opaque.
-        return GpuInput(q_vectors, dtype=self.dtype)
+        # Support pure python kernel call
+        if len(q_vectors) == 1 and callable(self.info['Iq']):
+            return PyInput(q_vectors, dtype=self.dtype)
+        elif callable(self.info['Iqxy']):
+            return PyInput(q_vectors, dtype=self.dtype)
+        else:
+            return GpuInput(q_vectors, dtype=self.dtype)
 
 # TODO: check that we don't need a destructor for buffers which go out of scope
Index: sasmodels/models/spherepy.py
===================================================================
--- sasmodels/models/spherepy.py	(revision be802cba5f7b2aa210e548392cb72633a33d2c7b)
+++ sasmodels/models/spherepy.py	(revision b3f6bc3b9123ab59a035e8360b67cc7d380ea695)
@@ -83,19 +83,25 @@
 
 
-# No volume normalization despite having a volume parameter
-# This should perhaps be volume normalized?
 def form_volume(radius):
-    return 1.333333333333333*pi*radius*radius*radius
+    return 1.333333333333333*pi*radius**3
 
 def Iq(q, sld, solvent_sld, radius):
+    #print "q",q
+    #print "sld,r",sld,solvent_sld,radius
     qr = q*radius
     sn, cn = sin(qr), cos(qr)
-    bes = 3 * (sn-qr*cn)/qr**3
+    # FOR VECTORIZED VERSION, UNCOMMENT THE NEXT TWO LINES
+    bes = 3 * (sn-qr*cn)/qr**3 # may be 0/0 but we fix that next line
     bes[qr==0] = 1
+    # FOR NON VECTORIZED VERSION, UNCOMMENT THE NEXT LINE
+    #bes = 3 * (sn-qr*cn)/qr**3 if qr>0 else 1
     fq = bes * (sld - solvent_sld) * form_volume(radius)
-    return 1.0e-4*fq*fq
+    return 1.0e-4*fq**2
+# FOR VECTORIZED VERSION, UNCOMMENT THE NEXT LINE
+Iq.vectorized = True
 
 def Iqxy(qx, qy, sld, solvent_sld, radius):
     return Iq(sqrt(qx**2 + qy**2), sld, solvent_sld, radius)
+Iqxy.vectorized = True
 
 def ER(radius):
Index: sasmodels/pykernel.py
===================================================================
--- sasmodels/pykernel.py	(revision b3f6bc3b9123ab59a035e8360b67cc7d380ea695)
+++ sasmodels/pykernel.py	(revision b3f6bc3b9123ab59a035e8360b67cc7d380ea695)
@@ -0,0 +1,209 @@
+import numpy as np
+from numpy import pi, sin, cos, sqrt
+
+from .generate import F32, F64
+
+class PyInput(object):
+    """
+    Make q data available to the gpu.
+
+    *q_vectors* is a list of q vectors, which will be *[q]* for 1-D data,
+    and *[qx, qy]* for 2-D data.  Internally, the vectors will be reallocated
+    to get the best performance on OpenCL, which may involve shifting and
+    stretching the array to better match the memory architecture.  Additional
+    points will be evaluated with *q=1e-3*.
+
+    *dtype* is the data type for the q vectors. The data type should be
+    set to match that of the kernel, which is an attribute of
+    :class:`GpuProgram`.  Note that not all kernels support double
+    precision, so even if the program was created for double precision,
+    the *GpuProgram.dtype* may be single precision.
+
+    Call :meth:`release` when complete.  Even if not called directly, the
+    buffer will be released when the data object is freed.
+    """
+    def __init__(self, q_vectors, dtype):
+        self.nq = q_vectors[0].size
+        self.dtype = dtype
+        self.is_2D = (len(q_vectors) == 2)
+        self.q_vectors = [np.ascontiguousarray(q,self.dtype) for q in q_vectors]
+        self.q_pointers = [q.ctypes.data for q in q_vectors]
+
+    def release(self):
+        self.q_vectors = []
+
+class PyKernel(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.info = info
+        self.input = input
+        self.res = np.empty(input.nq, input.dtype)
+        dim = '2d' if input.is_2D else '1d'
+        # Loop over q unless user promises that the kernel is vectorized by
+        # taggining it with vectorized=True
+        if not getattr(kernel, 'vectorized', False):
+            if dim == '2d':
+                def vector_kernel(qx, qy, *args):
+                    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])
+            self.kernel = vector_kernel
+        else:
+            self.kernel = kernel
+        fixed_pars = info['partype']['fixed-'+dim]
+        pd_pars = info['partype']['pd-'+dim]
+        vol_pars = info['partype']['volume']
+
+        # 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])
+        try: self.theta_index = pars.index('theta')+offset
+        except ValueError: self.theta_index = -1
+
+        # Caller needs fixed_pars and pd_pars
+        self.fixed_pars = fixed_pars
+        self.pd_pars = pd_pars
+
+    def __call__(self, fixed, pd, cutoff):
+        #print "fixed",fixed
+        #print "pd", pd
+        args = self.args[:]  # grab a copy of the args
+        form, form_volume = self.kernel, self.info['form_volume']
+        # First two fixed
+        scale, background = fixed[:2]
+        for index,value in zip(self.fixed_index, fixed[2:]):
+            args[index] = float(value)
+        res = _loops(form, form_volume, cutoff, scale, background,  args,
+                     pd, self.pd_index, self.vol_index, self.theta_index)
+
+        return res
+
+    def release(self):
+        self.input = None
+
+def _loops(form, form_volume, cutoff, scale, background,
+           args, pd, pd_index, vol_index, theta_index):
+    """
+    *form* is the name of the form function, which should be vectorized over
+    q, but otherwise have an interface like the opencl kernels, with the
+    q parameters first followed by the individual parameters in the order
+    given in model.parameters (see :mod:`sasmodels.generate`).
+
+    *form_volume* calculates the volume of the shape.  *vol_index* gives
+    the list of volume parameters
+
+    *cutoff* ignores the corners of the dispersion hypercube
+
+    *scale*, *background* multiplies the resulting form and adds an offset
+
+    *args* is the prepopulated set of arguments to the form function, starting
+    with the q vectors, and including slots for all the parameters.  The
+    values for the parameters will be substituted with values from the
+    dispersion functions.
+
+    *pd* is the list of dispersion parameters
+
+    *pd_index* are the indices of the dispersion parameters in *args*
+
+    *vol_index* are the indices of the volume parameters in *args*
+
+    *theta_index* is the index of the theta parameter for the sasview
+    spherical correction, or -1 if there is no angular dispersion
+    """
+
+    ##########################################################
+    #                                                        #
+    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
+    #   !!                                              !!   #
+    #   !!  KEEP THIS CODE CONSISTENT WITH GENERATE.PY  !!   #
+    #   !!                                              !!   #
+    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
+    #                                                        #
+    ##########################################################
+
+    # weight vector, to be populated by polydispersity loops
+    weight = np.empty(len(pd), 'd')
+
+    # identify which pd parameters are volume parameters
+    vol_weight_index = np.array([(index in vol_index) for index in pd_index])
+
+    # Sort parameters in decreasing order of pd length
+    Npd = np.array([len(pdi[0]) for pdi in pd], 'i')
+    order = np.argsort(Npd)[::-1]
+    stride = np.cumprod(Npd[order])
+    pd = [pd[index] for index in order]
+    pd_index = pd_index[order]
+    vol_weight_index = vol_weight_index[order]
+
+    fast_value = pd[0][0]
+    fast_weight = pd[0][1]
+
+    ret = np.zeros_like(args[0])
+    norm = np.zeros_like(ret)
+    vol = np.zeros_like(ret)
+    vol_norm = np.zeros_like(ret)
+    for k in range(stride[-1]):
+        # update polydispersity parameter values
+        fast_index = k%stride[0]
+        if fast_index == 0:  # bottom loop complete ... check all other loops
+            for i,index, in enumerate(k%stride):
+                args[pd_index[i]] = pd[i][0][index]
+                weight[i] = pd[i][1][index]
+        else:
+            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.
+        # Note: make sure this is consistent with the code in PY_LOOP_BODY!!
+        # Note: can precompute w1*w2*...*wn
+        w = np.prod(weight)
+        if w > cutoff:
+            I = form(*args)
+            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
+            spherical_correction = 1.0
+            ret += w*I*spherical_correction*positive
+            norm += w*positive
+
+            # 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 form_volume:
+                vol_args = [args[index] for index in vol_index]
+                vol_weight = np.prod(weight[vol_weight_index])
+                vol += vol_weight*form_volume(*vol_args)*positive
+                vol_norm += vol_weight*positive
+
+    positive = (vol*vol_norm != 0.0)
+    ret[positive] *= vol_norm[positive] / vol[positive]
+    result = scale*ret/norm+background
+    return result
