Index: doc/developer/calculator.rst
===================================================================
--- doc/developer/calculator.rst	(revision 6fd8de01432581dc71f1290cf5070acac4860fd7)
+++ doc/developer/calculator.rst	(revision 6fd8de01432581dc71f1290cf5070acac4860fd7)
@@ -0,0 +1,225 @@
+Calculator Interface
+====================
+
+The environment needs to provide the following #defines:
+
+- USE_OPENCL is defined if running in opencl
+- KERNEL declares a function to be available externally
+- KERNEL_NAME is the name of the function being declared
+- MAX_PD is the maximum depth of the polydispersity loop [model specific]
+- NPARS is the number of parameters in the kernel
+- PARAMETER_TABLE is the declaration of the parameters to the kernel::
+
+    Cylinder:
+
+        #define PARAMETER_TABLE \
+        double length; \
+        double radius; \
+        double sld; \
+        double sld_solvent
+
+    Note: scale and background are never included
+
+    Multi-shell cylinder (10 shell max):
+
+        #define PARAMETER_TABLE \
+        double num_shells; \
+        double length; \
+        double radius[10]; \
+        double sld[10]; \
+        double sld_solvent
+
+- CALL_IQ(q, i, var) is the declaration of a call to the kernel::
+
+    Cylinder:
+
+        #define CALL_IQ(q, i, var) Iq(q[i], \
+        var.length, \
+        var.radius, \
+        var.sld, \
+        var.sld_solvent)
+
+    Multi-shell cylinder:
+
+        #define CALL_IQ(q, i, var) Iq(q[i], \
+        var.num_shells, \
+        var.length, \
+        var.radius, \
+        var.sld, \
+        var.sld_solvent)
+
+    Cylinder2D:
+
+        #define CALL_IQ(q, i, var) Iqxy(q[2*i], q[2*i+1], \
+        var.length, \
+        var.radius, \
+        var.sld, \
+        var.sld_solvent, \
+        var.theta, \
+        var.phi)
+
+- CALL_VOLUME(var) is similar, but for calling the form volume::
+
+        #define CALL_VOLUME(var) \
+        form_volume(var.length, var.radius)
+
+- INVALID(var) is a test for model parameters in the correct range::
+
+    Cylinder:
+
+        #define INVALID(var) 0
+
+    BarBell:
+
+        #define INVALID(var) (var.bell_radius < var.radius)
+
+    Model with complicated constraints:
+
+        inline bool constrained(p1, p2, p3) { return expression; }
+        #define INVALID(var) constrained(var.p1, var.p2, var.p3)
+
+Our design supports a limited number of polydispersity loops, wherein
+we need to cycle through the values of the polydispersity, calculate
+the I(q, p) for each combination of parameters, and perform a normalized
+weighted sum across all the weights.  Parameters may be passed to the
+underlying calculation engine as scalars or vectors, but the polydispersity
+calculator treats the parameter set as one long vector.
+
+Let's assume we have 6 parameters in the model, with two polydisperse::
+
+    0: scale        {scl = constant}
+    1: background   {bkg = constant}
+    2: length       {l = vector of 30pts}
+    3: radius       {r = vector of 10pts}
+    4: sld          {s = constant/(radius**2*length)}
+    5: sld_solvent  {s2 = constant}
+
+This generates the following call to the kernel (where x stands for an
+arbitrary value that is not used by the kernel evaluator)::
+
+    NPARS = 4  // scale and background are in all models
+    problem {
+        pd_par = {3, 2, x, x}         // parameters *radius* and *length* vary
+        pd_length = {30, 10, 0, 0}    // *length* has more, so it is first
+        pd_offset = {10, 0, x, x}     // *length* starts at index 10 in weights
+        pd_stride = {1, 30, 300, 300} // cumulative product of pd length
+        pd_isvol = {True, True, x, x}       // true if weight is a volume weight
+        par_offset = {2, 3, 303, 313}  // parameter offsets
+        par_coord = {0, 3, 2, 1} // bitmap of parameter dependencies
+        fast_coord_index = {3, 5, x, x} // radius and sld have fast index
+        fast_coord_count = 2  // two parameters vary with *length* distribution
+        theta_var = -1   // no spherical correction
+        fast_theta = 0   // spherical correction angle is not pd 1
+    }
+
+    weight = { l0, .., l29, r0, .., r9} //length comes first as the longest vec
+    pars = { scl, bkg, l0, ..., l29, r0, r1, ..., r9,
+             s[l0,r0], ... s[l0,r9], s[l1,r0], ... s[l29,r9] , s2}
+             //where s[x,y] stands for material sld, s2 = solvent sld
+
+    nq = 130
+    q = { q0, q1, ..., q130, x, x }  # pad to 8 element boundary
+    result = {r1, ..., r130, norm, vol, vol_norm, x, x, x, x, x, x, x}
+
+
+The polydisperse parameters are stored in as an array of parameter
+indices, one for each polydisperse parameter, stored in pd_par[n].
+Non-polydisperse parameters do not appear in this array. Each polydisperse
+parameter has a weight vector whose length is stored in pd_length[n].
+The weights are stored in a contiguous vector of weights for all
+parameters, with the starting position for the each parameter stored
+in pd_offset[n].  The values corresponding to the weights are stored
+together in a separate weights[] vector, with offset stored in
+par_offset[pd_par[n]]. Polydisperse parameters should be stored in
+decreasing order of length for highest efficiency.
+
+We limit the number of polydisperse dimensions to MAX_PD (currently 4),
+though some models may have fewer if they have fewer polydisperse
+parameters. This cuts the size of the structure in half compared to
+allowing a separate polydispersity for each parameter.  This will
+help a little bit for models with large numbers of parameters, such
+as the onion model.
+
+Parameters may be coordinated.  That is, we may have the value of one
+parameter depend on a set of other parameters, some of which may be
+polydisperse.  For example, if sld is inversely proportional to the
+volume of a cylinder, and the length and radius are independently
+polydisperse, then for each combination of length and radius we need a
+separate value for the sld.  The caller must provide a coordination table
+for each parameter containing the value for each parameter given the
+value of the polydisperse parameters v1, v2, etc.  The tables for each
+parameter are arranged contiguously in a vector, with offset[k] giving the
+starting location of parameter k in the vector.  Each parameter defines
+par_coord[k] as a bit mask indicating which polydispersity parameters the
+parameter depends upon. Usually this is zero, indicating that the parameter
+is independent, but for the cylinder example given, the bits for the
+radius and length polydispersity parameters would both be set, the result
+being a (#radius x #length) table, or maybe a (#length x #radius) table
+if length comes first in the polydispersity table.
+
+NB: If we can guarantee that a compiler and OpenCL driver are available,
+we could instead create the coordination function on the fly for each
+parameter, saving memory and transfer time, but requiring a C compiler
+as part of the environment.
+
+In ordering the polydisperse parameters by decreasing length we can
+iterate over the longest dispersion weight vector first.  All parameters
+coordinated with this weight vector (the 'fast' parameters), can be
+updated with a simple increment to the next position in the parameter
+value table.  The indices of these parameters is stored in fast_coord_index[],
+with fast_coord_count being the number of fast parameters.  A total
+of NPARS slots is allocated to allow for the case that all parameters
+are coordinated with the fast index, though this will likely be mostly
+empty.  When the fast increment count reaches the end of the weight
+vector, then the index of the second polydisperse parameter must be
+incremented, and all of its coordinated parameters updated.  Because this
+operation is not in the inner loop, a slower algorithm can be used.
+
+If there is no polydispersity we pretend that it is polydisperisty with one
+parameter, pd_start=0 and pd_stop=1.  We may or may not short circuit the
+calculation in this case, depending on how much time it saves.
+
+The problem details structure can be allocated and sent in as an integer
+array using the read-only flag.  This allows us to copy it once per fit
+along with the weights vector, since features such as the number of
+polydisperity elements per pd parameter or the coordinated won't change
+between function evaluations.  A new parameter vector is sent for
+each I(q) evaluation.
+
+To protect against expensive evaluations taking all the GPU resource
+on large fits, the entire polydispersity will not be computed at once.
+Instead, a start and stop location will be sent, indicating where in the
+polydispersity loop the calculation should start and where it should
+stop.  We can do this for arbitrary start/stop points since we have
+unwound the nested loop.  Instead, we use the same technique as array
+index translation, using div and mod to figure out the i,j,k,...
+indices in the virtual nested loop.
+
+The results array will be initialized to zero for polydispersity loop
+entry zero, and preserved between calls to [start, stop] so that the
+results accumulate by the time the loop has completed.  Background and
+scale will be applied when the loop reaches the end.  This does require
+that the results array be allocated read-write, which is less efficient
+for the GPU, but it makes the calling sequence much more manageable.
+
+Scale and background cannot be coordinated with other polydisperse parameters
+
+Oriented objects in 2-D need a spherical correction on the angular variation
+in order to preserve the 'surface area' of the weight distribution.
+
+cutoff parameter limits integration area within polydispersity hypercude,
+which speeds calculations
+
+For accuracy we may want to introduce Kahan summation into the integration::
+
+
+    double accumulated_error = 0.0;
+    ...
+    #if USE_KAHAN_SUMMATION
+        const double y = next - accumulated_error;
+        const double t = ret + y;
+        accumulated_error = (t - ret) - y;
+        ret = t;
+    #else
+        ret += next;
+    #endif
Index: doc/developer/index.rst
===================================================================
--- doc/developer/index.rst	(revision 56fc97aa819cfe686101e6c6b3bac4c2c7cfa01e)
+++ doc/developer/index.rst	(revision b85be2da3b04a7ce8ef87e2a85a44572b2b1480c)
@@ -3,3 +3,7 @@
 ***************************
 
+.. toctree::
+   :numbered: 4
+   :maxdepth: 4
 
+   calculator.rst
Index: sasmodels/bumps_model.py
===================================================================
--- sasmodels/bumps_model.py	(revision 9217ef8174f5566df0d9259319dc1bde48372242)
+++ sasmodels/bumps_model.py	(revision 9217ef8174f5566df0d9259319dc1bde48372242)
@@ -81,9 +81,9 @@
     from bumps.names import Parameter
 
-    pars = {}
+    pars = {} # => floating point parameters
     for p in model_info['parameters']:
         value = kwargs.pop(p.name, p.default)
         pars[p.name] = Parameter.default(value, name=p.name, limits=p.limits)
-    for name in model_info['partype']['pd-2d']:
+    for name in model_info['par_type']['pd']:
         for xpart, xdefault, xlimits in [
                 ('_pd', 0., pars[name].limits),
@@ -95,6 +95,6 @@
             pars[xname] = Parameter.default(xvalue, name=xname, limits=xlimits)
 
-    pd_types = {}
-    for name in model_info['partype']['pd-2d']:
+    pd_types = {}  # => distribution names
+    for name in model_info['par_type']['pd']:
         xname = name + '_pd_type'
         xvalue = kwargs.pop(xname, 'gaussian')
Index: sasmodels/compare.py
===================================================================
--- sasmodels/compare.py	(revision 1671636e9c5f4ddecde92b003b0a89e86eaef0c8)
+++ sasmodels/compare.py	(revision 69aa451cd7078d77c255be72d487793b00b34eee)
@@ -309,6 +309,5 @@
         exclude = lambda n: False
     else:
-        partype = model_info['partype']
-        par1d = set(partype['fixed-1d']+partype['pd-1d'])
+        par1d = model_info['par_type']['1d']
         exclude = lambda n: n not in par1d
     lines = []
@@ -454,5 +453,6 @@
     """
     # initialize the code so time is more accurate
-    value = calculator(**suppress_pd(pars))
+    if Nevals > 1:
+        value = calculator(**suppress_pd(pars))
     toc = tic()
     for _ in range(max(Nevals, 1)):  # make sure there is at least one eval
@@ -877,11 +877,10 @@
         pars, pd_types = bumps_model.create_parameters(model_info, **opts['pars'])
         if not opts['is2d']:
-            active = [base + ext
-                      for base in model_info['partype']['pd-1d']
-                      for ext in ['', '_pd', '_pd_n', '_pd_nsigma']]
-            active.extend(model_info['partype']['fixed-1d'])
-            for k in active:
-                v = pars[k]
-                v.range(*parameter_range(k, v.value))
+            for p in model_info['parameters'].type['1d']:
+                for ext in ['', '_pd', '_pd_n', '_pd_nsigma']:
+                    k = p.name+ext
+                    v = pars.get(k, None)
+                    if v is not None:
+                        v.range(*parameter_range(k, v.value))
         else:
             for k, v in pars.items():
Index: sasmodels/core.py
===================================================================
--- sasmodels/core.py	(revision e79f0a543061a4d9cac5c2ffb7c6e6366e733570)
+++ sasmodels/core.py	(revision 69aa451cd7078d77c255be72d487793b00b34eee)
@@ -26,5 +26,5 @@
 __all__ = [
     "list_models", "load_model_info", "precompile_dll",
-    "build_model", "make_kernel", "call_kernel", "call_ER_VR",
+    "build_model", "call_kernel", "call_ER_VR",
 ]
 
@@ -167,11 +167,5 @@
 
 
-def make_kernel(model, q_vectors):
-    """
-    Return a computation kernel from the model definition and the q input.
-    """
-    return model(q_vectors)
-
-def get_weights(model_info, pars, name):
+def get_weights(parameter, values):
     """
     Generate the distribution for parameter *name* given the parameter values
@@ -181,11 +175,13 @@
     from the *pars* dictionary for parameter value and parameter dispersion.
     """
-    relative = name in model_info['partype']['pd-rel']
-    limits = model_info['limits'][name]
-    disperser = pars.get(name+'_pd_type', 'gaussian')
-    value = pars.get(name, model_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)
+    value = values.get(parameter.name, parameter.default)
+    if parameter.type not in ('volume', 'orientation'):
+        return np.array([value]), np.array([1.0])
+    relative = parameter.type == 'volume'
+    limits = parameter.limits
+    disperser = values.get(parameter.name+'_pd_type', 'gaussian')
+    npts = values.get(parameter.name+'_pd_n', 0)
+    width = values.get(parameter.name+'_pd', 0.0)
+    nsigma = values.get(parameter.name+'_pd_nsigma', 3.0)
     value, weight = weights.get_weights(
         disperser, npts, width, nsigma, value, limits, relative)
@@ -208,5 +204,5 @@
 def call_kernel(kernel, pars, cutoff=0, mono=False):
     """
-    Call *kernel* returned from :func:`make_kernel` with parameters *pars*.
+    Call *kernel* returned from *model.make_kernel* with parameters *pars*.
 
     *cutoff* is the limiting value for the product of dispersion weights used
@@ -216,13 +212,32 @@
     with an error of about 1%, which is usually less than the measurement
     uncertainty.
-    """
-    fixed_pars = [pars.get(name, kernel.info['defaults'][name])
-                  for name in kernel.fixed_pars]
+
+    *mono* is True if polydispersity should be set to none on all parameters.
+    """
     if mono:
-        pd_pars = [( np.array([pars[name]]), np.array([1.0]) )
-                   for name in kernel.pd_pars]
-    else:
-        pd_pars = [get_weights(kernel.info, pars, name) for name in kernel.pd_pars]
-    return kernel(fixed_pars, pd_pars, cutoff=cutoff)
+        active = lambda name: False
+    elif kernel.dim == '1d':
+        pars_1d = set(p.name for p in kernel.info['parameters'].type['1d'])
+        active = lambda name: name in pars_1d
+    elif kernel.dim == '2d':
+        pars_2d = set(p.name for p in kernel.info['parameters'].type['2d'])
+        active = lambda name: name in pars_1d
+    else:
+        active = lambda name: True
+
+    vw_pairs = [(get_weights(p, pars) if active(p.name) else ([p.default], [1]))
+                for p in kernel.info['parameters']]
+    values, weights = zip(*vw_pairs)
+
+    if max([len(w) for w in weights]) > 1:
+        details = generate.poly_details(kernel.info, weights)
+    else:
+        details = kernel.info['mono_details']
+
+    weights, values = [np.hstack(v) for v in (weights, values)]
+
+    weights = weights.astype(dtype=kernel.dtype)
+    values = values.astype(dtype=kernel.dtype)
+    return kernel(details, weights, values, cutoff)
 
 def call_ER_VR(model_info, vol_pars):
@@ -247,16 +262,17 @@
 
 
-def call_ER(info, pars):
-    """
-    Call the model ER function using *pars*.
-    *info* is either *model.info* if you have a loaded model, or *kernel.info*
-    if you have a model kernel prepared for evaluation.
-    """
-    ER = info.get('ER', None)
+def call_ER(model_info, values):
+    """
+    Call the model ER function using *values*. *model_info* is either
+    *model.info* if you have a loaded model, or *kernel.info* if you
+    have a model kernel prepared for evaluation.
+    """
+    ER = model_info.get('ER', None)
     if ER is None:
         return 1.0
     else:
-        vol_pars = [get_weights(info, pars, name)
-                    for name in info['partype']['volume']]
+        vol_pars = [get_weights(parameter, values)
+                    for parameter in model_info['parameters']
+                    if parameter.type == 'volume']
         value, weight = dispersion_mesh(vol_pars)
         individual_radii = ER(*value)
@@ -264,5 +280,5 @@
         return np.sum(weight*individual_radii) / np.sum(weight)
 
-def call_VR(info, pars):
+def call_VR(model_info, values):
     """
     Call the model VR function using *pars*.
@@ -270,10 +286,11 @@
     if you have a model kernel prepared for evaluation.
     """
-    VR = info.get('VR', None)
+    VR = model_info.get('VR', None)
     if VR is None:
         return 1.0
     else:
-        vol_pars = [get_weights(info, pars, name)
-                    for name in info['partype']['volume']]
+        vol_pars = [get_weights(parameter, values)
+                    for parameter in model_info['parameters']
+                    if parameter.type == 'volume']
         value, weight = dispersion_mesh(vol_pars)
         whole, part = VR(*value)
Index: sasmodels/direct_model.py
===================================================================
--- sasmodels/direct_model.py	(revision 02e70ffbd2f37e9a91d5d027972b30069bd05c6b)
+++ sasmodels/direct_model.py	(revision 48fbd50a25db17a9c6bf7770aeb061514b3e8879)
@@ -25,5 +25,4 @@
 import numpy as np
 
-from .core import make_kernel
 from .core import call_kernel, call_ER_VR
 from . import sesans
@@ -66,6 +65,4 @@
             self.data_type = 'Iq'
 
-        partype = model.info['partype']
-
         if self.data_type == 'sesans':
             
@@ -81,4 +78,5 @@
             q_mono = sesans.make_all_q(data)
         elif self.data_type == 'Iqxy':
+            partype = model.info['par_type']
             if not partype['orientation'] and not partype['magnetic']:
                 raise ValueError("not 2D without orientation or magnetic parameters")
@@ -154,9 +152,13 @@
     def _calc_theory(self, pars, cutoff=0.0):
         if self._kernel is None:
-            self._kernel = make_kernel(self._model, self._kernel_inputs)  # pylint: disable=attribute-dedata_type
-            self._kernel_mono = make_kernel(self._model, self._kernel_mono_inputs) if self._kernel_mono_inputs else None
+            self._kernel = self._model.make_kernel(self._kernel_inputs)
+            self._kernel_mono = (
+                self._model.make_kernel(self._kernel_mono_inputs)
+                if self._kernel_mono_inputs else None
+            )
 
         Iq_calc = call_kernel(self._kernel, pars, cutoff=cutoff)
-        Iq_mono = call_kernel(self._kernel_mono, pars, mono=True) if self._kernel_mono_inputs else None
+        Iq_mono = (call_kernel(self._kernel_mono, pars, mono=True)
+                   if self._kernel_mono_inputs else None)
         if self.data_type == 'sesans':
             result = sesans.transform(self._data,
Index: sasmodels/generate.py
===================================================================
--- sasmodels/generate.py	(revision 9ef9dd9f11c54d995dd94b650a4563230b868352)
+++ sasmodels/generate.py	(revision 69aa451cd7078d77c255be72d487793b00b34eee)
@@ -21,4 +21,8 @@
 
     *VR(p1, p2, ...)* returns the volume ratio for core-shell style forms.
+
+    #define INVALID(v) (expr)  returns False if v.parameter is invalid
+    for some parameter or other (e.g., v.bell_radius < v.radius).  If
+    necessary, the expression can call a function.
 
 These functions are defined in a kernel module .py script and an associated
@@ -65,39 +69,4 @@
 The constructor code will generate the necessary vectors for computing
 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.
-
-    *limits = [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*.  "sld" parameters
-    can automatically be promoted to magnetic parameters, each of which
-    will have a magnitude and a direction, which may be different from
-    other sld parameters.
-
-    *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:
 
@@ -212,23 +181,23 @@
 from __future__ import print_function
 
-# TODO: identify model files which have changed since loading and reload them.
+#TODO: determine which functions are useful outside of generate
+#__all__ = ["model_info", "make_doc", "make_source", "convert_type"]
 
 import sys
 from os.path import abspath, dirname, join as joinpath, exists, basename, \
-    splitext
+    splitext, getmtime
 import re
 import string
 import warnings
-from collections import namedtuple
 
 import numpy as np
 
-PARAMETER_FIELDS = ['name', 'units', 'default', 'limits', 'type', 'description']
-Parameter = namedtuple('Parameter', PARAMETER_FIELDS)
-
-#TODO: determine which functions are useful outside of generate
-#__all__ = ["model_info", "make_doc", "make_source", "convert_type"]
-
-C_KERNEL_TEMPLATE_PATH = joinpath(dirname(__file__), 'kernel_template.c')
+from .modelinfo import ModelInfo, Parameter, make_parameter_table
+
+# TODO: identify model files which have changed since loading and reload them.
+
+TEMPLATE_ROOT = dirname(__file__)
+
+MAX_PD = 4
 
 F16 = np.dtype('float16')
@@ -239,10 +208,4 @@
 except TypeError:
     F128 = None
-
-# Scale and background, which are parameters common to every form factor
-COMMON_PARAMETERS = [
-    ["scale", "", 1, [0, np.inf], "", "Source intensity"],
-    ["background", "1/cm", 1e-3, [0, np.inf], "", "Source background"],
-    ]
 
 # Conversion from units defined in the parameter table for each model
@@ -338,4 +301,5 @@
     raise ValueError("%r not found in %s" % (filename, search_path))
 
+
 def model_sources(model_info):
     """
@@ -346,17 +310,14 @@
     return [_search(search_path, f) for f in model_info['source']]
 
-# Pragmas for enable OpenCL features.  Be sure to protect them so that they
-# still compile even if OpenCL is not present.
-_F16_PRAGMA = """\
-#if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp16)
-#  pragma OPENCL EXTENSION cl_khr_fp16: enable
-#endif
-"""
-
-_F64_PRAGMA = """\
-#if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp64)
-#  pragma OPENCL EXTENSION cl_khr_fp64: enable
-#endif
-"""
+def timestamp(model_info):
+    """
+    Return a timestamp for the model corresponding to the most recently
+    changed file or dependency.
+    """
+    source_files = (model_sources(model_info)
+                    + model_templates()
+                    + [model_info['filename']])
+    newest = max(getmtime(f) for f in source_files)
+    return newest
 
 def convert_type(source, dtype):
@@ -369,5 +330,5 @@
     if dtype == F16:
         fbytes = 2
-        source = _F16_PRAGMA + _convert_type(source, "half", "f")
+        source = _convert_type(source, "float", "f")
     elif dtype == F32:
         fbytes = 4
@@ -375,5 +336,5 @@
     elif dtype == F64:
         fbytes = 8
-        source = _F64_PRAGMA + source  # Source is already double
+        # no need to convert if it is already double
     elif dtype == F128:
         fbytes = 16
@@ -418,28 +379,72 @@
 
 
-LOOP_OPEN = """\
-for (int %(name)s_i=0; %(name)s_i < N%(name)s; %(name)s_i++) {
-  const double %(name)s = loops[2*(%(name)s_i%(offset)s)];
-  const double %(name)s_w = loops[2*(%(name)s_i%(offset)s)+1];\
+_template_cache = {}
+def load_template(filename):
+    path = joinpath(TEMPLATE_ROOT, filename)
+    mtime = getmtime(path)
+    if filename not in _template_cache or mtime > _template_cache[filename][0]:
+        with open(path) as fid:
+            _template_cache[filename] = (mtime, fid.read(), path)
+    return _template_cache[filename][1]
+
+def model_templates():
+    # TODO: fails DRY; templates are listed in two places.
+    # should instead have model_info contain a list of paths
+    return [joinpath(TEMPLATE_ROOT, filename)
+            for filename in ('kernel_header.c', 'kernel_iq.c')]
+
+
+_FN_TEMPLATE = """\
+double %(name)s(%(pars)s);
+double %(name)s(%(pars)s) {
+    %(body)s
+}
+
+
 """
-def build_polydispersity_loops(pd_pars):
-    """
-    Build polydispersity loops
-
-    Returns loop opening and loop closing
-    """
-    depth = 4
-    offset = ""
-    loop_head = []
-    loop_end = []
-    for name in pd_pars:
-        subst = {'name': name, 'offset': offset}
-        loop_head.append(indent(LOOP_OPEN % subst, depth))
-        loop_end.insert(0, (" "*depth) + "}")
-        offset += '+N' + name
-        depth += 2
-    return "\n".join(loop_head), "\n".join(loop_end)
-
-C_KERNEL_TEMPLATE = None
+
+def _gen_fn(name, pars, body):
+    """
+    Generate a function given pars and body.
+
+    Returns the following string::
+
+         double fn(double a, double b, ...);
+         double fn(double a, double b, ...) {
+             ....
+         }
+    """
+    par_decl = ', '.join(p.as_argument() for p in pars) if pars else 'void'
+    return _FN_TEMPLATE % {'name': name, 'body': body, 'pars': par_decl}
+
+def _call_pars(prefix, pars):
+    """
+    Return a list of *prefix.parameter* from parameter items.
+    """
+    return [p.as_call_reference(prefix) for p in pars]
+
+_IQXY_PATTERN = re.compile("^((inline|static) )? *(double )? *Iqxy *([(]|$)",
+                           flags=re.MULTILINE)
+def _have_Iqxy(sources):
+    """
+    Return true if any file defines Iqxy.
+
+    Note this is not a C parser, and so can be easily confused by
+    non-standard syntax.  Also, it will incorrectly identify the following
+    as having Iqxy::
+
+        /*
+        double Iqxy(qx, qy, ...) { ... fill this in later ... }
+        */
+
+    If you want to comment out an Iqxy function, use // on the front of the
+    line instead.
+    """
+    for code in sources:
+        if _IQXY_PATTERN.search(code):
+            return True
+    else:
+        return False
+
 def make_source(model_info):
     """
@@ -460,174 +465,103 @@
     # for computing volume even if we allow non-disperse volume parameters.
 
-    # Load template
-    global C_KERNEL_TEMPLATE
-    if C_KERNEL_TEMPLATE is None:
-        with open(C_KERNEL_TEMPLATE_PATH) as fid:
-            C_KERNEL_TEMPLATE = fid.read()
-
-    # Load additional sources
-    source = [open(f).read() for f in model_sources(model_info)]
-
-    # Prepare defines
-    defines = []
-    partype = model_info['partype']
-    pd_1d = partype['pd-1d']
-    pd_2d = partype['pd-2d']
-    fixed_1d = partype['fixed-1d']
-    fixed_2d = partype['fixed-1d']
-
-    iq_parameters = [p.name
-                     for p in model_info['parameters'][2:]  # skip scale, background
-                     if p.name in set(fixed_1d + pd_1d)]
-    iqxy_parameters = [p.name
-                       for p in model_info['parameters'][2:]  # skip scale, background
-                       if p.name in set(fixed_2d + pd_2d)]
-    volume_parameters = [p.name
-                         for p in model_info['parameters']
-                         if p.type == 'volume']
-
-    # Fill in defintions for volume parameters
-    if volume_parameters:
-        defines.append(('VOLUME_PARAMETERS',
-                        ','.join(volume_parameters)))
-        defines.append(('VOLUME_WEIGHT_PRODUCT',
-                        '*'.join(p + '_w' for p in volume_parameters)))
-
-    # Generate form_volume function from body only
+    partable = model_info['parameters']
+
+    # Identify parameters for Iq, Iqxy, Iq_magnetic and form_volume.
+    # Note that scale and volume are not possible types.
+
+    # Load templates and user code
+    kernel_header = load_template('kernel_header.c')
+    kernel_code = load_template('kernel_iq.c')
+    user_code = [open(f).read() for f in model_sources(model_info)]
+
+    # Build initial sources
+    source = [kernel_header] + user_code
+
+    vol_parameters = partable.kernel_pars('volume')
+    iq_parameters = partable.kernel_pars('1d')
+    iqxy_parameters = partable.kernel_pars('2d')
+
+    # Make parameters for q, qx, qy so that we can use them in declarations
+    q, qx, qy = [Parameter(name=v) for v in ('q', 'qx', 'qy')]
+    # Generate form_volume function, etc. from body only
     if model_info['form_volume'] is not None:
-        if volume_parameters:
-            vol_par_decl = ', '.join('double ' + p for p in volume_parameters)
-        else:
-            vol_par_decl = 'void'
-        defines.append(('VOLUME_PARAMETER_DECLARATIONS',
-                        vol_par_decl))
-        fn = """\
-double form_volume(VOLUME_PARAMETER_DECLARATIONS);
-double form_volume(VOLUME_PARAMETER_DECLARATIONS) {
-    %(body)s
-}
-""" % {'body':model_info['form_volume']}
-        source.append(fn)
-
-    # Fill in definitions for Iq parameters
-    defines.append(('IQ_KERNEL_NAME', model_info['name'] + '_Iq'))
-    defines.append(('IQ_PARAMETERS', ', '.join(iq_parameters)))
-    if fixed_1d:
-        defines.append(('IQ_FIXED_PARAMETER_DECLARATIONS',
-                        ', \\\n    '.join('const double %s' % p for p in fixed_1d)))
-    if pd_1d:
-        defines.append(('IQ_WEIGHT_PRODUCT',
-                        '*'.join(p + '_w' for p in pd_1d)))
-        defines.append(('IQ_DISPERSION_LENGTH_DECLARATIONS',
-                        ', \\\n    '.join('const int N%s' % p for p in pd_1d)))
-        defines.append(('IQ_DISPERSION_LENGTH_SUM',
-                        '+'.join('N' + p for p in pd_1d)))
-        open_loops, close_loops = build_polydispersity_loops(pd_1d)
-        defines.append(('IQ_OPEN_LOOPS',
-                        open_loops.replace('\n', ' \\\n')))
-        defines.append(('IQ_CLOSE_LOOPS',
-                        close_loops.replace('\n', ' \\\n')))
+        pars = vol_parameters
+        source.append(_gen_fn('form_volume', pars, model_info['form_volume']))
     if model_info['Iq'] is not None:
-        defines.append(('IQ_PARAMETER_DECLARATIONS',
-                        ', '.join('double ' + p for p in iq_parameters)))
-        fn = """\
-double Iq(double q, IQ_PARAMETER_DECLARATIONS);
-double Iq(double q, IQ_PARAMETER_DECLARATIONS) {
-    %(body)s
-}
-""" % {'body':model_info['Iq']}
-        source.append(fn)
-
-    # Fill in definitions for Iqxy parameters
-    defines.append(('IQXY_KERNEL_NAME', model_info['name'] + '_Iqxy'))
-    defines.append(('IQXY_PARAMETERS', ', '.join(iqxy_parameters)))
-    if fixed_2d:
-        defines.append(('IQXY_FIXED_PARAMETER_DECLARATIONS',
-                        ', \\\n    '.join('const double %s' % p for p in fixed_2d)))
-    if pd_2d:
-        defines.append(('IQXY_WEIGHT_PRODUCT',
-                        '*'.join(p + '_w' for p in pd_2d)))
-        defines.append(('IQXY_DISPERSION_LENGTH_DECLARATIONS',
-                        ', \\\n    '.join('const int N%s' % p for p in pd_2d)))
-        defines.append(('IQXY_DISPERSION_LENGTH_SUM',
-                        '+'.join('N' + p for p in pd_2d)))
-        open_loops, close_loops = build_polydispersity_loops(pd_2d)
-        defines.append(('IQXY_OPEN_LOOPS',
-                        open_loops.replace('\n', ' \\\n')))
-        defines.append(('IQXY_CLOSE_LOOPS',
-                        close_loops.replace('\n', ' \\\n')))
+        pars = [q] + iq_parameters
+        source.append(_gen_fn('Iq', pars, model_info['Iq']))
     if model_info['Iqxy'] is not None:
-        defines.append(('IQXY_PARAMETER_DECLARATIONS',
-                        ', '.join('double ' + p for p in iqxy_parameters)))
-        fn = """\
-double Iqxy(double qx, double qy, IQXY_PARAMETER_DECLARATIONS);
-double Iqxy(double qx, double qy, IQXY_PARAMETER_DECLARATIONS) {
-    %(body)s
-}
-""" % {'body':model_info['Iqxy']}
-        source.append(fn)
-
-    # Need to know if we have a theta parameter for Iqxy; it is not there
-    # for the magnetic sphere model, for example, which has a magnetic
-    # orientation but no shape orientation.
-    if 'theta' in pd_2d:
-        defines.append(('IQXY_HAS_THETA', '1'))
-
-    #for d in defines: print(d)
-    defines = '\n'.join('#define %s %s' % (k, v) for k, v in defines)
-    sources = '\n\n'.join(source)
-    return C_KERNEL_TEMPLATE % {
-        'DEFINES': defines,
-        'SOURCES': sources,
-        }
+        pars = [qx, qy] + iqxy_parameters
+        source.append(_gen_fn('Iqxy', pars, model_info['Iqxy']))
+
+    # Define the parameter table
+    source.append("#define PARAMETER_TABLE \\")
+    source.append("\\\n".join(p.as_definition()
+                                  for p in model_info['parameters'][2:]))
+
+    # Define the function calls
+    if vol_parameters:
+        refs = _call_pars("v.", vol_parameters)
+        call_volume = "#define CALL_VOLUME(v) form_volume(%s)" % (",".join(refs))
+    else:
+        # Model doesn't have volume.  We could make the kernel run a little
+        # faster by not using/transferring the volume normalizations, but
+        # the ifdef's reduce readability more than is worthwhile.
+        call_volume = "#define CALL_VOLUME(v) 0.0"
+    source.append(call_volume)
+
+    refs = ["q[i]"] + _call_pars("v.", iq_parameters)
+    call_iq = "#define CALL_IQ(q,i,v) Iq(%s)" % (",".join(refs))
+    if _have_Iqxy(user_code):
+        # Call 2D model
+        refs = ["q[2*i]", "q[2*i+1]"] + _call_pars("v.", iqxy_parameters)
+        call_iqxy = "#define CALL_IQ(q,i,v) Iqxy(%s)" % (",".join(refs))
+    else:
+        # Call 1D model with sqrt(qx^2 + qy^2)
+        warnings.warn("Creating Iqxy = Iq(sqrt(qx^2 + qy^2))")
+        # still defined:: refs = ["q[i]"] + _call_pars("v", iq_parameters)
+        pars_sqrt = ["sqrt(q[2*i]*q[2*i]+q[2*i+1]*q[2*i+1])"] + refs[1:]
+        call_iqxy = "#define CALL_IQ(q,i,v) Iq(%s)" % (",".join(pars_sqrt))
+
+    # Fill in definitions for numbers of parameters
+    source.append("#define MAX_PD %s"%model_info['max_pd'])
+    source.append("#define NPARS %d"%(len(partable.kernel_pars())))
+
+    # TODO: allow mixed python/opencl kernels?
+
+    # define the Iq kernel
+    source.append("#define KERNEL_NAME %s_Iq"%model_info['name'])
+    source.append(call_iq)
+    source.append(kernel_code)
+    source.append("#undef CALL_IQ")
+    source.append("#undef KERNEL_NAME")
+
+    # define the Iqxy kernel from the same source with different #defines
+    source.append("#define KERNEL_NAME %s_Iqxy"%model_info['name'])
+    source.append(call_iqxy)
+    source.append(kernel_code)
+    source.append("#undef CALL_IQ")
+    source.append("#undef KERNEL_NAME")
+
+    return '\n'.join(source)
 
 def categorize_parameters(pars):
     """
-    Build parameter categories out of the the parameter definitions.
-
-    Returns a dictionary of categories.
-
-    Note: these categories are subject to change, depending on the needs of
-    the UI and the needs of the kernel calling function.
-
-    The categories are as follows:
-
-    * *volume* list of volume parameter names
-    * *orientation* list of orientation parameters
-    * *magnetic* list of magnetic parameters
-    * *<empty string>* list of parameters that have no type info
-
-    Each parameter is in one and only one category.
-
-    The following derived categories are created:
-
-    * *fixed-1d* list of non-polydisperse parameters for 1D models
-    * *pd-1d* list of polydisperse parameters for 1D models
-    * *fixed-2d* list of non-polydisperse parameters for 2D models
-    * *pd-d2* list of polydisperse parameters for 2D models
-    """
-    partype = {
-        'volume': [], 'orientation': [], 'magnetic': [], 'sld': [], '': [],
-        'fixed-1d': [], 'fixed-2d': [], 'pd-1d': [], 'pd-2d': [],
-        'pd-rel': set(),
-    }
-
-    for p in pars:
-        if p.type == 'volume':
-            partype['pd-1d'].append(p.name)
-            partype['pd-2d'].append(p.name)
-            partype['pd-rel'].add(p.name)
-        elif p.type == 'magnetic':
-            partype['fixed-2d'].append(p.name)
-        elif p.type == 'orientation':
-            partype['pd-2d'].append(p.name)
-        elif p.type in ('', 'sld'):
-            partype['fixed-1d'].append(p.name)
-            partype['fixed-2d'].append(p.name)
-        else:
-            raise ValueError("unknown parameter type %r" % p.type)
-        partype[p.type].append(p.name)
-
-    return partype
+    Categorize the parameters by use:
+
+    * *pd* list of polydisperse parameters in order; gui should test whether
+      they are in *2d* or *magnetic* as appropriate for the data
+    * *1d* set of parameters that are used to compute 1D patterns
+    * *2d* set of parameters that are used to compute 2D patterns (which
+      includes all 1D parameters)
+    * *magnetic* set of parameters that are used to compute magnetic
+      patterns (which includes all 1D and 2D parameters)
+    * *pd_relative* is the set of parameters with relative distribution
+      width (e.g., radius +/- 10%) rather than absolute distribution
+      width (e.g., theta +/- 6 degrees).
+    * *theta_par* is the index of the polar angle polydispersion parameter
+      or -1 if no such parameter exists
+    """
+    par_set = {}
 
 def process_parameters(model_info):
@@ -635,21 +569,105 @@
     Process parameter block, precalculating parameter details.
     """
-    # convert parameters into named tuples
-    for p in model_info['parameters']:
-        if p[4] == '' and (p[0].startswith('sld') or p[0].endswith('sld')):
-            p[4] = 'sld'
-            # TODO: make sure all models explicitly label their sld parameters
-            #raise ValueError("%s.%s needs to be explicitly set to type 'sld'" %(model_info['id'], p[0]))
-
-    pars = [Parameter(*p) for p in model_info['parameters']]
-    # Fill in the derived attributes
-    model_info['parameters'] = pars
-    partype = categorize_parameters(pars)
-    model_info['limits'] = dict((p.name, p.limits) for p in pars)
-    model_info['partype'] = partype
-    model_info['defaults'] = dict((p.name, p.default) for p in pars)
+    partable = model_info['parameters']
     if model_info.get('demo', None) is None:
-        model_info['demo'] = model_info['defaults']
-    model_info['has_2d'] = partype['orientation'] or partype['magnetic']
+        model_info['demo'] = partable.defaults
+
+    # Don't use more polydisperse parameters than are available in the model
+    # Note: we can do polydispersity on arbitrary parameters, so it is not
+    # clear that this is a good idea; it does however make the poly_details
+    # code easier to write, so we will leave it in for now.
+    model_info['max_pd'] = min(partable.num_pd, MAX_PD)
+
+def mono_details(model_info):
+    # TODO: move max_pd into ParameterTable?
+    max_pd = model_info['max_pd']
+    pars = model_info['parameters'].kernel_pars()
+    npars = len(pars)
+    par_offset = 5*max_pd
+    constants_offset = par_offset + 3*npars
+
+    details = np.zeros(constants_offset + 2, 'int32')
+    details[0*max_pd:1*max_pd] = range(max_pd)       # pd_par: arbitrary order; use first
+    details[1*max_pd:2*max_pd] = [1]*max_pd          # pd_length: only one element
+    details[2*max_pd:3*max_pd] = range(max_pd)       # pd_offset: consecutive 1.0 weights
+    details[3*max_pd:4*max_pd] = [1]*max_pd          # pd_stride: vectors of length 1
+    details[4*max_pd:5*max_pd] = [0]*max_pd          # pd_isvol: doens't matter if no norm
+    details[par_offset+0*npars:par_offset+1*npars] = range(2, npars+2) # par_offset: skip scale and background
+    details[par_offset+1*npars:par_offset+2*npars] = [0]*npars         # no coordination
+    #details[p+npars] = 1 # par_coord[0] is coordinated with the first par?
+    details[par_offset+2*npars:par_offset+3*npars] = 0 # fast coord with 0
+    details[constants_offset]   = 1     # fast_coord_count: one fast index
+    details[constants_offset+1] = -1    # theta_par: None
+    return details
+
+def poly_details(model_info, weights):
+    weights = weights[2:]
+
+    # TODO: move max_pd into ParameterTable?
+    max_pd = model_info['max_pd']
+    pars = model_info['parameters'].kernel_pars
+    npars = len(pars)
+    par_offset = 5*max_pd
+    constants_offset = par_offset + 3*npars
+
+    # Decreasing list of polydispersity lengths
+    # Note: the reversing view, x[::-1], does not require a copy
+    pd_length = np.array([len(w) for w in weights])
+    print (pd_length)
+    print (weights)
+    pd_offset = np.cumsum(np.hstack((0, pd_length)))
+    pd_isvol = np.array([p.type=='volume' for p in pars])
+    idx = np.argsort(pd_length)[::-1][:max_pd]
+    print (idx)
+    pd_stride = np.cumprod(np.hstack((1, pd_length[idx][:-1])))
+    par_offsets = np.cumsum(np.hstack((2, pd_length)))[:-1]
+
+    theta_par = -1
+    if 'theta_par' in model_info:
+        theta_par = model_info['theta_par']
+        if theta_par >= 0 and pd_length[theta_par] <= 1:
+            theta_par = -1
+
+    details = np.empty(constants_offset + 2, 'int32')
+    details[0*max_pd:1*max_pd] = idx             # pd_par
+    details[1*max_pd:2*max_pd] = pd_length[idx]
+    details[2*max_pd:3*max_pd] = pd_offset[idx]
+    details[3*max_pd:4*max_pd] = pd_stride
+    details[4*max_pd:5*max_pd] = pd_isvol[idx]
+    details[par_offset+0*npars:par_offset+1*npars] = par_offsets
+    details[par_offset+1*npars:par_offset+2*npars] = 0  # no coordination for most
+    details[par_offset+2*npars:par_offset+3*npars] = 0  # no fast coord with 0
+    coord_offset = par_offset+1*npars
+    for k,parameter_num in enumerate(idx):
+        details[coord_offset+parameter_num] = 2**k
+    details[constants_offset] = 1   # fast_coord_count: one fast index
+    details[constants_offset+1] = theta_par
+    print ("details",details)
+    return details
+
+def constrained_poly_details(model_info, weights, constraints):
+    # Need to find the independently varying pars and sort them
+    # Need to build a coordination list for the dependent variables
+    # Need to generate a constraints function which takes values
+    # and weights, returning par blocks
+    raise NotImplementedError("Can't handle constraints yet")
+
+
+def create_default_functions(model_info):
+    """
+    Autogenerate missing functions, such as Iqxy from Iq.
+
+    This only works for Iqxy when Iq is written in python. :func:`make_source`
+    performs a similar role for Iq written in C.
+    """
+    if model_info['Iq'] is not None and model_info['Iqxy'] is None:
+        partable = model_info['parameters']
+        if partable.type['1d'] != partable.type['2d']:
+            raise ValueError("Iqxy model is missing")
+        Iq = model_info['Iq']
+        def Iqxy(qx, qy, **kw):
+            return Iq(np.sqrt(qx**2 + qy**2), **kw)
+        model_info['Iqxy'] = Iqxy
+
 
 def make_model_info(kernel_module):
@@ -678,9 +696,5 @@
       model variants (e.g., the list of cases) or is None if there are no
       model variants
-    * *defaults* is the *{parameter: value}* table built from the parameter
-      description table.
-    * *limits* is the *{parameter: [min, max]}* table built from the
-      parameter description table.
-    * *partypes* categorizes the model parameters. See
+    * *par_type* categorizes the model parameters. See
       :func:`categorize_parameters` for details.
     * *demo* contains the *{parameter: value}* map used in compare (and maybe
@@ -699,8 +713,12 @@
       *model_info* blocks for the composition objects.  This allows us to
       build complete product and mixture models from just the info.
+    * *max_pd* is the max polydispersity dimension.  This is constant and
+      should not be reset.  You may be able to change it when the program
+      starts by setting *sasmodels.generate.MAX_PD*.
 
     """
     # TODO: maybe turn model_info into a class ModelDefinition
-    parameters = COMMON_PARAMETERS + kernel_module.parameters
+    #print("make parameter table", kernel_module.parameters)
+    parameters = make_parameter_table(kernel_module.parameters)
     filename = abspath(kernel_module.__file__)
     kernel_id = splitext(basename(filename))[0]
@@ -731,4 +749,9 @@
     functions = "ER VR form_volume Iq Iqxy shape sesans".split()
     model_info.update((k, getattr(kernel_module, k, None)) for k in functions)
+    create_default_functions(model_info)
+    # Precalculate the monodisperse parameters
+    # TODO: make this a lazy evaluator
+    # make_model_info is called for every model on sasview startup
+    model_info['mono_details'] = mono_details(model_info)
     return model_info
 
Index: sasmodels/kernel_header.c
===================================================================
--- sasmodels/kernel_header.c	(revision 5cf3c33c6d8945e61877718f81a29df7bbd1478d)
+++ sasmodels/kernel_header.c	(revision 5cf3c33c6d8945e61877718f81a29df7bbd1478d)
@@ -0,0 +1,108 @@
+#ifdef __OPENCL_VERSION__
+# define USE_OPENCL
+#elif defined(_OPENMP)
+# define USE_OPENMP
+#endif
+
+// If opencl is not available, then we are compiling a C function
+// Note: if using a C++ compiler, then define kernel as extern "C"
+#ifndef USE_OPENCL
+#  ifdef __cplusplus
+      #include <cstdio>
+      #include <cmath>
+      using namespace std;
+      #if defined(_MSC_VER)
+         #include <limits>
+         #include <float.h>
+         #define kernel extern "C" __declspec( dllexport )
+         static double trunc(double x) { return x>=0?floor(x):-floor(-x); }
+         static double fmin(double x, double y) { return x>y ? y : x; }
+         static double fmax(double x, double y) { return x<y ? y : x; }
+         static double isnan(double x) { return _isnan(x); }
+         #define NAN (std::numeric_limits<double>::quiet_NaN()) // non-signalling NaN
+         static double cephes_expm1(double x) {
+            // Adapted from the cephes math library.
+            // Copyright 1984 - 1992 by Stephen L. Moshier
+            if (x != x || x == 0.0) {
+               return x; // NaN and +/- 0
+            } else if (x < -0.5 || x > 0.5) {
+               return exp(x) - 1.0;
+            } else {
+               const double xsq = x*x;
+               const double p = (((
+                  +1.2617719307481059087798E-4)*xsq
+                  +3.0299440770744196129956E-2)*xsq
+                  +9.9999999999999999991025E-1);
+               const double q = ((((
+                  +3.0019850513866445504159E-6)*xsq
+                  +2.5244834034968410419224E-3)*xsq
+                  +2.2726554820815502876593E-1)*xsq
+                  +2.0000000000000000000897E0);
+               double r = x * p;
+               r =  r / (q - r);
+               return r+r;
+             }
+         }
+         #define expm1 cephes_expm1
+         typedef __int32 int32_t
+     #else
+         #define kernel extern "C"
+         #include <cstdint>
+     #endif
+     static void SINCOS(double angle, double &svar, double &cvar) { svar=sin(angle); cvar=cos(angle); }
+#  else
+     #include <inttypes.h>  // C99 guarantees that int32_t types is here
+     #include <stdio.h>
+     #include <tgmath.h> // C99 type-generic math, so sin(float) => sinf
+     // MSVC doesn't support C99, so no need for dllexport on C99 branch
+     #define kernel
+     #define SINCOS(angle,svar,cvar) do {const double _t_=angle; svar=sin(_t_);cvar=cos(_t_);} while (0)
+#  endif
+#  define global
+#  define local
+#  define constant const
+// OpenCL powr(a,b) = C99 pow(a,b), b >= 0
+// OpenCL pown(a,b) = C99 pow(a,b), b integer
+#  define powr(a,b) pow(a,b)
+#  define pown(a,b) pow(a,b)
+#else
+   typedef int int32_t;
+#  if defined(USE_SINCOS)
+#    define SINCOS(angle,svar,cvar) svar=sincos(angle,&cvar)
+#  else
+#    define SINCOS(angle,svar,cvar) do {const double _t_=angle; svar=sin(_t_);cvar=cos(_t_);} while (0)
+#  endif
+#endif
+
+// Standard mathematical constants:
+//   M_E, M_LOG2E, M_LOG10E, M_LN2, M_LN10, M_PI, M_PI_2=pi/2, M_PI_4=pi/4,
+//   M_1_PI=1/pi, M_2_PI=2/pi, M_2_SQRTPI=2/sqrt(pi), SQRT2, SQRT1_2=sqrt(1/2)
+// OpenCL defines M_constant_F for float constants, and nothing if double
+// is not enabled on the card, which is why these constants may be missing
+#ifndef M_PI
+#  define M_PI 3.141592653589793
+#endif
+#ifndef M_PI_2
+#  define M_PI_2 1.570796326794897
+#endif
+#ifndef M_PI_4
+#  define M_PI_4 0.7853981633974483
+#endif
+#ifndef M_E
+#  define M_E 2.718281828459045091
+#endif
+
+// Non-standard function library
+// pi/180, used for converting between degrees and radians
+// 4/3 pi for computing sphere volumes
+// square and cube for computing squares and cubes
+#ifndef M_PI_180
+#  define M_PI_180 0.017453292519943295
+#endif
+#ifndef M_4PI_3
+#  define M_4PI_3 4.18879020478639
+#endif
+static double square(double x) { return x*x; }
+static double cube(double x) { return x*x*x; }
+static double sinc(double x) { return x==0 ? 1.0 : sin(x)/x; }
+
Index: sasmodels/kernel_iq.c
===================================================================
--- sasmodels/kernel_iq.c	(revision 69aa451cd7078d77c255be72d487793b00b34eee)
+++ sasmodels/kernel_iq.c	(revision 69aa451cd7078d77c255be72d487793b00b34eee)
@@ -0,0 +1,218 @@
+
+/*
+    ##########################################################
+    #                                                        #
+    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
+    #   !!                                              !!   #
+    #   !!  KEEP THIS CODE CONSISTENT WITH KERNELPY.PY  !!   #
+    #   !!                                              !!   #
+    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
+    #                                                        #
+    ##########################################################
+*/
+
+#ifndef _PAR_BLOCK_ // protected block so we can include this code twice.
+#define _PAR_BLOCK_
+
+typedef struct {
+    int32_t pd_par[MAX_PD];     // id of the nth polydispersity variable
+    int32_t pd_length[MAX_PD];  // length of the nth polydispersity weight vector
+    int32_t pd_offset[MAX_PD];  // offset of pd weights in the value & weight vector
+    int32_t pd_stride[MAX_PD];  // stride to move to the next index at this level
+    int32_t pd_isvol[MAX_PD];   // True if parameter is a volume weighting parameter
+    int32_t par_offset[NPARS];  // offset of par values in the value & weight vector
+    int32_t par_coord[NPARS];   // polydispersity coordination bitvector
+    int32_t fast_coord_pars[NPARS]; // ids of the fast coordination parameters
+    int32_t fast_coord_count;   // number of parameters coordinated with pd 1
+    int32_t theta_par;          // id of spherical correction variable
+} ProblemDetails;
+
+typedef struct {
+    PARAMETER_TABLE;
+} ParameterBlock;
+#endif
+
+
+kernel
+void KERNEL_NAME(
+    int32_t nq,                 // number of q values
+    const int32_t pd_start,     // where we are in the polydispersity loop
+    const int32_t pd_stop,      // where we are stopping in the polydispersity loop
+    global const ProblemDetails *problem,
+    global const double *weights,
+    global const double *values,
+    global const double *q, // nq q values, with padding to boundary
+    global double *result,  // nq+3 return values, again with padding
+    const double cutoff     // cutoff in the polydispersity weight product
+    )
+{
+  // Storage for the current parameter values.  These will be updated as we
+  // walk the polydispersity cube.
+  local ParameterBlock local_values;  // current parameter values
+  double *pvec = (double *)(&local_values);  // Alias named parameters with a vector
+
+  local int offset[NPARS];  // NPARS excludes scale/background
+
+#if 0 // defined(USE_SHORTCUT_OPTIMIZATION)
+  if (problem->pd_length[0] == 1) {
+    // Shouldn't need to copy!!
+
+    for (int k=0; k < NPARS; k++) {
+      pvec[k] = values[k+2];  // skip scale and background
+    }
+
+    const double volume = CALL_VOLUME(local_values);
+    #ifdef USE_OPENMP
+    #pragma omp parallel for
+    #endif
+    for (int i=0; i < nq; i++) {
+      const double scattering = CALL_IQ(q, i, local_values);
+      result[i] = values[0]*scattering/volume + values[1];
+    }
+    return;
+  }
+  printf("falling through\n");
+#endif
+
+  printf("Entering polydispersity\n");
+  // Since we are no longer looping over the entire polydispersity hypercube
+  // for each q, we need to track the normalization values for each q in a
+  // separate work vector.
+  double norm;   // contains sum over weights
+  double vol; // contains sum over volume
+  double norm_vol; // contains weights over volume
+
+  // Initialize the results to zero
+  if (pd_start == 0) {
+    norm_vol = 0.0;
+    norm = 0.0;
+    vol = 0.0;
+
+    #ifdef USE_OPENMP
+    #pragma omp parallel for
+    #endif
+    for (int i=0; i < nq; i++) {
+      result[i] = 0.0;
+    }
+  } else {
+    //Pulling values from previous segment
+    norm = result[nq];
+    vol = result[nq+1];
+    norm_vol = result[nq+2];
+  }
+
+  // Location in the polydispersity hypercube, one index per dimension.
+  local int pd_index[MAX_PD];
+
+  // Trigger the reset behaviour that happens at the end the fast loop
+  // by setting the initial index >= weight vector length.
+  pd_index[0] = problem->pd_length[0];
+
+
+  // need product of weights at every Iq calc, so keep product of
+  // weights from the outer loops so that weight = partial_weight * fast_weight
+  double partial_weight = NAN; // product of weight w4*w3*w2 but not w1
+  double partial_volweight = NAN;
+  double weight = 1.0;        // set to 1 in case there are no weights
+  double vol_weight = 1.0;    // set to 1 in case there are no vol weights
+  double spherical_correction = 1.0;  // correction for latitude variation
+
+  // Loop over the weights then loop over q, accumulating values
+  for (int loop_index=pd_start; loop_index < pd_stop; loop_index++) {
+    // check if indices need to be updated
+    if (pd_index[0] >= problem->pd_length[0]) {
+
+      // RESET INDICES
+      pd_index[0] = loop_index%problem->pd_length[0];
+      partial_weight = 1.0;
+      partial_volweight = 1.0;
+      for (int k=1; k < MAX_PD; k++) {
+        pd_index[k] = (loop_index%problem->pd_length[k])/problem->pd_stride[k];
+        const double wi = weights[problem->pd_offset[0]+pd_index[0]];
+        partial_weight *= wi;
+        if (problem->pd_isvol[k]) partial_volweight *= wi;
+      }
+      for (int k=0; k < NPARS; k++) {
+        int coord = problem->par_coord[k];
+        int this_offset = problem->par_offset[k];
+        int block_size = 1;
+        for (int bit=0; bit < MAX_PD && coord != 0; bit++) {
+          if (coord&1) {
+              this_offset += block_size * pd_index[bit];
+              block_size *= problem->pd_length[bit];
+          }
+          coord /= 2;
+        }
+        offset[k] = this_offset;
+        pvec[k] = values[this_offset];
+      }
+      weight = partial_weight * weights[problem->pd_offset[0]+pd_index[0]];
+      if (problem->theta_par >= 0) {
+        spherical_correction = fabs(cos(M_PI_180*pvec[problem->theta_par]));
+      }
+      if (problem->theta_par == problem->pd_par[0]) {
+        weight *= spherical_correction;
+      }
+
+    } else {
+
+      // INCREMENT INDICES
+      pd_index[0] += 1;
+      const double wi = weights[problem->pd_offset[0]+pd_index[0]];
+      weight = partial_weight*wi;
+      if (problem->pd_isvol[0]) vol_weight *= wi;
+      for (int k=0; k < problem->fast_coord_count; k++) {
+         printf("fast loop %d coord %d idx %d ?%d offset %d %g\n",
+         k,problem->fast_coord_pars[k],pd_index[0],
+        problem->fast_coord_pars[k],
+            offset[problem->fast_coord_pars[k]],
+            values[offset[problem->fast_coord_pars[k]]]
+          );
+        pvec[problem->fast_coord_pars[k]]
+            = values[offset[problem->fast_coord_pars[k]]++];
+
+      }
+      if (problem->theta_par ==problem->pd_par[0]) {
+        weight *= fabs(cos(M_PI_180*pvec[problem->theta_par]));
+      }
+    }
+    printf("rad len %f %f\n",local_values.radius, local_values.length);
+    #ifdef INVALID
+    if (INVALID(local_values)) continue;
+    #endif
+
+    // Accumulate I(q)
+    // Note: weight==0 must always be excluded
+    if (weight > cutoff) {
+      norm += weight;
+      vol += vol_weight * CALL_VOLUME(local_values);
+      norm_vol += vol_weight;
+
+      #ifdef USE_OPENMP
+      #pragma omp parallel for
+      #endif
+      for (int i=0; i < nq; i++) {
+        const double scattering = CALL_IQ(q, i, local_values);
+        result[i] += weight*scattering;
+      }
+    }
+  }
+
+  //Makes a normalization avialable for the next round
+  result[nq] = norm;
+  result[nq+1] = vol;
+  result[nq+2] = norm_vol;
+
+  //End of the PD loop we can normalize
+  if (pd_stop >= problem->pd_stride[MAX_PD-1]) {
+    #ifdef USE_OPENMP
+    #pragma omp parallel for
+    #endif
+    for (int i=0; i < nq; i++) {
+      if (vol*norm_vol != 0.0) {
+        result[i] *= norm_vol/vol;
+      }
+      result[i] = values[0]*result[i]/norm + values[1];
+    }
+  }
+}
Index: sasmodels/kernelcl.py
===================================================================
--- sasmodels/kernelcl.py	(revision 8e0d974c44c6bf60f95c3ad8c308c4d0b1399656)
+++ sasmodels/kernelcl.py	(revision fec69dd09a0e9893e1cb253c163d9efbd4d4db7a)
@@ -75,4 +75,19 @@
 
 
+# Pragmas for enable OpenCL features.  Be sure to protect them so that they
+# still compile even if OpenCL is not present.
+_F16_PRAGMA = """\
+#if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp16)
+#  pragma OPENCL EXTENSION cl_khr_fp16: enable
+#endif
+"""
+
+_F64_PRAGMA = """\
+#if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp64)
+#  pragma OPENCL EXTENSION cl_khr_fp64: enable
+#endif
+"""
+
+
 ENV = None
 def environment():
@@ -142,10 +157,17 @@
         raise RuntimeError("%s not supported for devices"%dtype)
 
-    source = generate.convert_type(source, dtype)
+    source_list = [generate.convert_type(source, dtype)]
+
+    if dtype == generate.F16:
+        source_list.insert(0, _F16_PRAGMA)
+    elif dtype == generate.F64:
+        source_list.insert(0, _F64_PRAGMA)
+
     # Note: USE_SINCOS makes the intel cpu slower under opencl
     if context.devices[0].type == cl.device_type.GPU:
-        source = "#define USE_SINCOS\n" + source
+        source_list.insert(0, "#define USE_SINCOS\n")
     options = (get_fast_inaccurate_build_options(context.devices[0])
                if fast else [])
+    source = "\n".join(source)
     program = cl.Program(context, source).build(options=options)
     return program
@@ -220,7 +242,8 @@
         key = "%s-%s-%s"%(name, dtype, fast)
         if key not in self.compiled:
-            #print("compiling",name)
+            print("compiling",name)
             dtype = np.dtype(dtype)
-            program = compile_model(self.get_context(dtype), source, dtype, fast)
+            program = compile_model(self.get_context(dtype),
+                                    str(source), dtype, fast)
             self.compiled[key] = program
         return self.compiled[key]
@@ -314,9 +337,9 @@
         self.program = None
 
-    def __call__(self, q_vectors):
+    def make_kernel(self, q_vectors):
         if self.program is None:
             compiler = environment().compile_program
-            self.program = compiler(self.info['name'], self.source, self.dtype,
-                                    self.fast)
+            self.program = compiler(self.info['name'], self.source,
+                                    self.dtype, self.fast)
         is_2d = len(q_vectors) == 2
         kernel_name = generate.kernel_name(self.info, is_2d)
@@ -365,12 +388,22 @@
         # at this point, so instead using 32, which is good on the set of
         # architectures tested so far.
-        self.q_vectors = [_stretch_input(q, self.dtype, 32) for q in q_vectors]
+        if self.is_2d:
+            # Note: 17 rather than 15 because results is 2 elements
+            # longer than input.
+            width = ((self.nq+17)//16)*16
+            self.q = np.empty((width, 2), dtype=dtype)
+            self.q[:self.nq, 0] = q_vectors[0]
+            self.q[:self.nq, 1] = q_vectors[1]
+        else:
+            # Note: 33 rather than 31 because results is 2 elements
+            # longer than input.
+            width = ((self.nq+33)//32)*32
+            self.q = np.empty(width, dtype=dtype)
+            self.q[:self.nq] = q_vectors[0]
+        self.global_size = [self.q.shape[0]]
         context = env.get_context(self.dtype)
-        self.global_size = [self.q_vectors[0].size]
         #print("creating inputs of size", self.global_size)
-        self.q_buffers = [
-            cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=q)
-            for q in self.q_vectors
-        ]
+        self.q_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR,
+                             hostbuf=self.q)
 
     def release(self):
@@ -378,7 +411,7 @@
         Free the memory.
         """
-        for b in self.q_buffers:
-            b.release()
-        self.q_buffers = []
+        if self.q is not None:
+            self.q.release()
+            self.q = None
 
     def __del__(self):
@@ -406,11 +439,14 @@
     """
     def __init__(self, kernel, model_info, q_vectors, dtype):
+        max_pd = self.info['max_pd']
+        npars = len(model_info['parameters'])-2
         q_input = GpuInput(q_vectors, dtype)
+        self.dtype = dtype
+        self.dim = '2d' if q_input.is_2d else '1d'
         self.kernel = kernel
         self.info = model_info
-        self.res = np.empty(q_input.nq, q_input.dtype)
-        dim = '2d' if q_input.is_2d else '1d'
-        self.fixed_pars = model_info['partype']['fixed-' + dim]
-        self.pd_pars = model_info['partype']['pd-' + dim]
+        self.pd_stop_index = 4*max_pd-1
+        # plus three for the normalization values
+        self.result = np.empty(q_input.nq+3, q_input.dtype)
 
         # Inputs and outputs for each kernel call
@@ -418,48 +454,40 @@
         env = environment()
         self.queue = env.get_queue(dtype)
-        self.loops_b = cl.Buffer(self.queue.context, mf.READ_WRITE,
-                                 2 * MAX_LOOPS * q_input.dtype.itemsize)
-        self.res_b = cl.Buffer(self.queue.context, mf.READ_WRITE,
+
+        # details is int32 data, padded to an 8 integer boundary
+        size = ((max_pd*5 + npars*3 + 2 + 7)//8)*8
+        self.result_b = cl.Buffer(self.queue.context, mf.READ_WRITE,
                                q_input.global_size[0] * q_input.dtype.itemsize)
-        self.q_input = q_input
-
-        self._need_release = [self.loops_b, self.res_b, self.q_input]
-
-    def __call__(self, fixed_pars, pd_pars, cutoff):
+        self.q_input = q_input # allocated by GpuInput above
+
+        self._need_release = [ self.result_b, self.q_input ]
+
+    def __call__(self, details, weights, values, cutoff):
         real = (np.float32 if self.q_input.dtype == generate.F32
                 else np.float64 if self.q_input.dtype == generate.F64
                 else np.float16 if self.q_input.dtype == generate.F16
                 else np.float32)  # will never get here, so use np.float32
-
-        #print "pars", fixed_pars, pd_pars
-        res_bi = self.res_b
-        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.q_input.dtype)
-            loops = np.ascontiguousarray(loops.T, self.q_input.dtype).flatten()
-            #print("loops",Nloops, loops)
-
-            #import sys; print("opencl eval",pars)
-            #print("opencl eval",pars)
-            if len(loops) > 2 * MAX_LOOPS:
-                raise ValueError("too many polydispersity points")
-
-            loops_bi = self.loops_b
-            cl.enqueue_copy(self.queue, loops_bi, loops)
-            loops_l = cl.LocalMemory(len(loops.data))
-            #ctx = environment().context
-            #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.q_input.q_buffers + [res_bi, nq] + dispersed + fixed
+        assert details.dtype == np.int32
+        assert weights.dtype == real and values.dtype == real
+
+        context = self.queue.context
+        details_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR,
+                              hostbuf=details)
+        weights_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR,
+                              hostbuf=weights)
+        values_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR,
+                             hostbuf=values)
+
+        start, stop = 0, self.details[self.pd_stop_index]
+        args = [
+            np.uint32(self.q_input.nq), np.uint32(start), np.uint32(stop),
+            self.details_b, self.weights_b, self.values_b,
+            self.q_input.q_b, self.result_b, real(cutoff),
+        ]
         self.kernel(self.queue, self.q_input.global_size, None, *args)
-        cl.enqueue_copy(self.queue, self.res, res_bi)
-
-        return self.res
+        cl.enqueue_copy(self.queue, self.result, self.result_b)
+        [v.release() for v in details_b, weights_b, values_b]
+
+        return self.result[:self.nq]
 
     def release(self):
Index: sasmodels/kerneldll.py
===================================================================
--- sasmodels/kerneldll.py	(revision 6ad0e875ccffda2d8d0cb50c477d8c365bf0ffcd)
+++ sasmodels/kerneldll.py	(revision 69aa451cd7078d77c255be72d487793b00b34eee)
@@ -50,6 +50,5 @@
 import tempfile
 import ctypes as ct
-from ctypes import c_void_p, c_int, c_longdouble, c_double, c_float
-import _ctypes
+from ctypes import c_void_p, c_int32, c_longdouble, c_double, c_float
 
 import numpy as np
@@ -81,5 +80,5 @@
         COMPILE = "gcc -shared -fPIC -std=c99 -O2 -Wall %(source)s -o %(output)s -lm"
         if "SAS_OPENMP" in os.environ:
-            COMPILE = COMPILE + " -fopenmp"
+            COMPILE += " -fopenmp"
 else:
     COMPILE = "cc -shared -fPIC -fopenmp -std=c99 -O2 -Wall %(source)s -o %(output)s -lm"
@@ -141,7 +140,6 @@
  
     source = generate.convert_type(source, dtype)
-    source_files = generate.model_sources(model_info) + [model_info['filename']]
+    newest = generate.timestamp(model_info)
     dll = dll_path(model_info, dtype)
-    newest = max(os.path.getmtime(f) for f in source_files)
     if not os.path.exists(dll) or os.path.getmtime(dll) < newest:
         # Replace with a proper temp file
@@ -171,8 +169,4 @@
     return DllModel(filename, model_info, dtype=dtype)
 
-
-IQ_ARGS = [c_void_p, c_void_p, c_int]
-IQXY_ARGS = [c_void_p, c_void_p, c_void_p, c_int]
-
 class DllModel(object):
     """
@@ -197,9 +191,4 @@
 
     def _load_dll(self):
-        Nfixed1d = len(self.info['partype']['fixed-1d'])
-        Nfixed2d = len(self.info['partype']['fixed-2d'])
-        Npd1d = len(self.info['partype']['pd-1d'])
-        Npd2d = len(self.info['partype']['pd-2d'])
-
         #print("dll", self.dllpath)
         try:
@@ -212,13 +201,11 @@
               else c_double if self.dtype == generate.F64
               else c_longdouble)
-        pd_args_1d = [c_void_p, fp] + [c_int]*Npd1d if Npd1d else []
-        pd_args_2d = [c_void_p, fp] + [c_int]*Npd2d if Npd2d else []
+
+        # int, int, int, int*, double*, double*, double*, double*, double*, double
+        argtypes = [c_int32]*3 + [c_void_p]*5 + [fp]
         self.Iq = self.dll[generate.kernel_name(self.info, False)]
-        self.Iq.argtypes = IQ_ARGS + pd_args_1d + [fp]*Nfixed1d
-
         self.Iqxy = self.dll[generate.kernel_name(self.info, True)]
-        self.Iqxy.argtypes = IQXY_ARGS + pd_args_2d + [fp]*Nfixed2d
-        
-        self.release()
+        self.Iq.argtypes = argtypes
+        self.Iqxy.argtypes = argtypes
 
     def __getstate__(self):
@@ -229,5 +216,5 @@
         self.dll = None
 
-    def __call__(self, q_vectors):
+    def make_kernel(self, q_vectors):
         q_input = PyInput(q_vectors, self.dtype)
         if self.dll is None: self._load_dll()
@@ -272,36 +259,34 @@
     """
     def __init__(self, kernel, model_info, q_input):
+        self.kernel = kernel
         self.info = model_info
         self.q_input = q_input
-        self.kernel = kernel
-        self.res = np.empty(q_input.nq, q_input.dtype)
-        dim = '2d' if q_input.is_2d else '1d'
-        self.fixed_pars = model_info['partype']['fixed-' + dim]
-        self.pd_pars = model_info['partype']['pd-' + dim]
-
-        # In dll kernel, but not in opencl kernel
-        self.p_res = self.res.ctypes.data
-
-    def __call__(self, fixed_pars, pd_pars, cutoff):
+        self.dtype = q_input.dtype
+        self.dim = '2d' if q_input.is_2d else '1d'
+        self.result = np.empty(q_input.nq+3, q_input.dtype)
+
+    def __call__(self, details, weights, values, cutoff):
         real = (np.float32 if self.q_input.dtype == generate.F32
                 else np.float64 if self.q_input.dtype == generate.F64
                 else np.float128)
-
-        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.q_input.dtype).flatten()
-            p_loops = loops.ctypes.data
-            dispersed = [p_loops, cutoff] + loops_N
-        else:
-            dispersed = []
-        fixed = [real(p) for p in fixed_pars]
-        args = self.q_input.q_pointers + [self.p_res, nq] + dispersed + fixed
-        #print(pars)
+        assert details.dtype == np.int32
+        assert weights.dtype == real and values.dtype == real
+
+        max_pd = self.info['max_pd']
+        start, stop = 0, details[4*max_pd-1]
+        print(details)
+        args = [
+            self.q_input.nq, # nq
+            start, # pd_start
+            stop, # pd_stop pd_stride[MAX_PD]
+            details.ctypes.data, # problem
+            weights.ctypes.data,  # weights
+            values.ctypes.data,  #pars
+            self.q_input.q.ctypes.data, #q
+            self.result.ctypes.data,   # results
+            real(cutoff), # cutoff
+            ]
         self.kernel(*args)
-
-        return self.res
+        return self.result[:-3]
 
     def release(self):
Index: sasmodels/kernelpy.py
===================================================================
--- sasmodels/kernelpy.py	(revision a84a0cac83cecece77ff644df4f2552ef1452b62)
+++ sasmodels/kernelpy.py	(revision 48fbd50a25db17a9c6bf7770aeb061514b3e8879)
@@ -53,6 +53,11 @@
         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 self.q_vectors]
+        if self.is_2d:
+            self.q = np.empty((self.nq, 2), dtype=dtype)
+            self.q[:, 0] = q_vectors[0]
+            self.q[:, 1] = q_vectors[1]
+        else:
+            self.q = np.empty(self.nq, dtype=dtype)
+            self.q[:self.nq] = q_vectors[0]
 
     def release(self):
@@ -60,5 +65,5 @@
         Free resources associated with the model inputs.
         """
-        self.q_vectors = []
+        self.q = None
 
 class PyKernel(object):
Index: sasmodels/mixture.py
===================================================================
--- sasmodels/mixture.py	(revision 72a081dea6caa8fa493057133a6232f6b118c89a)
+++ sasmodels/mixture.py	(revision 69aa451cd7078d77c255be72d487793b00b34eee)
@@ -14,5 +14,5 @@
 import numpy as np
 
-from .generate import process_parameters, COMMON_PARAMETERS, Parameter
+from .modelinfo import Parameter, ParameterTable
 
 SCALE=0
@@ -34,19 +34,19 @@
 
     # Build new parameter list
-    pars = COMMON_PARAMETERS + []
+    pars = []
     for k, part in enumerate(parts):
         # Parameter prefix per model, A_, B_, ...
+        # Note that prefix must also be applied to id and length_control
+        # to support vector parameters
         prefix = chr(ord('A')+k) + '_'
-        for p in part['parameters']:
-            # No background on the individual mixture elements
-            if p.name == 'background':
-                continue
-            # TODO: promote Parameter to a full class
-            # this code knows too much about the implementation!
-            p = list(p)
-            if p[0] == 'scale':  # allow model subtraction
-                p[3] = [-np.inf, np.inf]  # limits
-            p[0] = prefix+p[0]   # relabel parameters with A_, ...
+        pars.append(Parameter(prefix+'scale'))
+        for p in part['parameters'].kernel_pars:
+            p = copy(p)
+            p.name = prefix+p.name
+            p.id = prefix+p.id
+            if p.length_control is not None:
+                p.length_control = prefix+p.length_control
             pars.append(p)
+    partable = ParameterTable(pars)
 
     model_info = {}
@@ -58,5 +58,5 @@
     model_info['docs'] = model_info['title']
     model_info['category'] = "custom"
-    model_info['parameters'] = pars
+    model_info['parameters'] = partable
     #model_info['single'] = any(part['single'] for part in parts)
     model_info['structure_factor'] = False
@@ -67,6 +67,4 @@
     # Remember the component info blocks so we can build the model
     model_info['composition'] = ('mixture', parts)
-    process_parameters(model_info)
-    return model_info
 
 
Index: sasmodels/model_test.py
===================================================================
--- sasmodels/model_test.py	(revision a84a0cac83cecece77ff644df4f2552ef1452b62)
+++ sasmodels/model_test.py	(revision 48fbd50a25db17a9c6bf7770aeb061514b3e8879)
@@ -51,5 +51,5 @@
 
 from .core import list_models, load_model_info, build_model, HAVE_OPENCL
-from .core import make_kernel, call_kernel, call_ER, call_VR
+from .core import call_kernel, call_ER, call_VR
 from .exception import annotate_exception
 
@@ -187,9 +187,9 @@
                 Qx, Qy = zip(*x)
                 q_vectors = [np.array(Qx), np.array(Qy)]
-                kernel = make_kernel(model, q_vectors)
+                kernel = model.make_kernel(q_vectors)
                 actual = call_kernel(kernel, pars)
             else:
                 q_vectors = [np.array(x)]
-                kernel = make_kernel(model, q_vectors)
+                kernel = model.make_kernel(q_vectors)
                 actual = call_kernel(kernel, pars)
 
Index: sasmodels/modelinfo.py
===================================================================
--- sasmodels/modelinfo.py	(revision 69aa451cd7078d77c255be72d487793b00b34eee)
+++ sasmodels/modelinfo.py	(revision 69aa451cd7078d77c255be72d487793b00b34eee)
@@ -0,0 +1,325 @@
+
+import numpy as np
+
+# TODO: turn ModelInfo into a proper class
+ModelInfo = dict
+
+COMMON_PARAMETERS = [
+    ["scale", "", 1, [0, np.inf], "", "Source intensity"],
+    ["background", "1/cm", 1e-3, [0, np.inf], "", "Source background"],
+]
+assert (len(COMMON_PARAMETERS) == 2
+        and COMMON_PARAMETERS[0][0]=="scale"
+        and COMMON_PARAMETERS[1][0]=="background"), "don't change common parameters"
+# assumptions about common parameters exist throughout the code, such as:
+# (1) kernel functions Iq, Iqxy, form_volume, ... don't see them
+# (2) kernel drivers assume scale is par[0] and background is par[1]
+# (3) mixture models drop the background on components and replace the scale
+#     with a scale that varies from [-inf, inf]
+# (4) product models drop the background and reassign scale
+# and maybe other places.
+# Note that scale and background cannot be coordinated parameters whose value
+# depends on the some polydisperse parameter with the current implementation
+
+def make_parameter_table(pars):
+    processed = []
+    for p in pars:
+        if not isinstance(p, list) or len(p) != 6:
+            raise ValueError("Parameter should be [name, units, default, limits, type, desc], but got %r"
+                             %str(p))
+        processed.append(parse_parameter(*p))
+    partable = ParameterTable(processed)
+    set_vector_length_from_reference(partable)
+    return partable
+
+def set_vector_length_from_reference(partable):
+    # Sort out the length of the vector parameters such as thickness[n]
+    for p in partable:
+        if p.length_control:
+           ref = partable[p.length_control]
+           low, high = ref.limits
+           if int(low) != low or int(high) != high or low<0 or high>20:
+               raise ValueError("expected limits on %s to be within [0, 20]"%ref.name)
+           p.length = low
+
+def parse_parameter(name, units='', default=None,
+                    limits=(-np.inf, np.inf), type='', description=''):
+    # Parameter is a user facing class.  Do robust type checking.
+    if not isstr(name):
+        raise ValueError("expected string for parameter name %r"%name)
+    if not isstr(units):
+        raise ValueError("expected units to be a string for %s"%name)
+    # if limits is a list of strings, then this is a choice list
+    # field, and limits are 1 to length of string list
+    if isinstance(limits, list) and all(isstr(k) for k in limits):
+        choices = limits
+        limits = [1, len(choices)]
+    else:
+        choices = None
+    # TODO: maybe allow limits of None for (-inf, inf)
+    try:
+        low, high = limits
+        if not isinstance(low, (int, float)):
+            raise TypeError("low is not numeric")
+        if not isinstance(high, (int, float)):
+            raise TypeError("high is not numeric")
+        if low >= high:
+            raise ValueError("require low < high")
+    except:
+        raise ValueError("invalid limits %s for %s"%(limits, name))
+
+    if not isinstance(default, (int, float)):
+        raise ValueError("expected default %r to be a number for %s"
+                         % (default, name))
+    if default < low or default > high:
+        raise ValueError("default value %r not in range for %s"
+                         % (default, name))
+
+    if type not in ("volume", "orientation", "sld", "magnetic", ""):
+        raise ValueError("unexpected type %r for %s" % (type, name))
+
+    if not isstr(description):
+        raise ValueError("expected description to be a string")
+
+
+    # Parameter id for name[n] does not include [n]
+    if "[" in name:
+        if not name.endswith(']'):
+            raise ValueError("Expected name[len] for vector parameter %s"%name)
+        pid, ref = name[:-1].split('[', 1)
+        ref = ref.strip()
+    else:
+        pid, ref = name, None
+
+
+    # automatically identify sld types
+    if type=='' and (pid.startswith('sld') or pid.endswith('sld')):
+        type = 'sld'
+
+    # Check if using a vector definition, name[k], as the parameter name
+    if ref:
+        if ref == '':
+            raise ValueError("Need to specify vector length for %s"%name)
+        try:
+            length = int(ref)
+            control = None
+        except:
+            length = None
+            control = ref
+    else:
+        length = 1
+        control = None
+
+    # Build the parameter
+    parameter = Parameter(name=name, units=units, default=default,
+                          limits=limits, type=type, description=description)
+
+    # TODO: need better control over whether a parameter is polydisperse
+    parameter.polydisperse = type in ('orientation', 'volume')
+    parameter.relative_pd = type in ('volume')
+    parameter.choices = choices
+    parameter.length = length
+    parameter.length_control = control
+
+    return parameter
+
+class Parameter(object):
+    """
+    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.
+
+    *limits = [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*.  "sld" parameters
+    can automatically be promoted to magnetic parameters, each of which
+    will have a magnitude and a direction, which may be different from
+    other sld parameters.
+
+    *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.
+
+    Additional values can be set after the parameter is created:
+
+    *length* is the length of the field if it is a vector field
+    *length_control* is the parameter which sets the vector length
+    *polydisperse* is true if the parameter accepts a polydispersity
+    *relative_pd* is true if that polydispersity is relative
+
+    In the usual process these values are set by :func:`make_parameter_table`
+    and :func:`parse_parameter` therein.
+    """
+    def __init__(self, name, units='', default=None, limits=(-np.inf, np.inf),
+                 type='', description=''):
+        self.id = name
+        self.name = name
+        self.default = default
+        self.limits = limits
+        self.type = type
+        self.description = description
+        self.choices = None
+
+        # Length and length_control will be filled in by
+        # set_vector_length_from_reference(partable) once the complete
+        # parameter table is available.
+        self.length = 1
+        self.length_control = None
+
+        # TODO: need better control over whether a parameter is polydisperse
+        self.polydisperse = False
+        self.relative_pd = False
+
+    def as_definition(self):
+        """
+        Declare space for the variable in a parameter structure.
+
+        For example, the parameter thickness with length 3 will
+        return "double thickness[3];", with no spaces before and
+        no new line character afterward.
+        """
+        if self.length == 1:
+            return "double %s;"%self.id
+        else:
+            return "double %s[%d];"%(self.id, self.length)
+
+    def as_function_argument(self):
+        """
+        Declare the variable as a function argument.
+
+        For example, the parameter thickness with length 3 will
+        return "double *thickness", with no spaces before and
+        no comma afterward.
+        """
+        if self.length == 1:
+            return "double %s"%self.id
+        else:
+            return "double *%s"%self.id
+
+    def as_call_reference(self, prefix=""):
+        # Note: if the parameter is a struct type, then we will need to use
+        # &prefix+id.  For scalars and vectors we can just use prefix+id.
+        return prefix + self.id
+
+    def __str__(self):
+        return "<%s>"%self.name
+
+    def __repr__(self):
+        return "P<%s>"%self.name
+
+class ParameterTable(object):
+    # scale and background are implicit parameters
+    COMMON = [Parameter(*p) for p in COMMON_PARAMETERS]
+
+    def __init__(self, parameters):
+        self.parameters = self.COMMON + parameters
+        self._name_table= dict((p.name, p) for p in parameters)
+        self._categorize_parameters()
+
+    def __getitem__(self, k):
+        if isinstance(k, (int, slice)):
+            return self.parameters[k]
+        else:
+            return self._name_table[k]
+
+    def __contains__(self, key):
+        return key in self._name_table
+
+    def __iter__(self):
+        return iter(self.parameters)
+
+    def kernel_pars(self, ptype=None):
+        """
+        Return the parameters to the user kernel which match the given type.
+
+        Types include '1d' for Iq kernels, '2d' for Iqxy kernels and
+        'volume' for form_volume kernels.
+        """
+        # Assumes background and scale are the first two parameters
+        if ptype is None:
+            return self.parameters[2:]
+        else:
+            return [p for p in self.parameters[2:] if p in self.type[ptype]]
+
+    def _categorize_parameters(self):
+        """
+        Build parameter categories out of the the parameter definitions.
+
+        Returns a dictionary of categories.
+
+        Note: these categories are subject to change, depending on the needs of
+        the UI and the needs of the kernel calling function.
+
+        The categories are as follows:
+
+        * *volume* list of volume parameter names
+        * *orientation* list of orientation parameters
+        * *magnetic* list of magnetic parameters
+        * *sld* list of parameters that have no type info
+        * *other* list of parameters that have no type info
+
+        Each parameter is in one and only one category.
+        """
+        pars = self.parameters
+
+        par_type = {
+            'volume': [], 'orientation': [], 'magnetic': [], 'sld': [], 'other': [],
+        }
+        for p in self.parameters:
+            par_type[p.type if p.type else 'other'].append(p)
+        par_type['1d'] = [p for p in pars if p.type not in ('orientation', 'magnetic')]
+        par_type['2d'] = [p for p in pars if p.type != 'magnetic']
+        par_type['magnetic'] = [p for p in pars]
+        par_type['pd'] = [p for p in pars if p.polydisperse]
+        par_type['pd_relative'] = [p for p in pars if p.relative_pd]
+        self.type = par_type
+
+        # find index of theta (or whatever variable is used for spherical
+        # normalization during polydispersity...
+        if 'theta' in par_type['2d']:
+            # TODO: may be an off-by 2 bug due to background and scale
+            # TODO: is theta always the polar coordinate?
+            self.theta_par = [k for k,p in enumerate(pars) if p.name=='theta'][0]
+        else:
+            self.theta_par = -1
+
+    @property
+    def defaults(self):
+        return dict((p.name, p.default) for p in self.parameters)
+
+    @property
+    def num_pd(self):
+        """
+        Number of distributional parameters in the model (polydispersity in
+        shape dimensions and orientational distributions).
+        """
+        return len(self.type['pd'])
+
+    @property
+    def has_2d(self):
+        return self.type['orientation'] or self.type['magnetic']
+
+
+def isstr(x):
+    # TODO: 2-3 compatible tests for str, including unicode strings
+    return isinstance(x, str)
+
Index: sasmodels/models/cylinder.c
===================================================================
--- sasmodels/models/cylinder.c	(revision 26141cba42ea53a3ca3b54e60b579f5525dd583d)
+++ sasmodels/models/cylinder.c	(revision 03cac087353e872f40154bf57ef8161eb5367fd4)
@@ -3,4 +3,6 @@
 double Iqxy(double qx, double qy, double sld, double solvent_sld,
     double radius, double length, double theta, double phi);
+
+#define INVALID(v) (v.radius<0 || v.length<0)
 
 double form_volume(double radius, double length)
@@ -15,5 +17,4 @@
     double length)
 {
-    // TODO: return NaN if radius<0 or length<0?
     // precompute qr and qh to save time in the loop
     const double qr = q*radius;
@@ -47,5 +48,4 @@
     double phi)
 {
-    // TODO: return NaN if radius<0 or length<0?
     double sn, cn; // slots to hold sincos function output
 
Index: sasmodels/models/gel_fit.c
===================================================================
--- sasmodels/models/gel_fit.c	(revision 30b4ddfb88f4c0a299757def2ee628c2a3d50a42)
+++ sasmodels/models/gel_fit.c	(revision 03cac087353e872f40154bf57ef8161eb5367fd4)
@@ -1,19 +1,3 @@
-double form_volume(void);
-
-double Iq(double q,
-          double guinier_scale,
-          double lorentzian_scale,
-          double gyration_radius,
-          double fractal_exp,
-          double cor_length);
-
-double Iqxy(double qx, double qy,
-          double guinier_scale,
-          double lorentzian_scale,
-          double gyration_radius,
-          double fractal_exp,
-          double cor_length);
-
-static double _gel_fit_kernel(double q,
+static double Iq(double q,
           double guinier_scale,
           double lorentzian_scale,
@@ -24,5 +8,5 @@
     // Lorentzian Term
     ////////////////////////double a(x[i]*x[i]*zeta*zeta);
-    double lorentzian_term = q*q*cor_length*cor_length;
+    double lorentzian_term = square(q*cor_length);
     lorentzian_term = 1.0 + ((fractal_exp + 1.0)/3.0)*lorentzian_term;
     lorentzian_term = pow(lorentzian_term, (fractal_exp/2.0) );
@@ -30,5 +14,5 @@
     // Exponential Term
     ////////////////////////double d(x[i]*x[i]*rg*rg);
-    double exp_term = q*q*gyration_radius*gyration_radius;
+    double exp_term = square(q*gyration_radius);
     exp_term = exp(-1.0*(exp_term/3.0));
 
@@ -37,39 +21,2 @@
     return result;
 }
-double form_volume(void){
-    // Unused, so free to return garbage.
-    return NAN;
-}
-
-double Iq(double q,
-          double guinier_scale,
-          double lorentzian_scale,
-          double gyration_radius,
-          double fractal_exp,
-          double cor_length)
-{
-    return _gel_fit_kernel(q,
-                          guinier_scale,
-                          lorentzian_scale,
-                          gyration_radius,
-                          fractal_exp,
-                          cor_length);
-}
-
-// Iqxy is never called since no orientation or magnetic parameters.
-double Iqxy(double qx, double qy,
-          double guinier_scale,
-          double lorentzian_scale,
-          double gyration_radius,
-          double fractal_exp,
-          double cor_length)
-{
-    double q = sqrt(qx*qx + qy*qy);
-    return _gel_fit_kernel(q,
-                          guinier_scale,
-                          lorentzian_scale,
-                          gyration_radius,
-                          fractal_exp,
-                          cor_length);
-}
-
Index: sasmodels/models/rpa.py
===================================================================
--- sasmodels/models/rpa.py	(revision aad336c963572c598c67e78f4e39f6c4a6b9298a)
+++ sasmodels/models/rpa.py	(revision 69aa451cd7078d77c255be72d487793b00b34eee)
@@ -86,29 +86,11 @@
 #   ["name", "units", default, [lower, upper], "type","description"],
 parameters = [
-    ["case_num", CASES, 0, [0, 10], "", "Component organization"],
+    ["case_num", "", 1, CASES, "", "Component organization"],
 
-    ["Na", "", 1000.0, [1, inf], "", "Degree of polymerization"],
-    ["Phia", "", 0.25, [0, 1], "", "volume fraction"],
-    ["va", "mL/mol", 100.0, [0, inf], "", "specific volume"],
-    ["La", "fm", 10.0, [-inf, inf], "", "scattering length"],
-    ["ba", "Ang", 5.0, [0, inf], "", "segment length"],
-
-    ["Nb", "", 1000.0, [1, inf], "", "Degree of polymerization"],
-    ["Phib", "", 0.25, [0, 1], "", "volume fraction"],
-    ["vb", "mL/mol", 100.0, [0, inf], "", "specific volume"],
-    ["Lb", "fm", 10.0, [-inf, inf], "", "scattering length"],
-    ["bb", "Ang", 5.0, [0, inf], "", "segment length"],
-
-    ["Nc", "", 1000.0, [1, inf], "", "Degree of polymerization"],
-    ["Phic", "", 0.25, [0, 1], "", "volume fraction"],
-    ["vc", "mL/mol", 100.0, [0, inf], "", "specific volume"],
-    ["Lc", "fm", 10.0, [-inf, inf], "", "scattering length"],
-    ["bc", "Ang", 5.0, [0, inf], "", "segment length"],
-
-    ["Nd", "", 1000.0, [1, inf], "", "Degree of polymerization"],
-    ["Phid", "", 0.25, [0, 1], "", "volume fraction"],
-    ["vd", "mL/mol", 100.0, [0, inf], "", "specific volume"],
-    ["Ld", "fm", 10.0, [-inf, inf], "", "scattering length"],
-    ["bd", "Ang", 5.0, [0, inf], "", "segment length"],
+    ["N[4]", "", 1000.0, [1, inf], "", "Degree of polymerization"],
+    ["Phi[4]", "", 0.25, [0, 1], "", "volume fraction"],
+    ["v[4]", "mL/mol", 100.0, [0, inf], "", "specific volume"],
+    ["L[4]", "fm", 10.0, [-inf, inf], "", "scattering length"],
+    ["b[4]", "Ang", 5.0, [0, inf], "", "segment length"],
 
     ["Kab", "", -0.0004, [-inf, inf], "", "Interaction parameter"],
Index: sasmodels/product.py
===================================================================
--- sasmodels/product.py	(revision 8e0d974c44c6bf60f95c3ad8c308c4d0b1399656)
+++ sasmodels/product.py	(revision 8e0d974c44c6bf60f95c3ad8c308c4d0b1399656)
@@ -99,6 +99,6 @@
         # a parameter map.
         par_map = {}
-        p_info = p_kernel.info['partype']
-        s_info = s_kernel.info['partype']
+        p_info = p_kernel.info['par_type']
+        s_info = s_kernel.info['par_type']
         vol_pars = set(p_info['volume'])
         if dim == '2d':
Index: sasmodels/resolution.py
===================================================================
--- sasmodels/resolution.py	(revision 8b935d12a0041bac8658cc0f3268fc22ca46a544)
+++ sasmodels/resolution.py	(revision 48fbd50a25db17a9c6bf7770aeb061514b3e8879)
@@ -477,5 +477,5 @@
     """
     from sasmodels import core
-    kernel = core.make_kernel(form, [q])
+    kernel = form.make_kernel([q])
     theory = core.call_kernel(kernel, pars)
     kernel.release()
@@ -502,7 +502,7 @@
     from scipy.integrate import romberg
 
-    if any(k not in form.info['defaults'] for k in pars.keys()):
-        keys = set(form.info['defaults'].keys())
-        extra = set(pars.keys()) - keys
+    par_set = set([p.name for p in form.info['parameters']])
+    if any(k not in par_set for k in pars.keys()):
+        extra = set(pars.keys()) - par_set
         raise ValueError("bad parameters: [%s] not in [%s]"%
                          (", ".join(sorted(extra)), ", ".join(sorted(keys))))
@@ -556,9 +556,10 @@
     from scipy.integrate import romberg
 
-    if any(k not in form.info['defaults'] for k in pars.keys()):
-        keys = set(form.info['defaults'].keys())
-        extra = set(pars.keys()) - keys
+    par_set = set([p.name for p in form.info['parameters']])
+    if any(k not in par_set for k in pars.keys()):
+        extra = set(pars.keys()) - par_set
         raise ValueError("bad parameters: [%s] not in [%s]"%
-                         (", ".join(sorted(extra)), ", ".join(sorted(keys))))
+                         (", ".join(sorted(extra)),
+                          ", ".join(sorted(pars.keys()))))
 
     _fn = lambda q, q0, dq: eval_form(q, form, pars)*gaussian(q, q0, dq)
@@ -693,5 +694,5 @@
     def _eval_sphere(self, pars, resolution):
         from sasmodels import core
-        kernel = core.make_kernel(self.model, [resolution.q_calc])
+        kernel = self.model.make_kernel([resolution.q_calc])
         theory = core.call_kernel(kernel, pars)
         result = resolution.apply(theory)
@@ -1062,5 +1063,5 @@
     model = core.build_model(model_info)
 
-    kernel = core.make_kernel(model, [resolution.q_calc])
+    kernel = model.make_kernel([resolution.q_calc])
     theory = core.call_kernel(kernel, pars)
     Iq = resolution.apply(theory)
Index: sasmodels/src
===================================================================
--- sasmodels/src	(revision 9c490bcf3c2afca690a8982acef8d9ab4213aff7)
+++ sasmodels/src	(revision 9c490bcf3c2afca690a8982acef8d9ab4213aff7)
@@ -0,0 +1,339 @@
+// GENERATED CODE --- DO NOT EDIT ---
+// Code is produced by sasmodels.gen from sasmodels/models/MODEL.c
+
+#ifdef __OPENCL_VERSION__
+# define USE_OPENCL
+#endif
+
+#define USE_KAHAN_SUMMATION 0
+
+// If opencl is not available, then we are compiling a C function
+// Note: if using a C++ compiler, then define kernel as extern "C"
+#ifndef USE_OPENCL
+#  ifdef __cplusplus
+      #include <cstdio>
+      #include <cmath>
+      using namespace std;
+      #if defined(_MSC_VER)
+         #include <limits>
+         #include <float.h>
+         #define kernel extern "C" __declspec( dllexport )
+         inline double trunc(double x) { return x>=0?floor(x):-floor(-x); }
+         inline double fmin(double x, double y) { return x>y ? y : x; }
+         inline double fmax(double x, double y) { return x<y ? y : x; }
+         inline double isnan(double x) { return _isnan(x); }
+         #define NAN (std::numeric_limits<double>::quiet_NaN()) // non-signalling NaN
+         static double cephes_expm1(double x) {
+            // Adapted from the cephes math library.
+            // Copyright 1984 - 1992 by Stephen L. Moshier
+            if (x != x || x == 0.0) {
+               return x; // NaN and +/- 0
+            } else if (x < -0.5 || x > 0.5) {
+               return exp(x) - 1.0;
+            } else {
+               const double xsq = x*x;
+               const double p = (((
+                  +1.2617719307481059087798E-4)*xsq
+                  +3.0299440770744196129956E-2)*xsq
+                  +9.9999999999999999991025E-1);
+               const double q = ((((
+                  +3.0019850513866445504159E-6)*xsq
+                  +2.5244834034968410419224E-3)*xsq
+                  +2.2726554820815502876593E-1)*xsq
+                  +2.0000000000000000000897E0);
+               double r = x * p;
+               r =  r / (q - r);
+               return r+r;
+             }
+         }
+         #define expm1 cephes_expm1
+     #else
+         #define kernel extern "C"
+     #endif
+     inline void SINCOS(double angle, double &svar, double &cvar) { svar=sin(angle); cvar=cos(angle); }
+#  else
+     #include <stdio.h>
+     #include <tgmath.h> // C99 type-generic math, so sin(float) => sinf
+     // MSVC doesn't support C99, so no need for dllexport on C99 branch
+     #define kernel
+     #define SINCOS(angle,svar,cvar) do {const double _t_=angle; svar=sin(_t_);cvar=cos(_t_);} while (0)
+#  endif
+#  define global
+#  define local
+#  define constant const
+// OpenCL powr(a,b) = C99 pow(a,b), b >= 0
+// OpenCL pown(a,b) = C99 pow(a,b), b integer
+#  define powr(a,b) pow(a,b)
+#  define pown(a,b) pow(a,b)
+#else
+#  if defined(USE_SINCOS)
+#    define SINCOS(angle,svar,cvar) svar=sincos(angle,&cvar)
+#  else
+#    define SINCOS(angle,svar,cvar) do {const double _t_=angle; svar=sin(_t_);cvar=cos(_t_);} while (0)
+#  endif
+#endif
+
+// Standard mathematical constants:
+//   M_E, M_LOG2E, M_LOG10E, M_LN2, M_LN10, M_PI, M_PI_2=pi/2, M_PI_4=pi/4,
+//   M_1_PI=1/pi, M_2_PI=2/pi, M_2_SQRTPI=2/sqrt(pi), SQRT2, SQRT1_2=sqrt(1/2)
+// OpenCL defines M_constant_F for float constants, and nothing if double
+// is not enabled on the card, which is why these constants may be missing
+#ifndef M_PI
+#  define M_PI 3.141592653589793
+#endif
+#ifndef M_PI_2
+#  define M_PI_2 1.570796326794897
+#endif
+#ifndef M_PI_4
+#  define M_PI_4 0.7853981633974483
+#endif
+#ifndef M_E
+#  define M_E 2.718281828459045091
+#endif
+
+// Non-standard function library
+// pi/180, used for converting between degrees and radians
+// 4/3 pi for computing sphere volumes
+// square and cube for computing squares and cubes
+#ifndef M_PI_180
+#  define M_PI_180 0.017453292519943295
+#endif
+#ifndef M_4PI_3
+#  define M_4PI_3 4.18879020478639
+#endif
+//inline double square(double x) { return pow(x,2.0); }
+//inline double square(double x) { return pown(x,2); }
+inline double square(double x) { return x*x; }
+inline double cube(double x) { return x*x*x; }
+inline double sinc(double x) { return x==0 ? 1.0 : sin(x)/x; }
+
+
+#define VOLUME_PARAMETER_DECLARATIONS void
+#define IQ_KERNEL_NAME bessel_Iq
+#define IQ_PARAMETERS ignored
+#define IQ_FIXED_PARAMETER_DECLARATIONS const double scale, \
+    const double background, \
+    const double ignored
+#define IQ_PARAMETER_DECLARATIONS double ignored
+#define IQXY_KERNEL_NAME bessel_Iqxy
+#define IQXY_PARAMETERS ignored
+#define IQXY_FIXED_PARAMETER_DECLARATIONS const double scale, \
+    const double background, \
+    const double ignored
+#define IQXY_PARAMETER_DECLARATIONS double ignored
+
+/*
+The wrapper for gamma function from OpenCL and standard libraries
+The OpenCL gamma function fails misserably on values lower than 1.0
+while works fine on larger values.
+We use gamma definition Gamma(t + 1) = t * Gamma(t) to compute
+to function for values lower than 1.0. Namely Gamma(t) = 1/t * Gamma(t + 1)
+*/
+
+
+inline double sas_gamma( double x) { return tgamma(x+1)/x; }
+
+double form_volume(VOLUME_PARAMETER_DECLARATIONS);
+double form_volume(VOLUME_PARAMETER_DECLARATIONS) {
+    
+    
+}
+
+
+double Iq(double q, IQ_PARAMETER_DECLARATIONS);
+double Iq(double q, IQ_PARAMETER_DECLARATIONS) {
+    
+    return sas_gamma(q);
+    
+}
+
+
+double Iqxy(double qx, double qy, IQXY_PARAMETER_DECLARATIONS);
+double Iqxy(double qx, double qy, IQXY_PARAMETER_DECLARATIONS) {
+    
+    // never called since no orientation or magnetic parameters.
+    //return -1.0;
+    
+}
+
+
+/*
+    ##########################################################
+    #                                                        #
+    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
+    #   !!                                              !!   #
+    #   !!  KEEP THIS CODE CONSISTENT WITH KERNELPY.PY  !!   #
+    #   !!                                              !!   #
+    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
+    #                                                        #
+    ##########################################################
+*/
+
+#ifdef IQ_KERNEL_NAME
+kernel void IQ_KERNEL_NAME(
+    global const double *q,
+    global double *result,
+    const int Nq,
+#ifdef IQ_OPEN_LOOPS
+  #ifdef USE_OPENCL
+    global double *loops_g,
+  #endif
+    local double *loops,
+    const double cutoff,
+    IQ_DISPERSION_LENGTH_DECLARATIONS,
+#endif
+    IQ_FIXED_PARAMETER_DECLARATIONS
+    )
+{
+#ifdef USE_OPENCL
+  #ifdef IQ_OPEN_LOOPS
+  // copy loops info to local memory
+  event_t e = async_work_group_copy(loops, loops_g, (IQ_DISPERSION_LENGTH_SUM)*2, 0);
+  wait_group_events(1, &e);
+  #endif
+
+  int i = get_global_id(0);
+  if (i < Nq)
+#else
+  #pragma omp parallel for
+  for (int i=0; i < Nq; i++)
+#endif
+  {
+    const double qi = q[i];
+#ifdef IQ_OPEN_LOOPS
+    double ret=0.0, norm=0.0;
+  #ifdef VOLUME_PARAMETERS
+    double vol=0.0, norm_vol=0.0;
+  #endif
+    IQ_OPEN_LOOPS
+    //for (int radius_i=0; radius_i < Nradius; radius_i++) {
+    //  const double radius = loops[2*(radius_i)];
+    //  const double radius_w = loops[2*(radius_i)+1];
+
+    const double weight = IQ_WEIGHT_PRODUCT;
+    if (weight > cutoff) {
+      const double scattering = Iq(qi, IQ_PARAMETERS);
+      // allow kernels to exclude invalid regions by returning NaN
+      if (!isnan(scattering)) {
+        ret += weight*scattering;
+        norm += weight;
+      #ifdef VOLUME_PARAMETERS
+        const double vol_weight = VOLUME_WEIGHT_PRODUCT;
+        vol += vol_weight*form_volume(VOLUME_PARAMETERS);
+        norm_vol += vol_weight;
+      #endif
+      }
+    //else { printf("exclude qx,qy,I:%g,%g,%g\n",qi,scattering); }
+    }
+    IQ_CLOSE_LOOPS
+  #ifdef VOLUME_PARAMETERS
+    if (vol*norm_vol != 0.0) {
+      ret *= norm_vol/vol;
+    }
+  #endif
+    result[i] = scale*ret/norm+background;
+#else
+    result[i] = scale*Iq(qi, IQ_PARAMETERS) + background;
+#endif
+  }
+}
+#endif
+
+
+#ifdef IQXY_KERNEL_NAME
+kernel void IQXY_KERNEL_NAME(
+    global const double *qx,
+    global const double *qy,
+    global double *result,
+    const int Nq,
+#ifdef IQXY_OPEN_LOOPS
+  #ifdef USE_OPENCL
+    global double *loops_g,
+  #endif
+    local double *loops,
+    const double cutoff,
+    IQXY_DISPERSION_LENGTH_DECLARATIONS,
+#endif
+    IQXY_FIXED_PARAMETER_DECLARATIONS
+    )
+{
+#ifdef USE_OPENCL
+  #ifdef IQXY_OPEN_LOOPS
+  // copy loops info to local memory
+  event_t e = async_work_group_copy(loops, loops_g, (IQXY_DISPERSION_LENGTH_SUM)*2, 0);
+  wait_group_events(1, &e);
+  #endif
+
+  int i = get_global_id(0);
+  if (i < Nq)
+#else
+  #pragma omp parallel for
+  for (int i=0; i < Nq; i++)
+#endif
+  {
+    const double qxi = qx[i];
+    const double qyi = qy[i];
+    #if USE_KAHAN_SUMMATION
+    double accumulated_error = 0.0;
+    #endif
+#ifdef IQXY_OPEN_LOOPS
+    double ret=0.0, norm=0.0;
+    #ifdef VOLUME_PARAMETERS
+    double vol=0.0, norm_vol=0.0;
+    #endif
+    IQXY_OPEN_LOOPS
+    //for (int radius_i=0; radius_i < Nradius; radius_i++) {
+    //  const double radius = loops[2*(radius_i)];
+    //  const double radius_w = loops[2*(radius_i)+1];
+
+    const double weight = IQXY_WEIGHT_PRODUCT;
+    if (weight > cutoff) {
+
+      const double scattering = Iqxy(qxi, qyi, IQXY_PARAMETERS);
+      if (!isnan(scattering)) { // if scattering is bad, exclude it from sum
+      //if (scattering >= 0.0) { // scattering cannot be negative
+        // TODO: use correct angle for spherical correction
+        // Definition of theta and phi are probably reversed relative to the
+        // equation which gave rise to this correction, leading to an
+        // attenuation of the pattern as theta moves through pi/2.  Either
+        // reverse the meanings of phi and theta in the forms, or use phi
+        // rather than theta in this correction.  Current code uses cos(theta)
+        // so that values match those of sasview.
+      #if defined(IQXY_HAS_THETA) // && 0
+        const double spherical_correction
+          = (Ntheta>1 ? fabs(cos(M_PI_180*theta))*M_PI_2:1.0);
+        const double next = spherical_correction * weight * scattering;
+      #else
+        const double next = weight * scattering;
+      #endif
+      #if USE_KAHAN_SUMMATION
+        const double y = next - accumulated_error;
+        const double t = ret + y;
+        accumulated_error = (t - ret) - y;
+        ret = t;
+      #else
+        ret += next;
+      #endif
+        norm += weight;
+      #ifdef VOLUME_PARAMETERS
+        const double vol_weight = VOLUME_WEIGHT_PRODUCT;
+        vol += vol_weight*form_volume(VOLUME_PARAMETERS);
+        norm_vol += vol_weight;
+      #endif
+      }
+      //else { printf("exclude qx,qy,I:%g,%g,%g\n",qi,scattering); }
+    }
+    IQXY_CLOSE_LOOPS
+  #ifdef VOLUME_PARAMETERS
+    if (vol*norm_vol != 0.0) {
+      ret *= norm_vol/vol;
+    }
+  #endif
+    result[i] = scale*ret/norm+background;
+#else
+    result[i] = scale*Iqxy(qxi, qyi, IQXY_PARAMETERS) + background;
+#endif
+  }
+}
+#endif
+
